苜蓿秸秆压缩仿真离散元模型参数标定

马彦华，宋春东，宣传忠，王昊毅，杨 帅，武 佩

（内蒙古农业大学机电工程学院，内蒙古自治区草业与养殖业智能装备工程技术研究中心，呼和浩特 010018）

0 引 言

苜蓿具有丰富的营养价值，是牲畜喜食的一种饲草料，尤其是将苜蓿压缩成颗粒饲料，不但可以增加密度，方便贮存和运输，而且还可以极大地提高牲畜的适口性，增加采食量[1-3]。因此苜蓿颗粒饲料具有广阔的市场前景。

苜蓿颗粒饲料是在模具中依靠压缩装置对揉碎后的苜蓿秸秆施压，克服物料变形阻力、物料间摩擦以及物料与模具内壁摩擦而成型。因此，研究模具中压力对物料的作用过程，有助于压缩工艺参数的合理选取。离散元仿真分析法在土壤、粉末等颗粒物料压缩试验中的应用[4-6]，为苜蓿秸秆致密成型过程应力分析提供了有效方法。离散元仿真分析需要较多参数，主要包括苜蓿秸秆物理、力学特性参数，苜蓿秸秆间、苜蓿秸秆与模具间的接触参数等[7-9]。客观、准确地获得仿真参数是离散元法应用于苜蓿秸秆压缩成型研究的前提。

有关苜蓿茎秆力学特性的研究，国内外学者主要围绕着拉伸强度、剪切强度、抗弯强度性能指标进行了大量试验。张涵[10]对紫花苜蓿茎秆进行了剪切试验，探讨了苜蓿茎秆剪切力、剪切能量的变化趋势，以及含水率与加载速度耦合作用对剪切力等参数的影响。赵春花等[11]测试了苜蓿茎秆的拉伸和剪切应力—应变规律，发现苜蓿茎秆的强度和刚度主要取决于机械组织的厚度、维管束的数量，并给出了低纤维苜蓿平均弹性模量和最大剪切力。Nazari 等[12]的研究表明，苜蓿茎秆的抗弯强度随含水率降低而逐渐增大，并给出苜蓿茎秆不同部位的弯曲应力。这些研究结果为苜蓿秸秆离散元模型的建立提供了基本力学参数。

目前关于苜蓿秸秆接触参数尚未有可借鉴的研究成果，但是接触参数可以依靠颗粒物料休止角试验测得。休止角表征物料流动、摩擦等特性，能够从宏观角度反映物料群综合作用结果[13-14]。国内外众多学者基于休止角试验对物料的接触参数进行了离散元参数标定。彭飞等[15]利用注入截面法测得颗粒饲料的休止角，并通过响应面法得到了大猪料最优的仿真参数组合。吴佳胜等[16]利用提升空心圆筒堆积试验法，并结合图像处理技术测量出前胡种子的休止角，通过Plackett-Burman 筛选试验、响应面方差分析法以及回归方程优化法标定了前胡种子的各物性参数。Lee 等[17]设计了4 个仿真摩擦试验，分别测定了物料与物料间、物料与容器间的静、滚动摩擦系数，并以休止角为响应值验证了所测系数的准确性。刘凡一等[18]采用圆筒提升法和MATLAB 图像处理技术获得小麦休止角，通过Plackett-Burman 试验、最陡爬坡以及Box-Behnken 试验对参数进行了离散元标定。

总结以上的研究发现，国内外众多学者大都参考通用颗粒材料模型数据库（Generic EDEM Material Model Database，GEMM）给定的参数区间，或对个别参数进行试验，来设置物料离散元参数的各优化仿真试验初始范围值，鲜有采用比较全面的物理试验对物料参数进行测定，并以实测参数为依据进行物料的离散元仿真参数标定。本文采用6 个接触参数测定试验，以此确定优化试验的参数区间，依次通过Plackett-Burman 试验、最陡爬坡试验以及Box-Behnken 试验对苜蓿秸秆的离散元仿真参数进行标定，并将仿真休止角与物理试验休止角利用T检验进行对比验证，以确定苜蓿秸秆离散元仿真参数值。通过以上试验得到的苜蓿秸秆离散元仿真参数，可为采用离散元法从微观上深入研究物料块的成型机理以及物料层中应力传递规律提供参数依据。

1 材料与方法

1.1 材料与休止角物理试验

试验原料选用2018 年9 月产自于呼和浩特市内蒙古农业大学自治区生物制造重点实验室试验田（106°42′E，37°35′N）的紫花苜蓿，清除杂质、晾晒，保证其含水率为18%。休止角物理试验时，将苜蓿原料用9RS-60 型揉碎机揉碎，采用10 目和12 目的标准筛对揉碎后的苜蓿秸秆进行筛分[19]，得到粒径为1.7～2.0 mm 的苜蓿秸秆。

根据文献[20]的方法进行苜蓿秸秆休止角试验，利用JA5003B 型电子天平称取8 g 筛选好的苜蓿秸秆，在漏斗上口中心位置使其自由下落，如图1。待所有苜蓿秸秆落在基板上静止后，用相机拍摄苜蓿秸秆堆，得到苜蓿秸秆堆的正视照片，应用Origin 2018 软件中的图像数字化工具进行处理。试验重复10 次，得到休止角的平均值为38.88°。

图1 FT-104B 休止角测定仪 Fig.1 Repose angle tester FT-104B

1.2 接触参数的测定

需要测定的接触参数包括：苜蓿秸秆和苜蓿秸秆之间、苜蓿秸秆和成型模具之间的碰撞恢复系数、静摩擦系数、滚动摩擦系数。测定的参数值为后续离散元仿真输入参数提供数据支持。

接触参数物理试验采用的苜蓿秸秆取自整株苜蓿秸秆的中间无茎节部分，其密度为256 kg/m3。压缩成型模具材料为45 钢，其泊松比为0.3，密度为7 800 kg/m，剪切模量为70 000 MPa[21]。接触参数试验时选材料为45 的圆钢，按照成型模具的工艺要求对其测试端面进行处理。

1.2.1 碰撞恢复系数的测定

碰撞恢复系数表征了被测物体在发生碰撞后，恢复到原来状态的能力，并定义为被测物体碰撞分离后的法向速度与碰撞前法向速度之比[22-23]。本文采用自由落体试验测定碰撞恢复系数，试验装置如图2 所示。

将高速摄像机放置在试验支架正前方，调节三脚架使其保持水平，高速摄像机与试验支架的距离为780 mm，如图2a 所示。为便于高速摄像分析软件TEMA 对苜蓿秸秆运动过程进行捕捉，需在试验支架上设置3 个参考点（图2b 中的点1、2、3）。由于单个苜蓿秸秆重量较小，在其自由下落过程中，为兼顾空气阻力的影响及碰撞反弹效果，要求苜蓿秸秆有一定的下落高度，经过多次反复试验确定1、2 两点距离为245 mm；1、3 两点的距离要大于苜蓿秸秆的运动范围，取其值为165 mm。

图2 碰撞恢复系数测定装置 Fig. 2 Measuring set-up for collision recovery coefficient

测定苜蓿秸秆间的碰撞恢复系数时，为了使苜蓿秸秆下落时能够与苜蓿秸秆碰撞，要保证试验区域具有一定的、均匀的接触面积，故由长度60 mm，直径2.06 mm的无茎节苜蓿秸秆组成苜蓿秸秆排。将单个苜蓿秸秆自苜蓿秸秆排正上方245 mm 处自由下落，当苜蓿秸秆与苜蓿秸秆排碰撞后快速弹起，利用Pco.dimax S 型高速摄像机对苜蓿秸秆运动的整个过程进行拍摄，同时将拍摄的苜蓿秸秆运动画面实时存入电脑。测定苜蓿秸秆与45 钢的碰撞恢复系数时，将苜蓿秸秆排去掉，其它过程同上。

应用TEMA 软件对拍摄的苜蓿秸秆运动过程进行分析，将分析结果导入Origin 2018 软件，得到如图3所示的苜蓿秸秆与苜蓿秸秆、苜蓿秸秆与45 钢碰撞运动位移-时间曲线。

假设苜蓿秸秆在下落过程中只受重力，由动能定理得到苜蓿秸秆在碰撞前的法向速度

碰撞分离后的法向速度

由碰撞恢复系数的定义，得到碰撞恢复系数为

式中1v 是碰撞前的法向速度，m/s；2v 是碰撞后的法向速度，m/s；1h 为碰撞前的下落高度，mm；2h 为碰撞后弹起的最大高度，mm。

图3 苜蓿秸秆与被测物体碰撞的运动位移-时间曲线 Fig.3 Displacement-time curves of collision between alfalfa straw and measured object

经多次试验，得到苜蓿秸秆间的碰撞恢复系数范围为0.1～0.3，苜蓿秸秆与45 钢的碰撞恢复系数范围为0.1～0.6。

1.2.2 静摩擦系数的测定

利用CNY-1 型斜面仪对苜蓿秸秆和苜蓿秸秆、苜蓿秸秆和45 钢间的静摩擦系数进行测定[24]。其中苜蓿秸秆选择长度为30 mm，直径为2.06 mm 的无茎节部分。试验前，将斜面仪底座水平放置，调节斜面仪测试平面使其指针指向0°刻度线。

苜蓿秸秆与45 钢的静摩擦系数测定试验如图4a 所示。将45 钢放置在斜面仪测试平面上，为了防止出现相对滑动，使用压敏胶带将45 钢下表面与斜面仪测试平面粘结在一起，然后将苜蓿秸秆沿斜面仪测试平面长度方向放置在45 钢的测定平面上。逆时针缓慢转动斜面仪测试平面，直至观察到苜蓿秸秆在45 钢测定平面上出现滑动，停止转动，记录此时斜面仪指针指示角度。通过式（4）计算得到苜蓿秸秆与45 钢平面的静摩擦系数。测定苜蓿秸秆之间的静摩擦系数时，只需要将45 钢替换为苜蓿秸秆即可。经多次试验，得到苜蓿秸秆间的静摩擦系数范围为0.3～0.6，苜蓿秸秆与45 钢的静摩擦系数范围为0.2～0.8。

式中1μ 为静摩擦系数，1φ 为苜蓿秸秆在测定平面上即将滑动时斜面仪指针指示角度(°)。

1.2.3 滚动摩擦系数的测定

本文从能量守恒原理出发，来进行滚动摩擦系数的测定[25-26]。

调节CNY-1 型斜面仪测试平面使其指针指向0°刻度线，将苜蓿秸秆沿斜面仪测试平面宽度方向放置在45钢的测定平面上，逆时针缓慢转动斜面仪测试平面，记录秸秆开始滚动时斜面仪指针指示角度，记为2φ ，如图 4b 所示。此时苜蓿秸秆在重力作用下，沿测试斜面滚落，滚落过程由能量守恒定律可得

式中W 为苜蓿秸秆在滚动中损失的能量，即滚动摩擦消耗的能量，J；pE 为苜蓿秸秆初始的重力势能，J；kE 为苜蓿秸秆在终了时刻的动能，J。

由于单个苜蓿秸秆质量非常小，其滚落到测试平面底端的动能很难获得，所以本试验采用初始滚动时刻内近似能量守恒获得的滚动摩擦系数来预估真实值。苜蓿秸秆在斜面仪上的受力分析如图4b 所示，当一重力为G的苜蓿秸秆滚动微小角度时，在45 钢测定平面上滚过的距离为δ ，苜蓿秸秆对测定平面的正压力滚动摩擦力滚动摩擦力做的功重力势能

由于在开始滚动的一微小角度内，速度由零增加，速度变化极小，此时假设苜蓿秸秆的动能为由以上分析，可得到滚动摩擦系数为

式中2μ 为滚动摩擦系数，2φ 为苜蓿秸秆在测定平面上开始滚动时斜面仪指针指示角度(°)。

图4 苜蓿秸秆与45 钢摩擦系数的测定试验 Fig.4 Tests of friction coefficient between alfalfa straw and 45 steel

测定苜蓿秸秆之间的滚动摩擦系数时，只需要将45钢替换为苜蓿秸秆排即可。经多次试验，得到苜蓿秸秆间的静摩擦系数范围为0～0.3，苜蓿秸秆与45 钢的静摩擦系数范围为0～0.3。

采用上述测试方法虽不能准确给出滚动摩擦系数，但能预估苜蓿秸秆-苜蓿秸秆、苜蓿秸秆-45 钢间滚动摩擦系数的最大值，为后续仿真试验中滚动摩擦系数区间的选择提供依据。

2 苜蓿秸秆离散元模型的建立及关键参数标定

在离散元仿真试验中，输入的参数不同，形成的休止角也各异，因此以物理试验得到的休止角为试验目标，设计优化仿真试验，对参数进行标定。

2.1 苜蓿秸秆离散元模型的建立

苜蓿秸秆近似长条状，经两层筛筛分后随机选取 100 根苜蓿秸秆，利用游标卡尺测量它们的长度和直径，得到其平均直径为2 mm，平均长度为10 mm。因此，参考EDEM 2018 软件自带的长条状模型，利用半径为1 mm的9 个基本球体（基本球体半径为Rmin）相贯组合成苜蓿秸秆离散元模型，如图5 所示。

图5 苜蓿秸秆离散元模型 Fig.5 Discrete element model of alfalfa straw

2.2 EDEM 软件工作参数设置

休止角仿真模型的建立依据休止角测定仪实际参数设置，如图6 所示。兼顾试验的可行性和软件运行的高效率，仿真中均采用固定尺寸的苜蓿秸秆离散元模型。设置仿真总时间7 s，时间步长为2.346 1×10-5s，网格尺寸为2.5Rmin。

图6 苜蓿秸秆休止角仿真模型 Fig.6 Simulation model of the repose angle for alfalfa straw

在漏斗上口位置建立动态生成工厂，每秒产生5 g 苜蓿秸秆模型，并在重力的作用下自由下落，经过7 s 的时间在基板上静止，形成苜蓿秸秆堆。苜蓿秸秆模型与苜蓿秸秆模型之间、苜蓿秸秆模型与漏斗之间分别采用EDEM 2018 软件内置的EEPA（Edinburgh Elasto-Plastic Adhesion Model）接触模型、Hertz-Mindlin（no slip）接触模型。

2.3 仿真模型的标定试验

2.3.1 Plackett-Burman 试验

本研究以物理试验测定值为依据，并查阅相关文献[10，27-28]，可知离散元仿真必要参数的取值范围。为减小试验次数，尽快得到确定的参数组合，应用Minitab 18 软件进行Plackett-Burman 试验设计，筛选出对试验结果影响显著的参数（表1）。对表1 中8 个不确定的参数进行Plackett-Burman 试验设计，将8 个参数的最小、最大值分别编码为-1、+1 水平。

表1 Plackett-Burman 试验参数表 Table 1 Test parameters of Plackett-Burman

该研究设置了1 个中心点（取高、低水平的中间值作为0 水平），共进行13 次试验。利用EDEM 2018 软件自带的量角工具在苜蓿秸秆堆的+X、+Y 两个方向测量休止角，如图7 所示，取其平均数，试验设计及结果如表2 所示。

图7 苜蓿秸秆离散元模型堆 Fig.7 Discrete element model pile for the alfalfa straw

表2 Plackett-Burman 试验设计及结果 Table 2 Design and results of Plackett-Burman test

利用Minitab 18 软件对试验结果进行方差分析，可得到各参数的显著性结果，如表3 所示。苜蓿秸秆-苜蓿秸秆静摩擦系数（X4）、苜蓿秸秆-苜蓿秸秆滚动摩擦系数（X5）、苜蓿秸秆-45 钢静摩擦系数（X7）对仿真休止角影响显著。

表3 Plackett-Burman 试验结果方差分析 Table 3 Variance analysis of Plackett-Burman test results

2.3.2 最陡爬坡试验

将Plackett-Burman 试验得到的3 个显著性参数，进行最陡爬坡试验，可快速逼近最优参数区域。各显著性参数根据设定好的步长在其取值范围内逐渐增加，进行休止角仿真试验，结果如表4 所示。随着显著性参数取值的逐步增大，仿真休止角和试验休止角的相对误差先减少后增大；在3 号试验对应的参数组合下，得到的相对误差最小，由此可知最优的参数水平范围在3 号试验附近。

表4 最陡爬坡试验设计及结果 Table 4 Design and results of steepest climbing test

仿真试验中苜蓿秸秆的泊松比取Plackett-Burman 试验的中间水平，即0.4；苜蓿秸秆剪切模量只影响仿真中苜蓿秸秆模型间的碰撞受力，对其运动状态影响较小[20]，为了提高仿真运行效率，苜蓿秸秆的剪切模量应尽可能小于实际数值，本研究中取值为5 MPa；其他非显著性接触参数取物理试验所测得的平均值：苜蓿与苜蓿碰撞恢复系数为0.11，苜蓿与45 钢碰撞恢复系数为0.16，苜蓿与45 钢滚动摩擦系数0.24。

2.3.3 显著接触参数的Box-Behnken 试验及分析

取最陡爬坡试验中的3 号试验参数为中间水平（0），4、2 号试验参数分别为高（+1）、低水平（-1），显著接触参数与水平如表5 所示，应用Design Expert 8 软件对其进行Box-Behnken 试验。其它非显著性参数的取值同最陡爬坡试验。显著接触参数的Box-Behnken 试验设计及结果如表6 所示。

表5 显著接触参数水平编码 Table 5 Levels of significant contact parameters

表6 显著接触参数的Box-Behnken 试验设计及结果 Table 6 Box-Behnken experimental design and results of significant contact parameters

利用Design-Expert 8 软件对表6 中的试验结果进行多元回归分析，得到苜蓿秸秆仿真堆积角与3 个显著性参数的二阶回归模型，方程为

该拟合模型的P< 0.0001，决定系数R2=0.981、校正决定系数R2adj=0.955，均接近1，且信噪比（Adeq precision）为22.393，说明该休止角回归模型极显著，可以对目标休止角进行预测；失拟项 P=0.1583>0.05，变异系数C.V.=1.07 % 较低，表明方程拟合程度很好。对该模型试验结果进行回归方差分析，结果如表7 所示，由该模型方差分析结果可知，苜蓿秸秆-苜蓿秸秆静摩擦系数（X4）、苜蓿秸秆-苜蓿秸秆滚动摩擦系数（X5）、苜蓿秸秆-45 钢的静摩擦系数（X7）及苜蓿秸秆-45 钢的静摩擦系数二次方项（X72）均对休止角影响极其显著；苜蓿秸秆-苜蓿秸秆滚动摩擦系数的二次方项（X52）对休止角影响显著；苜蓿秸秆-苜蓿秸秆静摩擦系数二次方项（X42）及其余交互项对休止角的影响不显著。

表7 Box-Behnken 试验设计回归模型方差分析 Table 7 Variation analysis of Box-Behnken design quadratic model

2.3.4 仿真最优参数组合的验证及确定

利用Design Expert 8 软件中的优化模块，以休止角试验值38.88°为目标值，对回归模型式（7）进行寻优，得到一组与物理试验数据相近的参数：苜蓿秸秆-苜蓿秸秆静摩擦系数为0.45，苜蓿秸秆-苜蓿秸秆滚动摩擦系数为0.08，苜蓿秸秆-45 钢的静摩擦系数为0.54，其余非显著性参数的取值同最陡爬坡试验。为检验最优参数组合是否满足仿真试验结果，将上述各参数值输入EDEM 2018 软件中进行仿真，得到苜蓿秸秆堆离散元模型，重复3 次仿真试验，所得苜蓿秸秆休止角分别为39.62°、37.83°、38.69°，平均值为38.71°。应用Origin 2018软件对该样本进行T 检验，得到P= 0.778>0.05，表明仿真休止角与物理试验值无显著性差异。最优参数组合下休止角仿真结果（38.71°）与实际物理休止角结果（38.88°）的相对误差为0.437%，试验对比如图8 所示。

图8 苜蓿秸秆堆休止角仿真试验与物理试验对比 Fig.8 Repose angle comparison between simulation and physical tests of alfalfa straw pile

2.3.5 参数确定

通过Plackett-Burman 试验，得到影响仿真休止角的显著性参数；接着利用最陡爬坡试验优化了显著性参数的取值范围；然后进行Box-Behnken 试验，以目标值对拟合的回归方程寻优，得到最优参数组合；再利用T 检验验证了仿真休止角与物理试验结果无显著性差异；最终得到各参数的确定数值，如表8 所示。

表8 各参数确定数值 Table 8 Determined values of each parameter

3 结 论

1）应用斜面仪和高速摄像机，通过物理试验得到了苜蓿秸秆间、苜蓿秸秆与45 钢间碰撞恢复系数的范围值分别为0.1～0.3、0.1～0.6；苜蓿秸秆间、苜蓿秸秆与45钢间静摩擦系数的范围值分别为0.3～0.6、0.2～0.8；苜蓿秸秆间、苜蓿秸秆与45 钢间滚动摩擦系数的范围值均为0～0.3。

2）以接触参数物理试验结果为依据，进行 Plackett-Burman 试验，方差分析结果表明，苜蓿秸秆-苜蓿秸秆静摩擦系数（X4）、苜蓿秸秆-苜蓿秸秆滚动摩擦系数（X5）、苜蓿秸秆-45 钢静摩擦系数（X7）对仿真休止角影响显著。

3）对显著性参数进行Box-Behnken 试验，建立了休止角和显著性参数的二阶回归模型，以休止角试验值38.88°为目标进行寻优，得到最优解：苜蓿秸秆-苜蓿秸秆静摩擦系数为0.45，苜蓿秸秆-苜蓿秸秆滚动摩擦系数为0.08，苜蓿秸秆-45 钢的静摩擦系数为0.54，优化参数组合下的休止角平均值为38.71°。

4）利用T 检验得到P= 0.778>0.05，表明仿真休止角与物理休止角无显著性差异，验证了最优参数组合的可靠性。