基于离散元的木薯种茎仿真参数标定方法研究

陈林涛，薛俊祥，牟向伟，马 旭，刘文杰，刘雨飞，向金山

(1.广西师范大学职业技术师范学院，广西桂林 541004； 2.华南农业大学工程学院，广东广州 510642；3.广东省水稻移栽机械装备工程技术研究中心，广东广州 510642)

研制木薯精密播种机需对切种、供种以及投种等多道工序进行结构设计优化[1-2]。现有木薯播种机采用预切种方式进行作业，利用切种机将种秆提前切成形状简单、长度约150 mm的种茎，由排种机构将种茎排至种沟[3-4]。众多研究表明，离散元法可应用于农业散体物料与机械装备相互作用的仿真分析[5-7]。应用离散元法研究木薯种茎离散元模型仿真参数，可加快木薯排种器研制周期，并改善其作业性能[8-10]。离散元模型仿真参数可分为材料本征参数和接触力学参数。材料本征参数大多可以通过标准化试验进行测定，但接触力学参数往往为微观参数，无法通过实测直接得到。因此，需对仿真参数进行重新标定[11-15]。

国内外学者在种子颗粒标定方面的研究主要集中在水稻、三七等植物上[14-15]，对木薯种茎离散元参数标定研究鲜有报道。但对于类似作物秸秆的离散元仿真开展过一系列研究。王宪良等对秸秆模型堆积角、恢复系数等参数进行了标定[16]。王云霞等通过在两段秸秆之间添加连接点的方法，建立了一种弹性离散元秸秆模型，模型赋予连接点弹性阻尼，使其可承受扭矩和弯力[17]。鹿芳媛等分别对番茄藤、茄子藤等藤茎类秸秆进行离散元模拟，分析不同参数刀具切割受力情况[18]。木薯种茎为一定长度的圆柱状木质茎秆，且茎秆表面具有凸起，相比于水稻、三七以及番茄藤等作物在质量和形态方面均具有较大差异，目前研究较少且无类似作物作仿真参数的参考，难以开展木薯机械化装备的离散仿真研究。

针对木薯机械化种植过程中切种与播种等环节离散元仿真缺乏准确模型等问题，笔者以桂热4号木薯为研究对象，通过离散元软件，采用Hertz-Mindlin 基本模型，结合堆积角测试物理试验，以及筛选显著性因素和响应曲面试验，确定碰撞恢复系数、静摩擦因数和滚动摩擦因数等基本接触参数；采用Hertz-Mindlin with bonding接触模型，利用种茎弯曲破坏试验和响应曲面确定法向与切向接触刚度、法向与切向应力等黏结参数[19-23]。

1 材料与方法

1.1 木薯种茎基本物理参数测定

1.1.1 基本参数测定 选用桂热4号(广西壮族自治区亚热带作物研究所提供)，随机选取100颗种茎(平均含水率为68.46%，百颗质量5 987.83 g)进行直径测量统计。经统计，种茎直径在29～31 mm间占比为59颗(图1)。仿真取种茎平均直径为 30 mm。通过体积和质量测出种茎的平均密度为 720 kg/m3。种茎泊松比为0.42，径向剪切模量 17.60 MPa，径向弹性模量50.01 MPa，平均径向压缩强度1.43 MPa、平均径向剪切强度为2.53 MPa[24]。

1.1.2 堆积角测定 木薯种茎堆积角通过圆筒提升法试验进行测定[25]。将木薯种茎预先切制成长度为150 mm的试件。将钢质无底圆筒(内径 400 mm×高600 mm)固定在万能试验机上，将圆筒以50 mm/s速度向上提升，使木薯种茎形成颗粒堆。测定试验设为10次，记录试验结果并取其均值，试验现场见图2。

为减小测量的人为误差，通过图像处理技术提取并测量试验结果。将10组照片导入MATLAB进行数据处理(图3)，计算木薯种茎堆积角，堆积角(θ0)平均值为27.7°。

1.1.3 切断破坏试验 通过万能试验机进行种茎切断破坏试验，观察切断过程测定切刀所受阻力(图4)。将木薯种茎试样分别固定在夹具上，对切刀进行加载，使其垂直向下运动切断种茎。计算剪切强度，试验重复10次。切刀材料为Q235钢(泊松比0.28，剪切模量8.20×104MPa)，切刀运行速度为1 mm/s。剪切强度如下[26]：

(1)

式中：τ为剪切强度，MPa；Fmax为最大切断力，N；A为横截面积，mm2。

通过种茎切断试验发现，切断过程中切刀受力-位移变化见图5，可知当刀具接触到种茎为位移零点时，随着切刀不断下压，切力受力逐渐增大；直到木薯种茎切断的临界点，切刀受力达到最大值，种茎此时被切断，切刀受力急速下降。

经过试验得出10组最大切断力，代入公式(1)，得出10组不同直径下木薯种茎相对应的剪切强度值(表1)。得出剪切强度平均值为0.19 MPa，反求出当木薯种茎的直径为30 mm时最大切断力为268.47 N。

表1 木薯种茎试样的剪切强度值

1.2 仿真模型的建立

1.2.1 接触模型选取 由于种茎颗粒表面没有黏附力，故颗粒与部件之间接触模型选择无滑动接触模型。

1.2.2 堆积角仿真模型 为真实模拟出种茎特性，采用多球聚合模型法建立木薯种茎的颗粒离散元模型[27-29]，用60个大小相异球体组合成种茎颗粒，木薯种茎主秆模型由25个大小相同的球体组成，木薯种茎主秆上的突起由35个大小不同的球体颗粒组成，Particles模块创建，建立的颗粒模型见图6。

参考相关藤类农业颗粒物料的离散元仿真参数[30-33]，得出与本研究相关参数见表2。在EDEM中建立颗粒工厂，采用动态生成种茎颗粒总量250颗。设定颗粒以1 000 mm/s初始速度下落，向上提升圆筒的速度设为50 mm/s，最终形成姿态稳定的颗粒堆(图7)。

表2 木薯种茎颗粒堆积角仿真模型参数

1.2.3 切断破坏仿真模型 由于Hertz-Mindlin整体颗粒模型难以适用于木薯切断模拟仿真。因此，本研究在Hertz-Mindlin模型基础上，选用Hertz-Mindlin with bonding模型，利用EDEM快速填充API进行颗粒排列填充(图8)，导出网格文件.msh。

接着导入颗粒工厂模板，将准备好的颗粒工厂模板API放到EDEM颗粒工厂导入路径，将前面导出的point.txt里的网格坐标数据复制并替换到颗粒工厂模板文件下的Block_Factory_Data.txt中。将颗粒材料参数、几何体材料参数、接触参数以及颗粒模型设置完成，导入使用UG建立好并导出的剪切模型，在Geometries下拉菜单的Plugin Factories下添加自定义颗粒工厂，选择准备好的API文件，设置完毕进入计算页面，设置好参数，颗粒会在0.01 s填充完成(图9)。

最后如图10所示为种茎离散元颗粒模型，模型由 15 000 个半径1 mm颗粒等间距黏结而成。随后删除圆柱体实体模型，将黏结颗粒漏出，设置刀具运动参数，运动速度设置1 mm/s。

黏结参数参考藤茎类仿真参数，取值范围见表3。

表3 仿真参数

1.3 种茎参数标定试验

1.3.1 针对堆积角仿真参数标定试验设计 (1)设计二因子试验，筛选对茎秆颗粒堆积角影响显著的参数；(2)设计最陡爬坡试验；(3)设计基本接触参数响应曲面试验，并基于以上试验结果，根据BBD设计原理，进行堆积角仿真试验。

1.3.2 针对种茎剪切仿真参数标定试验设计 根据CCD设计原理，利用表3参数设计响应曲面试验，开展种茎切断破坏仿真分析，黏结参数设定见表4。

表4 黏结模型参数编码

2 结果与分析

2.1 因子试验

以x1～x6作为因素，开展因子试验，试验设计及结果见表5。

表5 试验设计与结果

方差分析见表6，种茎间静摩擦系数x2、种茎间滚动摩擦系数x3和种茎与钢之间静摩擦系数x5对堆积角的影响较为显著。

表6 参数显著性分析

2.2 最陡爬坡试验

由于x2、x3和x5这3个显著性参数对堆积角效应均是正值，因此设计参数值逐增。将堆积角数值及相对误差进行记录，结果见表7。随着3个选定的试验因素数值增大，堆积角数值也逐渐增大，其中相对误差呈先减后增的趋势，3号试验中相对误差最小，所以后期试验的中心点将采用3号试验中的各参数值，2号、4号试验则作为低、高水平进行响应面试验。

表7 试验设计与结果

2.3 响应曲面试验

取种茎间静摩擦系数x2、种茎间滚动摩擦系数x3和种茎与钢之间静摩擦系数x53个基本接触参数，开展3因素3水平响应曲面试验，试验设计与结果见表8。

表8 试验设计与结果

进行拟合分析，模型显著(P=0.000 7)且失拟项不显著，方程为

(2)

表9 堆积角响应曲面二次全模型方差分析

在已知x2、x3和x5的情况下，利用响应曲面分析已生成的堆积角模型，得出颗粒堆积角，但在堆积角已知条件下，x2、x3和x5的具体参数是堆积角的等值面曲线(图11)，故采用软件将实际值与仿真值间相对误差作为响应值分析，方程为

(3)

二次多项式回归模型P值为0.006，除了种茎间静摩擦系数x2、种茎与钢的静摩擦系数x5的交互项以及种茎与钢之间的静摩擦系数x5的平方项P值大于0.05之外，其他因素的显著性较好。

通过DesignExpert 12.0对式(3)求误差最小极值点，得x2、x3和x5为0.20、0.09和0.68。将求解出的参数带入公式(2)，拟合值θc=28.56°，与θ0=27.7°误差为3.1%。建立堆积角模型开展仿真，种茎平均堆积角θs=27.29°，模拟结果与实际结果接近。

2.4 黏结参数标定试验

利用法向接触刚度x7、切向接触刚度x8、临界法向应力x9、临界切向应力x10等参数值，开展木薯种茎切断破坏仿真。试验设计及仿真结果见表10。

表10 试验设计与结果

对结果进行拟合分析。由方差分析(表11)可知，二次全模型P值小于0.01，且失拟项不显著；其中x7、x8及其二者平方项、交互项有显著影响。

表11 破坏力响应曲面二次全模型方差分析

剔除不显著项，得到如下回归方程

(4)

之后通过响应值取破坏力与实测力之间的误差ΔF，进行方差分析，结果见表12。

表12 破坏力误差响应曲面二次全模型方差分析

剔除x9、x10优化调整后建立破坏力误差ΔF与各显著性因素之间的二次多项式：

(5)

求最小值，得x7=2.5，x8=2.0。

将其带入公式(5)中得ΔF=0.912 N，相对误差为 0.3%；同理，将其带入公式(4)中得F=268.195 N，这时的相对误差为0.1%； 对计算所得参数进行仿真，得到最终平均破坏力为267.4 N，与实际试验破坏力之间的相对误差为0.40%，表明标定方法可靠、参数正确。

3 验证试验

为验证标定参数的准确性，模型参数按表13进行设置。采用侧壁坍塌法(图12-a)来模拟木薯种茎在种箱里的运动，物理试验和模拟试验各10次，记录休止角，将10次休止角进行对比。将实体模型、木薯种茎离散元模型和标定接触参数导入EDEM中进行仿真，仿真过程见图12-b。每组试验重复10次，取平均值。

表13 木薯种茎离散元仿真参数

对比各试验次数实测值和仿真值，可得10次休止角相对误差均小于2%，10次休止角相对误差平均值为0.69%，表明标定方法正确、标定的参数准确。

4 结论

木薯种茎间碰撞恢复系数、木薯种茎与钢之间的碰撞恢复系数及木薯种茎与钢间的滚动摩擦系数对堆积角影响较小，其取值确定为0.40、0.40和0.05；木薯种茎间静摩擦系数、木薯种茎间滚动摩擦系数、木薯种茎与钢间的静摩擦系数对堆积角的影响显著。木薯种茎间静摩擦系数为0.20，木薯种茎间滚动摩擦系数为0.09，木薯种茎与钢静摩擦系数为0.68。通过方差分析得出临界法向应力和临界切向应力对木薯种茎切断破坏力影响较小，取值分别确定为4.00和4.00；法向、切向接触刚度及其平方项、交互项有显著影响；法向接触刚度、切向接触刚度值分别为2.50和2.00；通过确定的黏结参数开展5次仿真试验，与实测值相对误差为0.40%。采用侧壁坍塌法模拟木薯种茎在种箱里运动，分别对比各试验次数的实测值和仿真值，10次休止角相对误差的平均值为0.69%，表明标定方法正确。本研究通过堆积角和弯曲破坏切断试验，确定了种茎切断离散元模型参数，为建立木薯种茎作物切断过程离散元仿真模型提供一种新方法，以期为精密播种器的研制提供技术参考。