基于JKR模型的辣椒籽离散元参数标定*

徐庄威，王士林，易中懿，潘健，吕晓兰

(1. 江苏省农业科学院农业设施与装备研究所,南京市,210014; 2. 南京林业大学机械电子工程学院,南京市,210037; 3. 农业农村部长江中下游设施农业工程重点试验室,南京市,210014)

0 引言

机械化播种是农业智能化发展的一种趋势,其对种子的质量与形状有一定的要求,因此小粒不规则种子的丸化加工便成为机械化播种的前提[1]。辣椒籽扁平、轻质,在造粒过程中,辣椒籽之间以及辣椒籽与釜壁之间的相互作用非常复杂,对于其涉及的一些接触参数难以定量,从而导致数值模拟时模型不精确的问题。

针对以上问题,目前学者主要以离散元法(DEM)研究散体颗粒的力学性质[2-3],以及对颗粒参数的标定[4-5]。石林榕等[6]以安息角为响应值,调节胡麻籽模拟模型时的滚动摩擦系数,最终得到胡麻籽颗粒真实滚动的摩擦系数。邢洁洁等[7]结合EDEM软件中Hertz-Mindlin with JKR模型研究海南砖红土壤的仿真参数,同样以休止角作为优化目标,最终获得仿真参数的最佳组合。侯占峰等[8]结合物理试验和无滑移模型标定了冰草籽颗粒表面参数。Coetzee等[9]为了精确颗粒在DEM仿真时参数的精确度,利用堆积角试验标定颗粒的摩擦系数。然而,关于辣椒籽的参数报告较少,且由于辣椒籽表面有粘结力的作用,简单的无滑移模型难以揭示丸化颗粒的运动机理。

因此,本文提出一种JKR粘力模型,引入表面能参数用于表征辣椒籽表面的粘结力,通过一系列物理试验和图像处理技术得到辣椒籽的基本参数,然后运用优化设计试验Plackett-Burman对所有影响因子进行显著性筛选,建立3个显著影响参数与响应值休止角的多元回归模型,以最小相对误差对模型进行优化,得到最佳参数组合。为辣椒籽丸粒化提供参数依据,提高离散元数值模拟的准确度和可靠性。

1 材料与方法

1.1 试验材料和物性参数测定

1.1.1 基本参数测定

随机选择1 000粒辣椒籽,使用精密电子天平测量其质量,将辣椒籽等分成5组,计算其质量均值。辣椒籽的规格尺寸、大小和密度由数字游标卡尺(精度为0.01 mm)和量筒(精度为1 mL)测量,含水率使用干燥法测得。共进行10次试验,结果如表1所示。

表1 辣椒籽基本物性参数Tab. 1 Basic physical parameters of pepper seeds

1.1.2 泊松比

辣椒籽的形状不同,在长度和厚度方向上也有许多差异,泊松比很难通过传统检测方法确定。使用专业质构仪(TMS-Touch),在20 ℃的试验温度下对辣椒籽进行压缩试验,并使用静态测量法获得辣椒籽弹性模量。

为避免重复,此处仅以编号为1的辣椒籽弹性模量的测试过程为例。将1号辣椒籽水平放置在试验的测试平台上,并选择直径20 mm的探针。使用Texture Lab Pro软件将传感器范围设置为0～250 N,将探头下降速度设置15 mm/min、运行时间设置5 s,选择90%形变量,开始后探头缓慢下降,得到时间—载荷变化曲线,试验结果如图1所示。

图1 时间—载荷散点图Fig. 1 Time-load scatter diagram

由图1可知,辣椒籽的时间—载荷曲线分为初始载荷段AB、线弹性段BC、峰值段CD和下降段DE。当质构仪探头刚开始挤压辣椒籽时,由于辣椒籽颗粒表面弯曲,竖直方向形变受压力的影响较大;当处于线弹性变形阶段时辣椒籽应力与应变比值为常数,如果解除应力,形变则会消失。当在BC段时,根据胡克定律,在线弹性变形下,辣椒籽的应变与应力成正比,且应力与应变的比例常数即为辣椒籽的弹性模量。

选取10次试验数据,分别对AB线弹性变形阶段进行分析。横坐标为时间变量,纵坐标为压力指数,截取3～4 s的数据拟合,获得高度为0.6 mm、横截面积为12.25 mm2辣椒籽压力随时间影响的折线图,如图2所示。

选取10组时间—载荷曲线的斜率的平均值代入式(1)计算可得到该类辣椒籽的弹性模量E为20.98 MPa。

(1)

式中: ΔF——压力增量,N;

A——厚度方向横截面的面积,mm;

Δt——时间增加量,s;

v——探头的移动速度,mm/s;

v·Δt——位移增加量,mm;

h——原厚度,mm;

k——线弹性变形阶段时辣椒籽的压力增加量与时间增加量的比值,即压缩时间—载荷曲线图上线弹性形变阶段的斜率。

泊松比在种子物性参数中极其重要,能够体现种子一定的物理性能。通过泊松比原理对其求解,过程如下。

(2)

式中:μ——泊松比;

εx——横向应变;

εy——纵向应变;

ΔL——横向绝对变形量,mm;

L——横向原长度,mm;

ΔH——厚度方向绝对变形量,mm;

H——厚度方向原长度,mm。

采用专业质构仪(TMS-Touch)作为加载装置进行泊松比试验。用游标卡尺量取单粒辣椒籽轴向和纵向的原尺寸。试验开始时,将探头下降速度设置成30 mm/min、运行时间设置成3 s,质构仪探头开始缓慢下降,给辣椒籽施加载荷直至停止。记录10次辣椒籽的纵向变形量并求均值。根据方程式(2)计算可得辣椒籽的泊松比为0.351±0.139,由式(3)可得辣椒籽的剪切模量为7.765 MPa。

(3)

式中:G——剪切模量,Pa;

E——弹性模量,Pa。

1.2 辣椒籽接触参数测定

当辣椒籽进行一些模拟试验时所涉及的参数有:辣椒籽间、辣椒籽与钢板之间的恢复系数、静摩擦系数和滚动摩擦系数。为了精确化所涉及的参数,以物理试验为参照,后用模拟试验标定。

1.2.1 碰撞恢复系数的测定试验

辣椒种子丸粒化作业时,颗粒间以及颗粒与转釜之间会碰撞接触,因此恢复系数的大小对丸粒化作业影响较大,能够直接影响辣椒种群的运动规律。根据牛顿碰撞定律,如果两个物体材质不变,碰撞后与碰撞前相对速度的比值便是恢复系数[10]。

(4)

式中:v1′、v2′——碰撞后物体1和2的速度,m/s;

v1、v2——碰撞前物体1和2的速度,m/s。

图3 恢复系数测试原理示意图Fig. 3 Schematic diagram of restoration coefficient test principle1.落种孔 2.刻度尺 3.承载架 4.底座

(5)

式中:g——重力加速度,m/s2。

H′——种子下落高度;

h′——种子反弹高度。

根据以上原理进行试验,利用UX100型高速摄像仪捕捉辣椒籽自由落体的过程,调节拍摄距离与试验台相距100 cm,整个系统主要有高速摄像仪、计算机和辅助光源。在本文中,辣椒籽均匀平铺粘合在钢板(转釜材料)上,进行辣椒籽碰撞试验,以确定辣椒籽与钢板之间、辣椒籽间的恢复系数。辣椒籽从20 cm高的承载架上自由落下,与被测物体碰撞后反弹,并使用高速摄像仪拍摄整个辣椒籽碰撞和移动过程。

利用Matlab软件将籽与钢板碰撞过程的视频进行分析,进行视频关键帧提取和三维灰度处理,通过旋转三维灰度曲面可以得到辣椒籽所在的区域以及环境的灰度值大概范围,通过阈值分割和降噪处理得到辣椒籽的区域图,接着运用region-props函数得到图像中辣椒籽区域的面积、质心,最后使用循环语句分析所有帧得到辣椒籽的运动轨迹和速度如图4所示。

(a) 关键帧提取

(b) 三维灰度处理

(c) 曲面旋转寻找阈值

(d) 阈值分割

(e) 质心标定

(f) 种子轨迹曲线图4 辣椒籽与钢板碰撞视频分析过程Fig. 4 Video analysis process of collision between seed and steel plate

图4(f)为辣椒籽颗粒和钢板接触碰撞后的路径轨迹,h′表示辣椒籽与钢板碰撞时反弹的高度,自由下落高度为20 cm,由式(5)可得到辣椒籽的恢复系数。

按照以上过程,进行10次试验,统计辣椒种间、辣椒籽与钢板间的恢复系数,计算其平均值,可知辣椒籽之间恢复系数0.377±0.039,辣椒籽与钢板之间恢复系数0.597±0.086。

1.2.2 静摩擦系数测定试验

静摩擦系数对辣椒籽丸粒化等作业影响显著。此处选用表面附有钢板的自制测试仪,如图5所示。

图5 摩擦角测量仪Fig. 5 Friction angle measuring instrument

将重量为G1的辣椒种群稳定放置在斜面仪上(此时测试仪为水平静止),摇动斜面仪的手柄,使种群所在平面的倾角α变大,当物体静止于斜面上时,可以得到

G1sinα=γG1cosα

(6)

当倾角增大,使物体表面正好有滑动的趋势时,即

tanα=γ

(7)

式中:γ——静摩擦系数。

此刻,倾角α即为辣椒籽的摩擦角。通过摩擦角与摩擦系数之间的关系便可得到辣椒籽在刚体上的静摩擦系数[11]。

辣椒籽间静摩擦系数的测量,先使用粘性材料将辣椒籽平铺固定在平面上,然后将辣椒籽群均匀置于该平面上,缓慢旋转抬起测试仪,当有少数辣椒籽在平面上滑动时,记录此刻测试仪角度,根据对应关系计算出辣椒籽间的静摩擦系数。辣椒籽与钢板间静摩擦系数原理相同,将钢板置于测量仪面上,操作一致。共进行10次试验,并统计记录摩擦系数均值。辣椒籽间的静摩擦系数均值为0.564±0.045,辣椒籽与钢板间的静摩擦系数均值为0.325±0.016。

1.2.3 滚动摩擦系数测定试验

当多数辣椒籽在倾斜提升时开始滚动,测量仪角度的正切值则是被测对象之间的滚动摩擦系数,原理与静摩擦系数测定相同。但是,由于辣椒籽较小,形状扁平,籽面粗糙,试验过程中滚动系数变化显著。经过10次试验得辣椒籽间的平均滚动摩擦系数为0.734±0.036,辣椒籽与钢板之间的平均滚动摩擦系数为0.513±0.036。但是由于物理试验与仿真试验之间存在或多或少的误差,物理试验测定的结果在仿真试验中难以具备说服力。为此,使用休止角进行模拟和标定滚动摩擦系数[12]。

1.3 物理试验休止角测定

由于辣椒籽规格较小,形状扁平轻质,在形成自然的堆垛时,其斜边是折线,休止角难以计算。因此,在累积时,使用无底圆筒填充材料,此为注射法[13]。试验材料包括一个圆柱体和一个载荷板,圆筒由亚克力板制成,根据辣椒籽大小定义尺寸。空心圆筒内径至少是辣椒籽尺寸的4～5倍,其高度应是直径得3倍[13]。为保证种子积累稳定,亚克力空心圆筒尺寸为:d=50 mm,H=150 mm。承载板由45钢制成,规格为400 mm×400 mm×3 mm。

水平拍摄试验堆积形成的自然休止角图像,截取种子堆的单侧图像通过MATLAB图像处理,经二值化、阈值分割以及膨胀腐蚀等处理,用边界函数提取边界信息,用最小二乘法调整边界[14],得到方程的斜率kθ,处理过程如图6所示,最终休止角θ表达式如式(8)所示。

(a) 原始图像

(b) 图像灰度化

(d) 图像腐蚀膨胀

(e) 边缘检测

(f) 线性拟合图6 休止角处理过程Fig. 6 Repose angle processing

θ=arctan|kθ|

(8)

利用无底圆筒法进行堆积试验,除了辣椒籽的大小、形状、自然参数和圆筒尺寸外,圆筒的上升速度也影响堆积角,为了防止堆积体过度散开并影响后续角度的准确性,圆通上升速度应小于0.07 m/s[15],上升速度越低,堆积体越稳定,但速度过低会延长仿真时间。鉴于上述情况,对圆筒不同上升速度下的辣椒籽开展了研究。如图7所示,模拟与试验中的休止角变化趋势一致。当上升速度低于0.02 m/s时,休止角与上升速度成反比;当上升速度由0.02 m/s增加到0.03 m/s时,休止角急速减小;当上升速度在0.03 m/s后,休止角几乎不变,趋于稳定。据上述分析,为提高模拟效率,在研究注入法模型在标定辣椒籽参数方面的有效性时,圆筒上升速度应为0.03 m/s。

图7 圆筒上升速度对休止角的影响Fig. 7 Influence of cylinder rising speed on angle of repose

试验时取1 000粒;将辣椒籽装入圆筒内,然后放入水平装料板上,以0.03 m/s的速度提高圆筒,在支撑板上形成稳定的堆垛,重复5次试验,求出休止角的平均值。对得到的堆垛进行图像处理,拟合结果,计算求得方程斜率kθ,得到自然休止角为28.88°、27.05°、26.77°、28.71°、28.59°,计算得均值为28.00°,标准偏差为1.01°。

2 基于JKR模型辣椒籽离散元参数标定

2.1 EDEM休止角仿真试验

通过初步物理试验测定了辣椒籽的主要物理参数,通过EDEM软件进行模拟,利用多球体重叠球模型[16-18],重叠球模型能有效减缓“自动阻塞”的发生,更好的拟合边界[19-20]。利用三维Catia建模软件,根据辣椒籽的实际全局尺寸建立几何模型(粒子模型)。然后将籽模型转换为IGES中性格式并导入EDEM作为颗粒模板。在EDEM软件中,辣椒籽颗粒模型填充独特的球形颗粒,建立辣椒籽颗粒的模拟模型,如图8所示。

图8 辣椒籽球填充模型Fig. 8 Filling model of pepper seed ball

在EDEM中选择接触模型“Hertz-Mindlin with JKR”进行模拟试验。模拟中参考物理试验圆筒的内径和高度建模,如图9所示。

(a) 种子生成

(c) 种子堆形成图9 辣椒籽堆积仿真过程Fig. 9 Chili seed accumulation simulation process

在圆筒顶端定义虚拟Polygon颗粒生成平面,选择dynamic生成模式,效率为12 000个/s,共10 000个。根据Rayleigh时间步长规则,固定时间步长在5%～30%为佳,为考虑效率选择30%固定步长,采取Euler积分,0.02 s记录一次数据。模拟共进行10 s,颗粒生成总时间为1.5 s,网格尺寸为最小颗粒直径的1.5倍[21]。辣椒籽稳定后,圆筒以0.03 m/s升起,待支撑板上辣椒籽堆稳定后,采用相同的图像处理方法,获得边界斜率kδ,并使用式(16)计算模拟辣椒籽休止角δ。

δ=arctan|kδ|

(9)

2.2 优化试验设计及分析

2.2.1 Plackett-Burman试验设计及显著参数筛选

试验设计模拟中使用的参数对丸粒化效果有不同的影响,并非所有使用的参数都对丸粒化试验有显著影响,具有显著影响的参数可用于模拟分析。由于影响参数太多,且与响应变量相关的参数的显著性尚未确定,因此本研究使用软件Design-Expert中Plackett-Burman Design设计模块过滤选择对辣椒籽运动轨迹具有较大影响的参数[22]。为避免试验分析过于繁琐,不计辣椒籽能直接精确测量的参数外,选取了以下7个参数:辣椒籽的表面能(A′)、种子—种子恢复系数(B′)、种子—种子静摩擦系数(C′)、种子—种子滚动摩擦系数(D′)、种子—钢板恢复系数(E′)、种子—钢板静摩擦系数(F′)、种子—钢板滚动摩擦系数(G′)等作为变量,考虑误差因素,在试验中添加了四个虚拟参数H′、J′、K′和L′。因此选用N=16的Plackett-Burman设计表,并以辣椒籽的休止角作为响应值。表2中参数设置高低2个水平,即+1和-1,共在试验的四个中心点进行了16次测试,试验项目和结果如表3所示。

表2 Plackett-Burman Design因素水平Tab. 2 Factors and levels of Plackett-Burman design

表3 Plackett-Burman方案及结果Tab. 3 Scheme and results of Plackett-Burman design

对结果作方差分析,获得所选参数的显著性,如表4所示。辣椒籽—辣椒籽滚动摩擦系数D′和辣椒籽表面能A′的P<0.01,表明有极大影响;辣椒籽—钢板滚动摩擦系数G′的P<0.05,表明有影响;其余仿真试验参数的P>0.05,表明影响忽略不计。

表4 Plackett-Burman试验参数显著性分析Tab. 4 Significance analysis of Plackett-Burman test parameters

Plackett-Burman试验参数显著性分析表明辣椒籽的表面能A′、辣椒籽—辣椒籽滚动摩擦系数D′、辣椒籽—钢板滚动摩擦系数G′对休止角的影响在0.05水平达显著,P值越小影响越大,对休止角的影响顺序是D′>A′>G′。

由图10可知,D′、A′和G′因素对休止角影响在界值2.365以上,影响显著。

图10 PB Design帕累托图Fig. 10 Pareto chart of PB Design

2.2.2 最陡爬坡试验设计

为更有效地进行响应面试验,找到最大梯度下降的中心点,以误差为指标,对范围进行寻优[19],爬坡试验设计及结果如表5所示。由表5可知,相对误差呈现先减小后增大的趋势,在水平3时,准确度最高,可知水平3左右为理想范围。因此,水平3被视为中心点,水平2、4分别为响应面试验的低、高水平。其中相对误差φ的表达式如式(10)所示。

表5 最陡爬坡试验设计及结果Tab. 5 Design and results of steepest climbing test

(10)

式中:δ′——休止角测量值;

θ′——休止角目标值。

2.2.3 Box-Behnken试验设计及回归模型分析

1) Box-Behnken响应曲面试验设计。在Box-Behnken试验中,由爬坡试验的结果可知,爬坡试验中水平3为中心水平,水平2、水平4被视为设计试验的-1和+1水平。根据BBD设计原则,由于有12个因素点以及5个中心点故选择3因子3水平响应面,如表6所示,进行17次试验,并通过分析获得二次回归方程以及参数的最优组合。

表6 Box-Behnken试验设计及结果Tab. 6 Box-Behnken test design and results

2) 休止角回归模型分析。对表6的试验结果再次分析,建立完全二次模型,表7为模型方差分析表。

表7 完全二次模型分析结果Tab. 7 Complete quadratic model analysis results

模型各项方差分析表明,辣椒籽的表面能A′、辣椒籽—辣椒籽滚动摩擦系数D′、辣椒籽—钢板滚动摩擦系数G′对休止角影响极其显著,同时因素A′的曲线效应也显著,交互项D′G′的影响较为显著,其余项均不显著。

仿真试验休止角二阶回归方程为

δ′=28.098+3.153A′+2.263D′+1.955G′-

0.103A′D′+0.688A′G′-1.283D′G′-

2.293A′2-0.653D′2-0.628G′2

(11)

3) 因素效应分析。如图11所示,随着表面能A′水平增加休止角先快速增加后缓慢下降,辣椒籽间滚动摩擦系数、籽与钢板间摩擦系数水平与休止角正相关。

(b) 辣椒籽—辣椒籽滚动摩擦系数

(c) 辣椒籽—钢板滚动摩擦系数图11 单因素影响趋势图Fig. 11 Effect-tendency chart of single factor

4) 不同因素间交互作用分析。响应面分析图是由休止角和3个显著因素构成的三维曲面图(图12),显示了当任一因素(A′、D′和G′)处于零水平时,其余两个因素对休止角的影响,倾角越大,交互作用越强[23]。等高线图是曲面上相同的因素值在底面上形成的曲线(图12),等高线轮廓圆扁程度体现了两因素之间的交互作用的强弱,等高线越密集,交互作用越强[23]。

图12 响应面交互作用分析图Fig. 12 Interaction analysis diagram of response surface

由图12可知,辣椒籽—辣椒籽滚动摩擦系数D′和辣椒籽—钢板滚动摩擦系数G′都与辣椒籽的表面能A′有交互作用。在表面能A′低水平时,交互作用较强;在表面能A′高水平时,交互作用较弱。

5) 参数寻优。本研究对休止角建立二次多项式模型,利用Design-Expert根据实际休止角28.00°为指标对回归模型进行优化,优化方案唯一,如表8所示。

表8 相对误差最优解Tab. 8 Optimal solution of relative error

由表8可知,在该参数组合下得到休止角为27.3°,与试验的相对误差为2.5%,可取性(Desirability)为1,说明模型可靠。因此辣椒种子丸化加工时的表面能A′为0.31 J/m2,辣椒籽间滚动摩擦系数为0.75,辣椒籽与钢板的滚动摩擦系数为0.60。

2.2.4 仿真模拟验证

为了验证标定参数的准确度,利用辣椒籽标定后的参数进行模拟,并用辣椒籽仿真的结果图像与试验图像验证,如图13所示。反复模拟三次得到辣椒籽自然休止角为25.90°、26.98°、27.25°,平均值26.71°,与模拟休止角误差为2.16 %,与试验休止角误差为4.61%,说明辣椒籽所标定的物理参数具有可靠性及参考意义。

(a) 试验积累

(b) 模拟积累图13 辣椒籽仿真模拟与试验堆积体图像对比Fig. 13 Comparison between simulation and experimental accumulation of chili seed

3 结论

1) 首先通过试验确定辣椒籽便于测量的基本参数(总粒径、千粒重、密度、含水率、弹性模量、剪切模量和泊松比);利用UX100型高速摄像仪和自制的斜面仪测试辣椒籽种间碰撞恢复系数、静摩擦系数、滚动摩擦系数的均值分别为0.377±0.039、0.564±0.045、0.734±0.036;辣椒籽与钢板间碰撞恢复系数、静摩擦系数、滚动摩擦系数的均值分别为0.597±0.086、0.325±0.016、0.513±0.036。

2) 以辣椒籽为研究对象,基于离散元仿真软件EDEM中的“Hertz-Mindlin with JKR”粘力模型,选择注入法,进行休止角累积试验,同时对辣椒籽在接触过程中的物理性能参数进行标定。

3) 利用Plackett-Burman试验过滤出对休止角影响较大的参数(辣椒籽种间的滚动摩擦系数、静摩擦系数、JKR表面能),而后通过最陡爬坡试验寻找显著项的最优区间,并通过响应面Box-Behnken试验以最佳休止角目标值(28.00°)进行显著项交互作用分析。最后通过模型优化得到最优拟合度(辣椒籽种间滚动摩擦为0.75、辣椒籽表面能为0.31 J/m2、辣椒籽—钢板滚动摩擦系数为0.60)及拟合方程,用于预测休止角。利用试验结果得出的参数组合进行验证试验,与试验休止角相对误差为4.61%,说明了模型的可靠性。