胡麻茎秆离散元柔性模型建立与接触参数试验验证

史瑞杰 戴 飞 赵武云 张锋伟 石林榕 郭军海

(甘肃农业大学机电工程学院， 兰州 730070)

0 引言

胡麻(LinumusitatissimumL.)即油用亚麻，是我国西北和华北黄土高原旱作农业区重要的油料作物，胡麻茎秆含有大量纤维，是提取纤维素的主要原料，胡麻籽可以榨油[1-3]。戴飞等[4-5]研制的自走式胡麻联合收割机、丘陵山地胡麻联合收割机填补了国内胡麻联合收获机具的空白；文献[6-9]研制的胡麻割晒机、全喂入式胡麻脱粒机、脱粒物料清选机完善了胡麻分段收获的机械化作业过程，使得胡麻机械化收获程度进一步提高。

随着离散元法在农业工程学科的广泛应用，物料间互作过程的研究取得了一系列研究成果，其中物料离散元模型及接触参数对仿真模型可靠性有重要影响[10-12]。近年来，JIA等[13]采用离散元法分析稻秸之间、稻秸与农机部件之间的相互作用关系，结果表明各标定参数下的仿真结果与实测值不具有显著差异；田辛亮等[14]对黑土区玉米秸秆-土壤混料进行离散元模型建立与参数标定，通过仿真试验对最优参数组合进行对比验证，结果表明标定参数真实可靠；石林榕等[12,15]对玉米籽粒、胡麻籽粒进行离散元仿真参数标定与排种试验验证，结果表明标定结果与试验验证结果误差较小；郝建军等[16]利用三维扫描逆向建模技术与EDEM软件建立油葵籽粒离散元模型，通过物理试验与虚拟仿真试验对仿真参数进行标定，结果表明油葵籽粒模型和标定所得的离散元仿真参数具有可靠性；侯占峰等[17]通过物理试验和仿真试验相结合的方法进行冰草种子物理性状参数测定与离散元仿真参数标定，结果表明标定所得最优参数可用于冰草种子丸化包衣过程的离散元仿真试验；韩树杰等[18]以不同参数组合下的堆积角为响应值确定新疆果园散体厩肥离散元仿真参数，得到显著性参数最优值并进行物理试验验证，结果表明标定参数可靠。精准农业的深入实施和快速发展对作物机械化收获提出了更高的要求，离散元法为现代农业装备数字化设计提供了新手段[19-20]。同时离散元仿真参数标定多针对谷物籽粒，针对谷物茎秆的研究较少，而对胡麻茎秆的离散元模型建立与接触参数标定鲜见报道。

针对胡麻茎秆含有大量纤维素，茎秆韧性强，联合收获动态仿真过程缺乏参数等现象，本文以胡麻茎秆为研究对象，拟通过胡麻茎秆生物力学特性试验确定其离散元法bonding模型建模参数，并以胡麻茎秆各部本征参数与接触参数试验值为高低水平，通过Plackett-Burman试验和Central Composite试验确定胡麻茎秆之间、茎秆与收获装备之间的接触参数，通过胡麻茎秆剪切试验与堆积角试验验证模型可靠性。

1 胡麻茎秆本征参数测定及离散元柔性模型建立

1.1 本征参数测定

1.1.1几何尺寸及分布规律

以旱地密植陇亚14号胡麻茎秆为试验材料，将胡麻茎秆分为根部、中部、颈部，分别随机截取任意长度的茎秆作为试验材料，使用精度0.01 mm数显游标卡尺分别测其样品直径和长度。茎秆样品与直径分布如图1所示。其中胡麻根部茎秆平均直径为2.13 mm，长度为16.98 mm；中部茎秆平均直径1.96 mm，长度17.7 mm；颈部茎秆平均直径1.53 mm，长度18.12 mm。

图1 胡麻茎秆样品与直径分布Fig.1 Sample and diameter distribution of flax stem

1.1.2密度与泊松比

使用卤素含水率测量仪进行胡麻茎秆样品含水率测定，根部含水率为3.87%，中部含水率3.72%，颈部含水率3.69%。使用精度0.001 g的电子天平对胡麻茎秆样品多次称量后取平均值，同时采用比重瓶测试法测量胡麻茎秆样品体积，重复多次取平均值，通过计算得根部密度为410.9 kg/m3、中部密度为485.1 kg/m3、颈部密度为465.6 kg/m3，平均密度为453.87 kg/m3。

随机选取胡麻茎秆样品中小样10根，在甘肃农业大学机电工程学院力学实验室使用Stable Micro Systems质构仪进行胡麻茎秆生物力学特性试验。试验时沿胡麻茎秆样品横向方向施加压力至茎秆不再变形为止，加载速度为20 mm/min。胡麻茎秆横向变形量由试验机获得，纵向变形量由数显游标卡尺测得，胡麻茎秆泊松比计算公式为

(1)

式中ε——胡麻茎秆泊松比

e′——胡麻茎秆纵向变形量，mm

e——胡麻茎秆横向变形量，mm

W1——加载前胡麻茎秆纵向长度，mm

W2——加载后胡麻茎秆纵向长度，mm

L1——加载前胡麻茎秆横向长度，mm

L2——加载后胡麻茎秆横向长度，mm

经过10次试验后由式(1)计算得胡麻茎秆泊松比平均值根部为0.082，中部为0.085，颈部为0.089，平均值为0.085。

使用CMT2502型电子万能材料试验机对胡麻茎秆样品进行多次拉伸试验，该机最大载荷为500 N，精度为±0.001 N，速度为1～500 mm/min，由计算机控制，试验结果取平均值。结果表明：拉伸弹性模量根部为2 848.23 MPa，中部为6 416.67 MPa，颈部为4 674.69 MPa[21]。胡麻茎秆剪切模量计算公式为

(2)

式中G——胡麻茎秆剪切模量，Pa

E——胡麻茎秆弹性模量，Pa

由式(2)计算得胡麻茎秆根部剪切模量为1 316.19 MPa，中部为2 956.99 MPa，颈部为2 146.32 MPa，平均值为2 139.83 MPa。

1.2 离散元柔性模型建立

采用离散元软件EDEM中Hertz-Mindlin with bonding模型建立胡麻茎秆柔性模型，在该模型中表征胡麻茎秆的相邻两颗粒在接触点处发生平行粘结，形成粘结键。粘结键的作用效果相当于胡麻茎秆内部作用力，胡麻茎秆受外力作用时通过粘结键破裂情况反映其相关力学特性。为保证能够获得接近于真实状态的胡麻茎秆柔性模型，本文采用双峰分布颗粒群建模方法。该方式颗粒半径服从正态分布，通过大颗粒占据空间位置，小颗粒提高填充密度，降低模型孔隙率，粘结力更加牢固，使颗粒群力学特性与实际情况更接近，且可以减小计算机仿真负荷[22-24]。采用Hertz-Mindlin with bonding模型时颗粒间粘结键主要参数有法向刚度Kn、切向刚度Ks、法向临界应力σ、切向临界应力γ和粘结半径Rj[25-26]，计算公式为

(3)

式中εa、εb——颗粒泊松比

Ea、Eb——颗粒弹性模量，MPa

ra、rb——颗粒半径，mm

F——临界压力，N

R——压缩面半径，mm

c——茎秆内聚力，MPa

φ——内摩擦角，(°)

根据1.1.2节及课题组前期胡麻茎秆力学压缩、剪切试验结果[21]，取c=4 MPa，φ=34°，将试验数据代入式(3)，得出胡麻根部茎秆颗粒之间Kn=1.2×109N/m3，Ks=6×108N/m3，σ=13.4 MPa，γ=13.04 MPa，Rj=0.25 mm；胡麻中部茎秆颗粒之间Kn=1.3×109N/m3，Ks=6.5×108N/m3，σ=5.6 MPa，γ=7.78 MPa，Rj=0.25 mm；胡麻颈部茎秆颗粒之间Kn=9×108N/m3，Ks=4.5×108N/m3，σ=1 MPa，γ=4.67 MPa，Rj=0.25 mm。将胡麻茎秆各部位参数取平均值，得出整株胡麻茎秆柔性模型参数，即Kn=1.13×109N/m3，Ks=5.6×108N/m3，σ=6.67 MPa，γ=8.5 MPa，Rj=0.25 mm；由此建立整株胡麻茎秆离散元柔性模型，其中根部离散元模型共产生粘结键203 095个，中部178 330个，颈部107 258个，整株270 060个，粘结效果较好。如表1所示。

表1 胡麻茎秆离散元bonding模型Tab.1 Discrete element bonding model of flax stem

2 离散元接触参数测定

胡麻联合收获仿真计算需要测定胡麻茎秆-接触部件、胡麻茎秆-胡麻茎秆间的动摩擦因数、静摩擦因数、恢复系数等，其中胡麻茎秆-接触部件采用自制动静摩擦因数测量仪进行试验得出，胡麻茎秆-胡麻茎秆间的动静摩擦因数采用圆筒提升法试验得出。

2.1 动摩擦因数与静摩擦因数

图2 胡麻根部茎秆动静摩擦因数试验Fig.2 Dynamical and static friction coefficient test of flax root stem1.机架 2.胡麻茎秆样品 3、8.支撑板 4.钢板 5.传感器 6.数显量角器 7.旋转轴 9.根部茎秆

胡麻茎秆-接触部件动静摩擦因数测量时采用自制动静摩擦因数测量仪(图2)，该仪器包括机架、支撑杆、支撑板、传感器、数显量角器、电动机、旋转轴等。工作时由电动机驱动支撑杆将支撑板升起，绕旋转轴做圆周运动，当支撑板达到一定角度时物料开始下滑，当物料滑至传感器中间位置时，传感器向电动机发出信号，电动机停止工作，支撑板停止旋转，数显量角器显示当前角度。测量动摩擦因数时物料开始发生滚动趋势时支撑板停止旋转，开始读数；测量静摩擦因数时物料滑至传感器时支撑板自动停止旋转，开始读数。

利用斜面倾角等于斜面上物体的静摩擦角这一原理可得到研究对象的静摩擦因数[27]，计算公式为

μ=tanθ

(4)

式中μ——胡麻茎秆-钢制部件间静摩擦因数

θ——支撑板停止工作时与水平面夹角，(°)

多次测量后根据式(4)计算动静摩擦因数后取平均值，可得到胡麻茎秆-钢制部件动静摩擦因数如表2所示。

表2 胡麻茎秆-钢制部件动静摩擦因数Tab.2 Dynamical and static friction coefficient of flax culm-steel parts

2.2 恢复系数

恢复系数采用自由落体试验确定[15]，试验时胡麻茎秆由350 mm高处自由下落至钢板表面，使用摄像机记录胡麻茎秆弹起高度(图3)，恢复系数f计算公式为

(5)

式中v1——胡麻茎秆碰撞前速度，m/s

v2——钢板碰撞前速度，m/s，碰撞前钢板静止不动，v2=0

v′1——胡麻茎秆碰撞后速度，m/s

v′2——钢板碰撞后速度，m/s，碰撞后钢板静止不动，v′2=0

h——碰撞前胡麻茎秆高度，mm

h′——碰撞后胡麻茎秆弹起高度，mm

图3 胡麻根部茎秆恢复系数试验Fig.3 Test of recovery coefficient of flax root stem1.摄像机 2、7.钢板 3.支架 4.白纸 5.直尺 6.根部茎秆

多次试验后取平均值得到胡麻茎秆不同部位与钢制接触部件间恢复系数根部为0.334，中部为0.335，颈部为0.302。

2.3 堆积角试验

堆积角试验可直观表达散体物料颗粒间的摩擦作用、流动特性等，目前学者针对不同物料采用的方法有注入法、倾斜法、圆筒提升法等[16-18]。依据胡麻茎秆流动特性，选择圆筒提升法进行胡麻茎秆堆积角试验，该试验在甘肃农业大学机电工程学院力学实验室进行，试验装置由CMT2502型电子万能材料试验机、主机、显示器、文具夹、圆筒和托盘组成，如图4所示。试验时设定试验机向上提升速度为20 mm/min，圆筒内径为30 mm,圆筒和试验机传感器间由文具夹刚性连接，圆筒和托盘均为不锈钢材质，圆筒由试验机提升后胡麻茎秆落入托盘中央，形成堆积角。每组试验重复进行10次后由摄像机拍摄可得胡麻茎秆各部位堆积角图像。

图4 堆积角试验装置与根部堆积角Fig.4 Stacking angle test device and flax roots stacking angle1.托盘 2.圆筒 3.文具夹 4.试验机 5.操作按钮 6.显示器 7.主机

堆积角测量试验结束后将试验图像优化出堆积角区域，并使用Python计算机编程语言进行灰度、二值化、边界轮廓提取等工作，可得到胡麻茎秆样品堆积角试验轮廓图(图5)，再使用Origin数据处理软件进行轮廓线坐标数据拟合，得到各部位堆积角数值为根部36.02°，中部38.78°，颈部40.93°。

图5 胡麻茎秆根部堆积角试验图像处理过程Fig.5 Processing process of flax root stem accumulation angle experiment pictures

3 离散元接触参数优化

胡麻茎秆作为整体建立离散元模型时只需1组参数，为保证胡麻茎秆整株建模真实性，本文由前期胡麻茎秆各部位接触参数研究拓展至整株接触参数研究，以胡麻茎秆各部位接触参数的区间值作为各试验的水平值进行试验与仿真研究。

3.1 Plackett-Burman爬坡试验

基于Plackett-Burman的最陡爬坡试验通过比较试验因素两个高低水平的差异与整体的差异筛选出各因子的显著性[28]。鉴于胡麻茎秆Hertz-Mindlin with bonding模型颗粒数量较多，仿真过程计算机负荷较大，本文通过颗粒堆积方法建立胡麻茎秆离散元模型，采用Plackett-Burman试验方法，以胡麻茎秆堆积角α为响应值，以胡麻茎秆各部位与自身间、接触元件间的离散元接触参数为试验因素，以各因素的试验值和文献[2-3,15,29]试验值作为水平值，试验因素编码见表3。Plackett-Burman试验时设置3个中心点，共进行15组试验，其中第6、12、13组为零水平组，试验方案和结果见表4，表中X1～X11为编码值。试验结果方差分析见表5。

表3 Plackett-Burman试验因素编码Tab.3 Plackett-Burman test factors and coding

根据表4的Plackett-Burman试验结果，对堆积角进行回归模型显著性分析，结果如表5所示，得到α回归模型为

α=38.56-0.087X1+0.075X2+0.083X3+0.56X4+ 1.57X5+1.00X6-0.19X7+0.67X8+0.67X9

(6)

表4 Plackett-Burman试验方案与结果Tab.4 Design and results of Plackett-Burman test scheme

表5 Plackett-Burman试验结果显著性分析Tab.5 Analysis of significance of parameters in Plackett-Burman test

由表5可知，该模型P<0.01、R2=0.99，表明主要效应模型显著，说明模型所拟合的回归方程与实际情况相符合，能够表示因素X1～X9对响应值堆积角的影响程度，可知X4、X5、X6、X8、X9极显著，其他项不显著。根据模型回归方程式(6)一次项系数得到各因素对堆积角的影响主次顺序为X5、X6、X8、X9、X4、X3、X2、X1、X7，结合显著性分析贡献率，X5、X6对响应值影响较大，其中X5的贡献率为51.33%，X6的贡献率为20.84%，X8、X9的贡献率较为接近，但和其他因素贡献率同样均小于10%。选取对堆积角影响程度和贡献率较大的因素X5、X6进行最陡爬坡试验和Central Composite试验，因素X5、X6依据标准化效应值将水平值依次增加，堆积角以根部、颈部、中部试验平均值38.57°为目标值，因素X8、X9值采用Plackett-Burman试验优化解0.28、0.14，其余因素水平值均采用该试验零水平值。根据最陡爬坡试验堆积角δ结果与真实试验堆积角β结果的相对误差η作为最终Central Composite试验最优解的水平值。η计算公式为

(7)

由表6可知，当X5、X6逐渐增大时，胡麻茎秆堆积角逐渐增大，相对误差先减小、后增大，在第4组试验时相对误差最小，为精确得出胡麻茎秆间接触参数，以第4组最陡爬坡试验参数为中心值，以第3、5组试验数据为水平值进行Central Composite试验，以堆积角为响应值，寻求最优解。

表6 最陡爬坡试验方案与结果Tab.6 Test scheme and results of the steepest climb

3.2 Central-Composite试验

根据最陡爬坡试验结果进行胡麻茎秆接触参数Central-Composite试验，得出Central-Composite试验堆积角ζ和相对误差η，探究Central-Composite试验堆积角ζ对胡麻茎秆-胡麻茎秆静摩擦因数X5、胡麻茎秆-胡麻茎秆滚动摩擦因数X6间的响应面影响效果。本次仿真试验中，其他接触参数按照最陡爬坡试验参数设置，设计方案及结果如表7、8所示。

表7 Central-Composite试验因素编码Tab.7 Central-Composite test factors and coding

试验结束后对表8试验结果进行二元回归分析，建立堆积角ζ与胡麻茎秆-胡麻茎秆静摩擦因数X5、胡麻茎秆-胡麻茎秆滚动摩擦因数X6间的二元回归模型

表8 Central-Composite试验方案与结果Tab.8 Test design scheme and results of Central-Composite

(8)

3.3 响应面寻优试验

对Central-Composite试验结果堆积角ζ进行显著性分析，结果如表9所示。

表9 响应面寻优试验显著性分析Tab.9 Significance analysis of response surface optimization test

图6 胡麻茎秆堆积角响应面Fig.6 Response surface of stem stacking angle of flax

根据Central-Composite结果及二次回归方程，以试验所得堆积角相对误差η最小为目标，对试验因子X5、X6进行最优解分析，设定目标函数及约束条件为

(9)

最终得到胡麻茎秆-胡麻茎秆静摩擦因数为0.508、胡麻茎秆-胡麻茎秆滚动摩擦因数为0.033，以此为最佳参数进行胡麻茎秆接触参数验证试验。

4 试验验证

验证试验包括剪切试验和堆积角试验，其中剪切试验对胡麻茎秆离散元bonding模型进行验证，堆积角试验对胡麻茎秆接触参数进行验证。

4.1 剪切试验

胡麻茎秆生物力学特性剪切试验采用2.3节仪器设备，使用试验机上端夹具水平夹持刀片垂直切割剪切试验台样品。试验机上端夹具运行速度为20 mm/min，多次试验后得到胡麻茎秆不同部位剪切最大载荷。同时，在EDEM软件中将胡麻茎秆bonding模型水平放置在支撑平面上，在bonding模型上方建立竖直几何平面模拟刀片切割胡麻茎秆样品，设置竖直向下方向运动速度为20 mm/min，试验过程如图7所示。

图7 胡麻根部茎秆剪切试验Fig.7 Test of flax root stem shear

4.2 堆积角试验

将1.1节胡麻茎秆样品混合均匀进行混杂状态下胡麻茎秆堆积角试验，试验过程采取3.1节方法，试验结束后采用图像处理办法获取胡麻茎秆堆积角数值，多次试验后取平均值。在EDEM软件中采用3.1节方法将茎秆模型简化后进行仿真试验，多次测量后同样采用图像处理方法获取堆积角仿真值，并与试验值对比分析，试验过程如图8所示。

图8 胡麻茎秆堆积角试验Fig.8 Test of stalk stacking angle of flax stem

4.3 试验结果分析

剪切验证试验结束后，在万能材料试验机中导出茎秆剪切最大载荷，在EDEM软件后处理中导出剪切平面所受压力；堆积角验证试验结束后将得到的图像采用图像处理办法获取胡麻茎秆堆积角。多次试验后取平均值，试验结果如表10所示。试验结果表明，胡麻茎秆离散元bonding模型与实际物理模型较为接近，误差较小，其中万能材料试验机得到根部、中部、颈部剪切最大载荷分别为45.97、29.42、19.64 N，离散元仿真得到根部、中部、颈部剪切最大载荷分别为46.74、30.33、20.71 N，相对误差分别为1.67%、3.09%、5.44%，最大剪切载荷差值、相对误差由根部向颈部递增，说明胡麻茎秆物理特性由根部向颈部变化明显，茎秆差异性显著。堆积角试验得出胡麻茎秆平均堆积角为38.58°，仿真试验得出胡麻茎秆堆积角为38.7°，差值为0.12°，相对误差为0.31%，差值较小。可见胡麻茎秆离散元柔性模型与接触参数和实际情况较为相符，可表征胡麻茎秆物理特性，并为离散元仿真提供柔性模型及其接触参数。

表10 胡麻茎秆离散元模型与接触参数验证试验结果Tab.10 Discrete element model and contact parameter verification test results of flax stem

5 结论

(1)以胡麻根部茎秆、中部茎秆、颈部茎秆为研究对象，以其本征参数为研究基础，计算得胡麻茎秆各部建模参数，采用离散元法bonding模型建模方法构建胡麻茎秆柔性模型，并以胡麻茎秆各部分本征参数与接触参数试验值为高低水平，通过Plackett-Burman试验和Central-Composite试验确定胡麻茎秆之间、茎秆与收获装备之间的接触参数，通过胡麻茎秆剪切试验与堆积角试验验证模型可靠性。

(2)结果表明，胡麻植株离散元柔性模型参数中法向刚度Kn为1.13×109N/m3，切向刚度Ks为5.6×108N/m3，法向临界应力σ为6.67 MPa，切向临界应力γ为8.5 MPa，粘结半径Rj为0.25 mm；胡麻茎秆-钢制部件间恢复系数、静摩擦因数、滚动摩擦因数最优值分别为0.33、0.28、0.14，胡麻茎秆-胡麻茎秆间恢复系数、静摩擦因数、滚动摩擦因数最优值分别为0.3、0.508、0.033。

(3)剪切与堆积角验证试验结果表明，剪切试验中根部、中部、颈部剪切最大载荷分别为45.97、29.42、19.64 N，离散元仿真中根部、中部、颈部剪切最大载荷分别为46.74、30.33、20.71 N，相对误差分别为1.67%、3.09%、5.44%；堆积角试验中胡麻茎秆平均堆积角为38.58°，仿真试验中胡麻茎秆堆积角为38.7°，差值为0.12°，相对误差为0.31%，差值较小，可为胡麻茎秆离散元仿真提供参考。