蓝亚麻蒴果离散元模型仿真参数标定

徐福龙，贺俊林，王月华，李姣姣，刘少华，吴楠

(山西农业大学农业工程学院，山西 太谷 030801)

亚麻一般可分为纤维用亚麻、油用亚麻和油纤兼用亚麻三种类型，亚麻因其优良的营养和医用、饲用价值,已经被应用于食品、医药、饲料等较多领域，有着良好的开发前景[1]，蓝亚麻别名宿根亚麻(LinumperenneL.var.sibiricumPlanch.)，为亚麻科(Linaceae)亚麻属多年生草本花卉[2-3].

目前对亚麻(胡麻)相关参数测定和标定的研究主要有石林榕等[4]进行的胡麻籽粒离散元仿真参数标定与排种试验验证，对胡麻籽粒的基本物理参数(三轴尺寸、质量密度、体积密度、泊松比、千粒质量、含水率、弹性模量)和接触力学参数(碰撞恢复系数、静摩擦系数和动摩擦系数)进行了标定；史瑞杰等[5]进行的胡麻茎秆生物力学特性试验，对胡麻茎杆的含水率、平均直径、抗拉强度、最大拉伸力、抗压强度、抗弯强度、剪切强度和剪切最大载荷进行了试验测定.而针对亚麻(胡麻)蒴果的仿真参数测定和标定还未见相关的研究报道.

基于离散单元法的数值模拟仿真软件EDEM已经在农业工程领域中得到了广泛应用[6-8].亚麻(胡麻)收获脱粒过程属于颗粒流体运动问题，可通过离散单元法对亚麻(胡麻)收获脱粒装置进行脱粒性能优化.

本试验通过物理试验法、查阅文献法、仿真优化设计法和计算法等方法获得蓝亚麻蒴果的千粒质量、含水率、密度、三轴尺寸、泊松比、弹性模量、剪切模量、蒴果之间及蒴果与接触材料之间的碰撞恢复系数、静摩擦系数和滚动摩擦系数等，对于亚麻(胡麻)脱粒装备的改进与优化具有重要的参考价值.

1 材料与方法

1.1 材料的制备

试验所用材料为存放4个月的蓝亚麻的蒴果，由于成熟期蓝亚麻蒴果室背开列，蒴果头部会炸裂变成锯齿状的圆形缺口，会导致测试结果不准确，因此在试验前首先挑选室背还未开裂且籽粒饱满大小均匀的蓝亚麻蒴果，用镊子将其整齐排列粘附到塑料板上，尽可能的减小蒴果间的空隙，并利用大小不同的蒴果和特别挑选的蒴果壳碎片尽可能的使蒴果粘附板表面平整，制作蒴果横向排列粘附(蒴果腰部与塑料板粘附)板2个，如图1-A、B所示，制作蒴果纵向排列粘附(蒴果头部和尾部分别与塑料板粘附)板4个，如图1-C、D、E、F所示.

图1 制备的试验材料Figure 1 Test materials by prepared

1.2 设备与仪器

试验设备与仪器主要有：自制的接触力学参数测试装置，该装置主要由步进电动机、滑轮、Q235钢斜面板、水平仪、步进电机驱动器、调速器等组成，如图2所示，试验用接触底板为40钢板和Q235钢板；游标卡尺(精度0.02 mm)，桂林量具刃具有限责任公司；电热鼓风干燥箱，天津市三水科学仪器有限公司；电子天平(精度0.01 g),上海越平科学仪器(苏州)制造有限公司.

1:制备的试验材料;2:40钢板;3:步进电机;4:钢丝拉索;5:滑轮;6:Q235钢板;7:水平仪;8:摄像机:9:步进电机驱动器;10:调速器.1:Test materials prepared;2:40 Steel plate;3:Stepping motor;4:Wire and cable ;5:The pulley ;6:Q235 Steel plate;7:Level gauge;8:Video camera;9：Stepper motor driver;10:Speed controller.图2 接触力学参数测试装置Figure 2 Testing device for contact mechanics parameters

1.3 试验方法

1.3.1 千粒质量、含水率的测量方法 蒴果的含水率用直接干燥法测量，称取1 000粒蒴果样品，记录其质量为m1,然后将蒴果样品放入电热鼓风干燥箱中，设定温度为105 ℃，放置72 h，每隔3 h测量一次样品的质量，直到样品的质量不再变化时为止，记录此时样品的质量为m2.则千粒质量为m1，含水率为：(m1-m2)/m1×100%.

1.3.2 三轴尺寸、密度的测量方法 蒴果的三轴尺寸运用游标卡尺测量，随机选取成熟期蓝亚麻蒴果100颗，用游标卡尺(精度0.02 mm)进行三轴尺寸测定,运用MATLAB 2016软件对所测的蒴果三轴尺寸数据进行处理与分析，得到蒴果三轴尺寸的概率密度分布规律，蒴果的三轴尺寸是指蒴果腰部横截面直径(Wd)和室背高度(Td)，如图3所示.由蒴果三轴尺寸便可求得蒴果的体积，利用上一步测得的蒴果质量，进而求得蒴果的密度.

1.3.3 泊松比、剪切模量的确定方法 泊松比和剪切模量通过查阅文献法得到.

1.3.4 实际休止角的测量方法 蓝亚麻蒴果实际休止角的测量用无底圆筒法，如图4-A所示，试验时将Q235钢制的无底圆筒(内壁直径32 mm×高度70 mm)垂直放置于Q235钢板的平面上，圆筒内注入500 粒蒴果，待圆筒内的蒴果稳定后使用二相四线式混合步进电机以30 mm/s的速度匀速向上提升圆筒，此时圆筒内的蒴果颗粒在重力作用下自然下落堆积，形成的圆锥体底角即蒴果休止角.

图3 蓝亚麻蒴果三轴尺寸图Figure 3 Triaxial dimensions of capsule of blue flax

1:Q235钢制无底圆筒;2:步进电机;3:蓝亚麻蒴果堆;4:Q235钢基板.1:Bottomless cylinder of Q235 steel;2:Stepping motor;3：Capsule heap of blue flax;4：Base plate of Q235 steel.图4 蓝亚麻蒴果休止角测定试验Figure 4 Capsule angle of repose determination of blue flax

如 图4-B、C、D所示，待蒴果堆完全静止稳定后，分别测量蓝亚麻蒴果堆的高度H、蒴果堆在X轴和Y轴的宽度W1和W2，则蒴果堆的休止角在+X、+Y方向上的2个角度值分别为ɑ1=arctan(2H/W2)和ɑ2=arctan(2H/W1)，2个方向上休止角的平均值即为每次试验测定蒴果的休止角，重复进行5次试验，试验结果取平均值.

1:塑料刻度尺;2:Q235钢制接触底板;3:摄像机.1:Plastic scale;2:Contact plate of Q235 steel;3：The camera.图5 碰撞恢复系数的测量过程Figure 5 Measurement process of collision recovery coefficient

1.3.6 静摩擦系数的测量方法 本研究采用斜面法测量静摩擦系数[10]，试验时利用水平仪选择一个理想水平面并将接触力学参数测试装置放置在水平面上，先后分别将制作好的蒴果横向排列粘附板和蒴果纵向排列粘附板作为滑板放置于测试装置Q235钢板表面上，测试装置Q235钢板一端固定，利用二相四线式混合步进电动机缓慢匀速拉升测试装置的Q235钢板另一端，待其上的蒴果粘附板开始向下滑动时停止拉升并在垂直于测试装置斜面板倾角方向上进行拍照，将图片导入AUTOCAD2007软件中，标识Q235钢板与水平面之间的夹角，如图6所示.每个粘附板重复试验15次，根据f=tanɑ计算得到蒴果横向表面和纵向表面分别与Q235钢板之间的静摩擦系数，试验结果取蒴果2个表面静摩擦系数的平均值.测量蒴果与40钢板间静摩擦系数时，将40钢板作为接触底板固定在测试装置的Q235钢板上面；测量蒴果间静摩擦系数时，只需将40钢板替换为蒴果粘附板.蒴果与其他材料及蒴果之间的静摩擦系数测定过程分别如图6～7所示.

图6 蓝亚麻蒴果与其他材料间的静摩擦系数的测量过程Figure 6 Measurement process of static friction coefficient between capsule and other materials

1:下滑板;2:上滑板.1:Bottom slide;2:Upper slide.图7 蓝亚麻蒴果间的静摩擦系数的测量过程Figure 7 Measurement process of static friction coefficient between capsules of blue flax

1.3.7 滚动摩擦系数的测量方法 本研究通过斜面法测量滚动摩擦系数，如图8所示，试验时随机选取10颗蓝亚麻蒴果.每次试验时从10颗蓝亚麻蒴果中任选一颗放置在接触力学参数测试装置的Q235钢板表面，缓慢匀速提升Q235钢板的另一端，当蒴果开始滚动时停止提升，每颗蒴果重复5次试验，每次试验结束时用相机垂直于测试装置Q235钢板与水平面形成的角度方向拍照，将照片导入AUTOCAD2007软件中利用直线标注工具进行角度的标注.测量蒴果与40钢板之间的滚动摩擦系数时，将40钢板作为接触底板固定在测试装置的Q235钢板上面；测试蒴果间的滚动摩擦系数时，只需将40钢板替换为蒴果粘附板.

1:蓝亚麻蒴果;2:蒴果粘附板;3:40钢;4:Q235钢.1:Capsule of blue flax;2:Capsule adherent plate;2:40 steel;4:Q235 steel.图8 蓝亚麻蒴果间动摩擦系数的粗测过程Figure 8 Rough measurement of the rolling friction coefficient between capsules of blue flax

1.3.8 仿真逼近预测法 通过上述物理试验方法得到的泊松比、剪切模量、静摩擦系数和动摩擦系数仅仅是一个范围值，为了得到更为精确的数值，还需运用仿真逼近预测法进行进一步的试验.在离散元仿真试验中，输入的参数不同，形成的休止角也各异，因此以物理试验得到的各基本物理参数和接触参数为影响因数，以休止角为试验目标，设计优化仿真试验，通过仿真试验休止角逼近实际休止角的方法来得到各因数的最优参数组合，进而对参数进行标定，确定各基本物理参数和接触参数的理想数值.

1.4 数据分析

运用MATLAB 2016软件进行数据处理与绘图，运用Minitab 18软件进行Plackett-Burman试验设计，运用SAS 9.1软件进行试验数据的处理与分析.

2 结果与分析

2.1 基本物理参数的确定

首先测得本试验所用蓝亚麻蒴果的密度、千粒质量和含水率分别为230 kg/m3、17.9 g和8.60%，实际休止角的平均值为28.12 °.利用MATLAB 2016软件对所测的蒴果三轴尺寸数据进行处理与分析，得到蓝亚麻蒴果三轴尺寸的概率密度分布图，如图9所示，由概率密度分布图可知蓝亚麻蒴果的三轴尺寸基本呈正态分布规律，其中Wd～N(5.56,0.472),Td～N(4.52,0.302).

实际蒴果体积测量困难，可借助蒴果三轴尺寸计算分析出蒴果体积分布的概率密度及标准差，计算依据球体体积计算公式V=(4πWd3)/24.蓝亚麻蒴果三维体积分布规律如图9所示，由图可得蓝亚麻蒴果三维体积基本均呈正态分布，V～N(77.81,20.212).

表1 蓝亚麻蒴果三轴尺寸及体积

图9 蓝亚麻蒴果三轴尺寸概率密度分布图Figure 9 Probability density distribution of triaxial dimensions of capsule of blue flax

2.2 接触力学参数的预测

在用斜面法测量蓝亚麻蒴果与Q235钢、40钢之间的静、滚动摩擦系数时，蒴果的受力分析如图10所示，利用摩擦因数f与斜面倾斜角度ɑ的关系f=tanɑ，便可较精确地测量出蒴果与Q235钢、40钢之间的静、滚动摩擦系数.而用斜面法测量蒴果间的静、滚动摩擦系数时由于蒴果粘附板表面不平整，所以其受力较复杂，受力分析如图11所示，所以摩擦系数f与斜面倾斜角度ɑ的关系只能近似用f=tanɑ来表示，因此用斜面法测量蒴果间的静、滚动摩擦系数时得到是一个范围值，但为后续进行的参数标定仿真试验提供了参数取值范围.在用自由跌落法测蒴果间的碰撞恢复系数时，蒴果下落后接触蒴果粘附板，并反作用于蒴果向上的弹力，由于蒴果粘附排上各个蒴果对下落的蒴果的综合弹力方向近似竖直向上，如图12所示，所以测得的蒴果间及蒴果与Q235钢、40钢间的碰撞恢复系数都较为精确.

1:蒴果粘附板;2:单个蒴果;3:Q235或40钢板1:Capsule adherent plate;2:A single capsule;3:Q235 or 40 steel plate.图10 斜面法测量蓝亚麻蒴果与Q235或40钢板的静、滚动摩擦系数时蒴果的受力分析示意图Figure 10 Schematic diagram of stress analysis of capsule for measuring static and rolling friction coefficients of blue flax capsule and Q235 or 40 steel plate by inclined plane method

1:上滑板(蒴果粘附板);2单个蒴果;3:下滑板(蒴果粘附板);4:Q235钢板.1:Upper slide;2:A single capsule;3:Bottom slide;4:Q235 steel plate.图11 斜面法测量蓝亚麻蒴果间的静、滚动摩擦系数时蒴果的受力分析示意图Figure 11 Schematic diagram of stress analysis of capsule when measuring static and rolling friction coefficients between capsules of blue flax by inclined plane method

图12 用自由跌落法测量蓝亚麻蒴果间及与Q235或40钢板的碰撞恢复系数时蒴果的受力分析示意图Figure 12 Schematic diagram of stress analysis of the capsule when the collision recovery coefficient between blue flax capsules and Q235 or 40 steel plate was measured by free drop method

通过上述分析，可以确定蒴果与Q235钢材之间、蒴果与40钢材、蒴果之间的静摩擦系数分别为0.68、0.62、0.39～0.88；蒴果与Q235钢板之间、蒴果与40钢板之间、蒴果之间的滚动摩擦系数分别为0.23、0.28和0.16～0.36；蒴果与Q235钢板、蒴果与40钢板及蒴果之间的碰撞恢复系数分别为0.32、0.31、0.28.

2.3 泊松比和剪切模量的预测

由于蓝亚麻蒴果颗粒较小，基于现有的试验设备及条件难以直接测量其泊松比和剪切模量，本研究通过查阅文献法来初步预测其数值，为后续的参数标定仿真试验提供参数取值范围.蓝亚麻蒴果内部的籽粒围绕蒴果中心呈纵向圆周排列，且成熟期饱满的蒴果内部籽粒与蒴果壳紧密接触，因此蓝亚麻蒴果的泊松比和剪切模量与其籽粒有很大程度的相似性，又蓝亚麻蒴果与胡麻蒴果具有很好的相似性，根据文献[4]测定的胡麻籽粒的泊松比和剪切模量可初步确定蓝亚麻蒴果的泊松比的取值范围为0.348～0.458,蓝亚麻蒴果的剪切模量取值范围为131.27～197.97 MPa.

2.4 离散元模型的建立与参数标定

2.4.1 蓝亚麻蒴果模型的建立 蓝亚麻蒴果近似圆球状，根据蒴果的三轴统计尺寸和形状，利用三维绘图建模软件SolidWorks对蓝亚麻蒴果的三维结构进行建模，并将蒴果三维结构模型导入离散元仿真软件EDEM中去,用多个不等直径的球体对其进行填充，直到蒴果三维模型被紧密填充、无可填充空间为止.蓝亚麻蒴果离散元模型采用1个半径为2.32 mm、10个半径为1.8 mm和10个半径为2 mm 3种不同粒径的球体颗粒填充，得到蒴果的离散元模型腰部横截面直径为5.08 mm、室背高度为4.60 mm.蓝亚麻蒴果的实物、三维结构模型和离散元模型分别如图13-A、B、C所示.

A：蒴果实物；B：蒴果三维模型；C：蒴果离散元模型.A:Capsule physical；B：Three-dimensional model of capsule；C：Capsule discrete element model.图13 蓝亚麻蒴果的实物、三维模型和离散元模型Figure 13 Physical,3D and discrete element model of capsules

2.4.2 EDEM软件仿真工作参数设置 由文献[4]可知，利用EDEM 2018软件进行蓝亚麻蒴果休止角仿真试验，仿真时间步长的选取至关重要，过大会导致颗粒发生爆炸式发散，过小会使计算量成倍增加.时间步长计算如以下公式所示.

式中：Δt为时间步长，s；R为蓝亚麻蒴果的球体直径，mm；vs为蓝亚麻蒴果下落过程中的最大速度，mm/s；ρ为蓝亚麻蒴果的密度，kg/mm3；G为蓝亚麻蒴果的剪切模量，MPa.

休止角仿真模型的建立依据休止角测定装置实际参数设置，如图14所示.兼顾仿真试验的可行性和软件运行的高效率，仿真中均采用固定尺寸的蒴果离散元模型.设置仿真总时间5 s，时间步长为1×10-5s，网格尺寸为2.5 R min.

设置直径36 mm、高度70 mm的Q235钢制无底圆筒用于形成蒴果休止角,以Q235钢板为蒴果堆积基板，在圆筒上口位置建立蒴果颗粒动态生成工厂，蒴果模型与蒴果模型之间、蒴果模型与圆筒之间都采用EDEM2018软件内置的Hertz-Mindlin(no slip)接触模型.每秒产生500颗蒴果模型，并在重力的作用下自由下落，经过1.1 s的时间，当容器内填充完500粒蒴果达到稳定状态后,设置提升速度为30 mm/s的速度匀速提升圆筒，待蓝亚麻蒴果种群失去圆筒的束缚作用后自然散开形成堆积，在基板上静止后堆积形成的圆锥体底角便是蓝亚麻蒴果堆积休止角.

图14 蓝亚麻蒴果休止角仿真试验模型Figure 14 Simulation model of the angle of repose for capsule

2.4.3 Plackett-Burman试验 蓝亚麻蒴果休止角仿真试验以物理试验法测定和查阅文献法得到的蓝亚麻蒴果实际参数值为依据.为减小试验次数，缩小试验时间，应用Minitab18软件进行Plackett-Burman试验设计，筛选出对蓝亚麻蒴果仿真试验休止角影响显著的参数.对表2中4个不确定的参数进行Plackett-Burman试验设计，将4个参数的最小、最大值分别编码为-1、+1 水平.

表2 Plackett-Burman试验参数表

试验设置了1个中心点(取高、低水平的中间值作为0水平)，共需进行13次蓝亚麻蒴果休止角仿真试验.利用EDEM 2018软件自带的量角工具在蓝亚麻蒴果离散元模型堆的+X、+Y2个方向测量休止角，如图15所示，取其平均数，试验设计及结果如表3所示.

利用SAS 9.1数据分析软件对试验结果进行方差分析，可得到各参数的显著性结果，如表4所示.方差分析结果显示蒴果间的滚动摩擦系数(X4)对仿真休止角影响显著.

2.4.4 最陡爬坡试验 将Plackett-Burman试验得到的显著性参数，进行最陡爬坡试验，结果如表5所示.首先进行蒴果堆积休止角仿真试验的预试验Ⅰ，通过预试验可快速逼近显著性参数最优区间，显著性参数根据设定好的步长在其取值范围内逐渐增加，随着显著性参数取值的逐步增大，仿真试验休止角和实际休止角的相对误差呈现增大的趋势，蒴果实际休止角(28.12 °)处在1号和2号试验对应的仿真试验休止角之间；1号与2号试验之间的显著性参数继续按照一定的步长取值，进行二次爬坡休止角仿真试验Ⅱ，随着显著性参数的取值逐步增大，仿真试验休止角和实际休止角的相对误差呈现先增大后减小的趋势，显著性参数取0.18时仿真试验休止角与实际休止角的相对误差较小，由此可知蒴果间滚动摩擦系数较优取值为0.18.

图15 蓝亚麻蒴果离散元模型堆Figure 15 Discrete element model pile of capsule of blue flax

表3 Plackett-Burman试验设计及结果

表4 Plackett-Burman试验结果方差分析

表5 最陡爬坡试验及结果

在上述休止角仿真试验中蓝亚麻蒴果的泊松比取Plackett-Burman试验的中间水平，即0.40；蓝亚麻蒴果的剪切模量只影响仿真中蒴果模型间的碰撞受力，对其运动状态影响较小，为了提高仿真运行效率，蓝亚麻蒴果的剪切模量应尽可能小于实际数值，本研究中取值为164 MPa；其他非显著性接触参数取物理试验所测得的平均值：蒴果间的碰撞恢复系数为0.28,蒴果与Q235钢碰撞恢复系数为0.32，蒴果间的静摩擦系数为0.64，蒴果与Q235钢的静摩擦系数为0.68，蒴果与Q235钢滚动摩擦系数0.24.

2.4.5 仿真试验最优参数组合的验证 通过上述进行的Plackett-Burman试验和最陡爬坡试验，获得了当仿真试验休止角与实际休止角的相对误差较小状态下的各仿真参数值.利用仿真试验获得的最优参数值进行3次重复的休止角仿真试验，测得仿真试验的休止角分别为28.39 °、28.24 °、28.57 °，3次仿真试验结果的平均值为28.40 °，与实际休止角28.12 °的相对误差为0.95%.蓝亚麻蒴果实际休止角与仿真试验休止角的对比如图16所示.

图15 蓝亚麻蒴果实际休止角与仿真休止角的对比Figure 15 Comparison between actual and simulated angles of repose blue flax capsule

2.5 离散元模型仿真参数的确定

本试验通过物理试验测定法和基于离散元仿真软件EDEM 2018的仿真逼近预测法确定的蓝亚麻蒴果离散元模型各个仿真参数如表6所示.

表6 蓝亚麻蒴果离散元仿真各参数确定的数值

3 结论

1) 通过物理试验法测得了蓝亚麻蒴果基本参数.蓝亚麻蒴果的三轴尺寸呈正态分布，密度为230 kg/m3，蒴果之间、蒴果与40钢、Q235钢之间的碰撞恢复系数分别为0.28、0.31、0.32，静摩擦系数分别为0.39～0.88、0.62、0.68.滚动摩擦系数分别为0.16～0.36、0.28、0.23.通过查阅文献法得到蓝亚麻蒴果的泊松比为0.35～0.46，剪切模量为131～197 MPa.

2) 以物理试验法测得的蓝亚麻蒴果基本参数为依据，基于离散元仿真软件EDEM 2018对基本参数中4个不确定参数进行Plackett-Burman 试验，方差分析结果表明，蒴果间滚动摩擦系数对仿真休止角影响显著.通过仿真逼近预测法对蒴果间滚动摩擦系数进行最陡爬坡试验，在仿真试验休止角与实际休止角相对误差较小的情况下得到蓝亚麻蒴果间的滚动摩擦系数为0.18.

3) 为验证所确定的蓝亚麻蒴果基本参数的可靠性，进行了休止角仿真试验验证，对蒴果仿真试验休止角与实际休止角进行了对比，结果显示两者的相对误差为0.94%；本研究标定的蓝亚麻蒴果离散元模型仿真参数对于利用EDEM软件优化亚麻(胡麻)脱粒装备的工作性能参数具有重要的参考价值.