Calibration and test of discrete element simulation parameters of fresh edamame pods

本研究所用鲜食毛豆采自安徽省滁州市某超市，购买后置于4 ℃冷藏保存，3 d内完成试验。试验前，样品恢复至实验室室温环境，同时剔除存在机械损伤、病虫害及明显畸形的豆荚后，采用随机抽样方法选取试样。经统计分析，单粒、双粒和三粒豆荚的占比分别为16%、65%、19%。按此分类各随机抽取100个豆荚样本，测量其长(L)、宽(B)、厚(T)，其中：L为果柄端到底端的最长距离，B为垂直L方向正投影中最大横向距离，T为垂直B且垂直L方向的最大距离，以三粒豆荚为例，测量位置如图1所示。测定豆荚的质量和体积并计算密度；采用干燥法获取豆荚的含水率^[17]。测量结果取平均值(表1)。

Figure 1.  Three-dimensional dimensions of fresh edamame pods

本研究采用压缩试验方法测定豆荚的泊松比^[18]：采用万能试验机，选用直径为100 mm的圆柱形压头。豆荚平放受载，设定加载前后压头移速为60 mm·min⁻¹，加载时移速为5 mm·min⁻¹，最大变形量控制为50%。测量压缩前后豆荚B、T方向的变化量，利用式(1)计算泊松比并取平均值。

式(1)中：μ为泊松比；ε_B、ε_T分别为豆荚的B方向、T方向的应变；ΔB、ΔT分别为豆荚B方向、T方向的绝对变形量(mm)；T₀、B₀分别为豆荚破裂前B、T方向的尺寸(mm)。经计算得到泊松比为0.24。

依托万能试验机，选用直径为4 mm的圆柱形压头。试验时，调整压头至刚好接触豆荚表面，以5 mm·min⁻¹的速率稳定加载至50%最大变形量。通过式(2)和式(3)计算弹性模量和剪切模量，结果取平均值。

式(2)~(3)中：E为弹性模量(MPa)；σ为豆荚的压应力(MPa)；G为剪切模量(MPa)。经计算，豆荚的弹性模量和剪切模量分别为19.94和8.04 MPa。

豆荚离散元模型的构建采用颗粒填充法^[19]。首先于4 mm方格背景下正交采集豆荚L、B、T 3个维度的图像，利用Blender描绘轮廓，生成STL模型(图2A)，并将其导入EDEM软件。以半径为1 mm、平滑度为3的球形颗粒进行填充，建立单粒、双粒和三粒豆荚的离散元模型(图2B)。本研究主要研究豆荚群体堆积与输送过程中的宏观运动行为，试验过程中未观察到豆荚之间存在明显的团聚及黏连现象，具有典型的散体流动特征，因此采用Hertz-Mindlin (no slip) 模型描述颗粒间接触作用。参考农业生物材料离散元参数标定方法^[7−13]，将高含水率、柔软组织及表面茸毛等因素的综合接触效应等效纳入恢复系数、静摩擦系数和滚动摩擦系数等接触参数体系，并结合堆积角试验和输送试验验证模型可靠性。

Figure 2.  Modeling of discrete element model of fresh edamame

为确保离散元仿真结果的准确性，需要合理设定物料的本征参数与接触参数。本征参数直接采用试验实测值；接触参数则通过物理试验与仿真相结合的方法展开标定。针对豆荚采运装备的主体材质，本研究主要对豆荚与304不锈钢、食品级聚氯乙烯(PVC)及豆荚相互间的接触参数进行标定。在标定豆荚间的接触参数时，采用自制的“豆荚板”作为接触面(图3)。制作豆荚板时，选取颗粒饱满、无机械损伤、厚度及弯曲程度相近的豆荚，沿长度方向紧密排列于不锈钢钢板表面，以热熔胶固定豆荚厚度方向的一侧于基板上，保证最大宽度方向与基板平行且相邻豆荚无明显间隙。胶体固化后，剔除局部高度差较大的试样，以确保接触表面平整。

Figure 3.  Fresh edamame contact plate

本研究选取的304不锈钢和PVC材料的参数可由文献和厂家手册获取。304不锈钢的泊松比、剪切模量和密度分别为0.30、28.54 GPa和7 930 kg·m⁻³，PVC的泊松比、剪切模量和密度分别为0.35、1.14 GPa和1 350 kg·m⁻³。

本研究采用自由落体测定碰撞恢复系数^[20]，试验原理如图4A。测定过程中，将豆荚放置于距接触面上方H处自由释放，并记录最大回弹高度h，通过式(4)计算碰撞恢复系数。

Figure 4.  Coefficient of restitution calibration test

式(4)中：H为豆荚自由释放高度(mm)；h为豆荚最大回弹高度(mm)；e为碰撞恢复系数；v₀、v₁分别为碰撞前、后碰撞接触点的瞬时速度(mm·s⁻¹)，g为重力加速度(9.8 m·s⁻²)。

根据上述原理，将豆荚从距接触面高H=230 mm处自由释放，采用相机(镜头距测试面40 cm)记录下落全过程。随后，将视频导入Tracker-6.2.0软件进行运动轨迹追踪，提取豆荚碰撞前后的位移曲线及最大回弹高度(h)，如图4B所示。每个过程测试60次。基于式(4)的计算结果表明：豆荚与304不锈钢、PVC表面的碰撞恢复系数分别为0.433、0.374，此时的最大回弹高度平均值分别为43.2、32.2 mm。

此外，基于图4A原理进行EDEM仿真验证。仿真中将其余接触参数归0处理，设定豆荚与接触面的碰撞恢复系数在0.1~1.0区间内按0.1梯度取值。随后，以设定的恢复系数为自变量，仿真获取的最大回弹高度为因变量绘制散点图，并进行回归拟合，得出非线性方程(图5)。

Figure 5.  Coefficient of restitution and maximum rebound height fit curve

将试验所得系数代入拟合方程中，求得豆荚与PVC、304不锈钢的理论反弹高度分别为32.47和44.42 mm。该预测值与实测值相对误差分别为0.84%、2.82%。由此最终标定：豆荚与PVC、304不锈钢材料的碰撞恢复系数分别为0.374、0.433。豆荚间碰撞恢复系数采用与豆荚-304不锈钢材料碰撞恢复系数相同的方法测定，以自制豆荚板代替不锈钢碰撞板。将单个豆荚从距豆荚板230 mm处自由释放，利用Tracker软件提取碰撞轨迹及最大回弹高度，根据式(4)计算碰撞恢复系数。考虑到豆荚外形不规则，且不同释放姿态和接触位置下的碰撞形式存在差异，碰撞恢复系数具有一定离散性，因此通过多次重复试验确定其取值范围为0.250~0.400，并作为后续P-B试验的参数取值范围。

本研究采用斜面法测定静摩擦系数^[21]。豆荚在斜面上的受力如图6A所示。将豆荚置于接触面一侧并调节斜面倾角，在豆荚即将发生滑动时记录倾角值，如图6B所示。代入式(5)计算求得对应的静摩擦系数。

Figure 6.  Static friction coefficient calibration test

式(5)中：μ₀为静摩擦系数；f为静摩擦力(N)；m为豆荚质量；F_N为支持力(N)；θ为接触面与水平面夹角(°)。每组试验重复60次，得到豆荚与304不锈钢、PVC之间的静摩擦系数分别为0.505、0.463，此时的摩擦角分别为26.80°、24.85°。

在EDEM仿真中，保持已标定参数不变且其余参数归0，将静摩擦系数设定为0.1~1.0区间内以0.1为梯度取值。提取各组静摩擦系数与对应的斜面临界倾角构建散点图，并进行回归拟合，得出非线性方程。结果如图7所示。

Figure 7.  Fitting curve of static friction coefficient and angle

将实测系数代入拟合公式中，求得豆荚与PVC、304不锈钢的理论滑动倾角分别为25.85°和27.36°。对比表明：该结果与实测值相对误差分别为4.02%、2.09%，由此最终标定豆荚与PVC、304不锈钢的静摩擦系数分别为0.463、0.505。同理，测得豆荚与豆荚间的静摩擦系数为0.55~0.85。

利用固定斜面滚动法测定滚动摩擦系数^[22]，试验原理如图8A所示。假设豆荚全程仅克服滚动摩擦做功，将其从倾角为θ的斜面上静止释放，在斜面与水平的滚动距离分别为S与L_s，依据能量守恒定律建立式(6)。

Figure 8.  Rolling friction coefficient calibration test

式(6)中：μ_s为滚动摩擦系数；S、L_s分别为在斜面、水平的滚动距离(mm)。

实际试验时，斜面与水平面均铺设相同接触材质，设定斜面倾角为45°，初始释放点对应的斜面距离为300 mm以保证充足的水平滑行余量，如图8B所示。通过式(6)计算得出：豆荚与304不锈钢、PVC及豆荚自身间的滚动摩擦系数分别为0.08~0.18、0.05~0.15和0.05~0.25。

鉴于豆荚的形状不规则且不易发生滚动，碰撞反弹微弱，直接测量豆荚间的接触参数以及豆荚与PVC、304不锈钢的滚动摩擦系数的误差较大，因此，本研究以堆积角为指标，依次开展P-B试验、最陡爬坡试验及Box-Behnken响应面法分析，以标定并求解上述接触参数的最优组合。

试验所用装置分为内径90 mm，高度200 mm的PVC柱形圆筒和进口内径120 mm、出口内径90 mm、高度220 mm的304不锈钢锥形桶。试验中缓慢上提圆筒，使物料受重力滑落形成稳定堆积体。

为消除传统人工测量带来的主观误差，借助Matlab软件对原始图像统一处理：依次进行灰度化、二值化及内部填充；进而提取边界像素坐标拟合直线，通过求取斜率的反正切值计算倾角大小。经多次试验取平均值，得出豆荚在PVC、304不锈钢表面的堆积角分别为37.15°、39.15°。

在EDEM仿真中采用与物理试验一致的304不锈钢圆锥筒模型。在锥筒上方设定颗粒工厂并按实测比例生成单粒、双粒、三粒豆荚，颗粒静置稳定后，以0.05 m·s⁻¹匀速向上提拉锥筒，待物料散落堆积稳定后停止仿真，并通过Matlab软件对图像进行处理，得到仿真堆积角。

为了提高参数标定效率，借助Design-Expert 13软件开展了12组P-B筛选试验。为筛除次要因素，试验以堆积角为响应值，集中提取对该指标有显著作用的核心控制参数。各测试因素的水平编码设置如表2所示。根据前期单因素仿真试验结果，仿真并测量堆积角。P-B试验结果如表3所示。

基于试验数据开展方差分析，构建了各参数标准化效应的半正态概率分布及帕累托图像(图9)，因素A、C与D不仅显著脱离了残差对角线，并表现出突出的t值响应，确定因素A、C、D对堆积角的影响显著，因此对这3个因素进行最陡爬坡试验，而因素B影响较小，取均值0.325。

Figure 9.  Plackett-Burman experimental analysis

为进一步确定显著因素的最佳取值区间，基于前期P-B试验筛选出的因素A、C、D实施最陡爬坡试验设计，以仿真堆积角与实测结果之间的相对误差为参数的评价指标。具体试验方案与结果如表4所示。

试验结果显示：第4组相对误差最小，为1.88%，因此，以第4组参数为中心水平，相邻第3、5组为低、高水平，进一步开展Box-Behnken响应面设计，对回归模型进行显著性检验，进而寻优求解出最佳的接触参数组合^[23]。

选用第4组处理条件作为Box-Behnken响应面设计的中心0点^[24]，构建了三因素三水平的试验模型。具体的试验方案布置与实测数据见表5。

对数据进行二次多元回归拟合后，得到的仿真堆积角(θ′)的回归模型如式(7)所示：

如表6所示：因素A、C、D对堆积角影响极显著；CD、C²对堆积角影响显著；AC、AD、A²、D²对堆积角影响不显著；在模型有效性方面，回归模型整体P值低于0.01，说明响应值与自变量间的线性关系高度显著，同时，失拟检验值为0.439 2，排除了模型失真的风险。此外，决定系数(R²)为0.983 5，校正值$R_{\mathrm{adj}}^2 $=0.962 2且精密度为24.556 9，充分论证了该二次回归方程具备优异的预测性能与数据解释能力。

为进一步探明参数间的交互关系，依托Design-Expert 13软件生成了反映3项显著因子交互特征的堆积角响应曲面，如图10所示。依托Design-Expert 13的响应面分析，将仿真堆积角设定为39.15°进行了参数寻优。得出的最佳拟合参数配置为：304不锈钢与豆荚间的滚动摩擦系数0.143，豆荚间静摩擦及滚动摩擦系数分别为0.747和0.171。将上述寻优获得的结果输入EDEM软件，开展重复性验证试验，进一步确认离散元模型的有效性。3次重复试验的堆积角分别为39.67°、40.32°、39.18°，均值为39.72°，与物理试验结果(39.15°)的相对误差为1.46%，证明所建离散元模型及接触参数具有较高可靠性。

Figure 10.  Response surface of factor interactions on angle of repose

为确定豆荚与PVC之间的滚动摩擦系数，在仿真中保持其余已测参数不变，将该系数在0.05~0.15范围内以0.02为梯度取值，开展PVC圆筒提拉堆积仿真，每种工况重复3次。提取各组系数与对应的仿真堆积角构建散点图，随后进行非线性回归拟合得出方程，结果如图11所示。将实测的37.15°堆积角代入方程求解，最终确定豆荚与PVC的滚动摩擦系数为0.091。

Figure 11.  Fitting curve between rolling friction coefficient and angle of repose

为检验标定参数在真实工况下的有效性，本研究开展了台架与仿真对比试验。在EDEM中构建与实体样机一致的仿真模型，赋予前文标定的接触参数。测试工况设定为：输送倾角依次为10°、15°和20°，304不锈钢料斗初始载量为0.3 kg，输送带速度为0.2 m·s⁻¹(图12)。

Figure 12.  Experimental validation

当输送倾角为10°、15°和20°时，运输量仿真值与实测值的相对偏差依次为5.88%、4.09%和3.37% (图13)。这一结果表明：所标定参数具有较高可靠性，为相关装备研发提供了核心仿真依据。

Figure 13.  Comparison of test results

鲜食毛豆豆荚属于典型的高含水率柔性农产品，其接触特性受材料力学性质、颗粒形状及表面结构等因素共同影响。本研究标定得到的豆荚间静摩擦系数和滚动摩擦系数分别为0.747和0.171，整体高于多数规则农业颗粒。

已有研究表明，农产品静摩擦系数通常随含水率增加而增大^[25]。较高的含水率使豆荚接触时更易发生局部变形，从而增大真实接触面积；同时，豆荚表面茸毛提高了界面粗糙度，因此导致静摩擦特性的增大。豆荚间滚动摩擦系数较大，主要与其几何形状有关。豆荚呈弯曲荚果状，形状不规则且质量分布不均，使其在运动过程中易发生滑移、翻滚及姿态调整。COETZEE等^[26]指出：滚动摩擦系数可等效表征颗粒形状、表面粗糙度及未显式建模微观结构对群体运动行为的综合影响。因此，较大的滚动摩擦系数有助于更准确地描述豆荚的群体运动特性。

仿真验证试验中，当输送倾角由10°增加至20°时，仿真与台架试验的相对误差由5.88%减小至3.37%，其原因可能在于：较低倾角下，鲜食毛豆豆荚在输送带表面易发生随机滚动、滑移及局部堆积，颗粒间接触状态变化频繁，运动过程受颗粒形状、表面茸毛及局部柔性变形等因素影响较大，因此仿真结果与试验结果存在一定偏差^[27−28]。随着输送倾角增大，重力沿输送方向的作用增强，颗粒运动方向更加一致，随机碰撞和局部滞留现象减少，群体运动趋于稳定，从而减弱了未建模因素对输送过程的影响，提高了离散元模型的预测精度^{[27, 29]}。类似规律也见于甘蔗Saccharum officinarum输送及倾斜带式输送机颗粒输送研究中^{[27, 30]}，即随着输送倾角增加，颗粒流动稳定性提高，仿真结果与试验结果的一致性增强。

总体而言，本研究标定得到的离散元参数能够表征鲜食毛豆豆荚的接触特性，经堆积角试验和台架试验验证，可为鲜食毛豆豆荚低损伤输送装备的设计与优化提供可靠的参数基础。

离散元模型构建与参数量化。针对鲜食毛豆豆荚含水率高、形态不规则及表面茸毛多的复杂生物学特性，构建了豆荚离散元模型，并结合物理试验与离散元仿真测定了其本征属性和多物料接触参数，为构建离散元模型提供了底层数据支撑。

参数寻优与宏微观映射。以堆积角为宏观响应指标，综合采用P-B试验筛选、最陡爬坡及Box-Behnken响应面优化方法，建立了接触参数与物料群体运动行为之间的非线性映射关系，明确了豆荚间静摩擦系数、滚动摩擦系数及豆荚-304不锈钢滚动摩擦系数的最优组合，实现了微观接触机理向宏观运动特性的表征。

多工况动态输送工程验证。在输送倾角为10°、15°和20°工况下，仿真输送量与实测值的相对偏差分别为5.88%、4.09%和3.37%。所建立模型能够较准确表征鲜食毛豆豆荚输送过程中的动态运动规律，为后续采运装备的结构优化提供了可靠的工程仿真依据。