谷壳垫料颗粒离散元模型参数标定

陶志影，孙家伟，韩志恒，刘 橙，张 林，叶鼎承，何金成

（1.福建农林大学 机电工程学院，福州350002；2.福建省农业科学院畜牧兽医研究所，福州350013）

异位发酵床技术是我国7种主要的畜禽粪污处理模式之一[1]，对发酵床垫料进行翻抛作业是异位发酵床运行必不可少的过程。建立异位发酵床垫料颗粒离散元仿真模型，对翻抛刀具的设计有着重要意义。异位发酵床中垫料的原材料主要是锯末、谷壳、玉米秸秆芯和麦麸等[2]，本研究选取谷壳作为发酵床垫料，对谷壳部分物理特性进行测定，并进行谷壳颗粒离散元模型参数标定，建立异位发酵床垫料颗粒模型，为后期利用离散元法建立耙刀与垫料床的模型，对翻抛刀具进行仿真设计研究提供理论基础。

离散元仿真模型的确立需要泊松比、密度和剪切模量，基本接触参数包括碰撞恢复系数、静摩擦系数和滚动摩擦系数。在EDEM 中建立的颗粒模型与真实颗粒外形存在差异和接触特征复杂，使得接触参数与真实值存在误差，需要对仿真模型的接触参数重新标定[3]。目前大多学者对仿真参数的标定是通过堆积角模拟试验的方法进行离散元模型接触参数标定[4-6]。刘文政等[7]根据马铃薯的物性参数，建立了微型马铃薯离散元参数获取模型；向伟等[8]基于土壤堆积试验，构建了南方黏壤土精准的离散元模型；王云霞等[9]对玉米种子颗粒之间的静摩擦系数和滚动摩擦系数进行标定，提高了离散元仿真结果的准确性；贺一鸣等[10]对南方地区黏弹性壤土颗粒的离散元接触模型的本征物理参数和接触力学参数进行标定，建立了离散元仿真的接触模型。

本研究主要对谷壳三轴尺寸、千粒重、含水率、堆积密度、谷壳之间摩擦系数、碰撞恢复系数、谷壳与钢板之间的摩擦系数、碰撞恢复系数及谷壳颗粒自然堆积角等物理特性进行测定研究。运用EDEM 软件来对谷壳进行离散元仿真建模，选用Hertz-Mindlin (no slip)模型作为谷壳颗粒之间的接触模型[11-12]，以堆积角作为评价指标，通过单因素试验与三因素三水平正交试验相结合对垫料谷壳模型参数进行标定，来找出最佳参数组合，使得模拟试验中颗粒的堆积角与现实中颗粒的自然堆积角一样，从而建立准确的离散元仿真模型。

1 材料与方法

1.1 材料及仪器

试验样品选用安徽淮南老刘农产品直销创业店的谷壳。试验仪器有电子游标卡尺(上海美耐特实业有限公司制造，分辨率为0.01mm，量程为0～150mm)、电子天平（上海精科天美仪器有限公司生产的，型号为JA5003B，精度0.001g，量程为0～500g）、电热鼓风干燥箱中（上海贺德实验设备有限公司生产，型号为DGH-9140A）、直剪仪（南京土壤仪器有限公司制造，ZJ 型）、双划杆刻度量角器（杭州阳鑫工五金具厂制造，型号为EJDG-300）、钢质圆筒和钢制圆台。

1.2 垫料谷壳颗粒物理参数的测定

1.2.1 垫料谷壳颗粒三轴尺寸与千粒重的测定 为了表示物料全部颗粒的尺寸，随机从100kg 样本中抽取50粒谷壳，用电子游标卡尺测量每个谷壳的基本尺寸：长度L(mm）、宽度W(mm）、厚度T(mm），长是指平面投影图中最大的尺寸，宽是指垂直于长度方向的最大尺寸，厚是指垂直于长、宽方向的直线尺寸，计算谷壳的三轴算术平均粒径、几何平均粒径和球度[13]。测量数据结果为:长度（L）的均值为9.37mm、中位数9.37mm，有96%分布在8.3～10.6mm 之间；宽度（W）的均值为2.38mm、中位数2.42mm，有90% 分布在1.9～3.3mm 之间；厚度（T）的均值为1.60mm，中位数为1.61mm，有90%分布在1.2～2.1mm之间。每粒谷壳的三轴算术平均粒径均值为4.45，几何平均粒径的均值为3.28，球度的均值为0.35。从100kg 样本中随机数出5 组1000 粒的种子，分别放在电子天平上称量，取均值，得出谷壳在该含水率下的千粒重。测得5 组千粒重的均值为3.76g，每粒谷壳的重量约为0.00376g。

1.2.2 垫料谷壳颗粒含水率的测定 本研究采用常压恒温烘干法，结合国标物料湿度（含水率）的测定标准（JB/T 9014.5-1999）测量垫料谷壳颗粒含水率。从样本中称取试样1000g，置于预先恒重的称量皿中，将其放入电热鼓风干燥箱中，温度调至120℃烘干2h，取出称量皿，再盖上皿盖，冷却至室温称量，反复烘30min，再冷却称量，直至两次称量的质量之差不大于1g。进行5次测量，测得谷壳含水率的均值为14.48%。

1.2.3 垫料谷壳颗粒堆积密度的测定 堆积密度的测量填充方式有中心流法、压实法和雨流法3 种[14]。本研究选用中心流法作为测试方法，通过管道引流将颗粒填充至量筒内部，量筒带有刻线，每次的体积控制在0.001m3，放在电子天平上称量，重复5次，5次测量的堆积密度均值为103.08kg·m-3。

1.3 垫料谷壳颗粒力学参数的测定

1.3.1 垫料谷壳颗粒之间的摩擦系数测定 本研究的垫料成分主要为谷壳，谷壳的力学参数关系到模型的正确与否，主要包括内摩擦角、内聚力、剪切模量及泊松比等。借鉴土壤的方法来测量谷壳的抗剪强度[15]。土壤抗剪强度常用的测定方法主要有三轴压缩试验及直接剪切试验[16]。本研究选择直接剪切试验的方法来测量谷壳的内摩擦角、内聚力及内摩擦系数，根据剪切过程垂直载荷与抗剪强度的关系，计算得出土壤内聚力、内摩擦角及内摩擦系数等参数。利用直剪仪测量谷壳的内摩擦角、内聚力及内摩擦系数时，试样的某一面上的抗剪强度就是该面两侧的试样发生滑动时所受到的最大阻力，该阻力主要由内摩擦角和内聚力所组成，可以近似地用库伦公式为：

式中：τ为抗剪强度(kPa)；σ为承受的垂直压力(kPa)；tgφ为内摩擦角(°)；c为内聚力(kPa)。

初始含水率的谷壳分别在垂直压力为50，100，150，200 kPa 作用下各进行8 次重复试验，剪切盒的上部分以0.8mm·mim-1的速度相对下部分匀速移动。相应的垂直压力的抗剪强度值取平均值，根据得到的数据绘制成图（图1）。

图1 垂直载荷与抗剪强度的关系曲线Figure 1 Relation curve between vertical load and shear strength

根据以上数据，求得拟合直线方程为y= 0.3094x+ 1.45，R²为0.9984，接近1，说明趋势线拟合程度可靠，可知谷壳的内摩擦系数约为0.31，内聚力约为1.45kPa。

1.3.2 垫料谷壳颗粒与铁板的摩擦系数测定 斜面滑动法是测量静摩擦系数较为常用的方法[7]。谷壳在斜面上受到自身的重力、斜面对它的支撑力和摩擦力，当斜面缓缓升起的时候，斜面与地面存在一个倾斜角β；当倾斜角小于滑动临界角时，谷壳始终静止在斜面上，随着斜面不断缓慢的升起，倾斜角也随之不断增大，直至超过滑动临界角，谷壳将沿斜面下滑。静摩擦系数μ和倾角β的关系为：μ= tanβ。

斜面仪如图2，试验时，钢板处于水平位置，将谷壳放在离钢板上端边缘5cm 处的中心点，通过摇动手柄使钢板匀速缓慢的升起，当谷壳开始微微滑动时，停止摇动手柄并保持钢板与水平面的夹角不变，通过游标卡尺测量夹角β，并计算谷壳与钢板的静摩擦系数，随机抽取20 粒谷壳作为测试样本，每个谷壳按同样的方法测试3次，将3次的数值取均值，测量结果如表1，20个样本的谷壳与钢板的静摩擦系数的均值为0.57。

图2 自制斜面仪Figure 2 Self-made inclined plane meter

表1 摩擦系数测量结果Table 1 Friction coefficient measurement result

1.3.3 垫料谷壳颗粒与钢板的碰撞恢复系数测定 利用谷壳颗粒的碰撞弹跳试验对谷壳与钢板的碰撞恢复系数进行测量[7,17]，具体测量方法如图3，图中A、B、C 表示不同的平面，A 为谷壳的下落起始面，B 为钢板所在位置平面，C 为谷壳与钢板碰撞后回弹的最高点所在的平面。谷壳与钢板碰撞恢复系数计算公式为:

图3 恢复系数测试示意图Figure 3 Schematic diagram of recovery coefficient Test

式中：v1为颗粒碰撞前的接触速度；v2为颗粒碰撞后的法向方向的分离速度；H为颗粒下落的高度（200mm）；h为颗粒碰撞后回弹的最高点与钢板的垂直距离。

共做10组试验，结果见表2，e1取值区间在0.18～0.20之间，取平均值为0.19。

表2 谷壳与钢板碰撞恢复系数测量结果Table 2 Measurement result of collision recovery coefficient between rice husk and steel plate

1.3.4 垫料谷壳颗粒之间的碰撞恢复系数测定 谷壳之间的碰撞恢复系数测量的方法和谷壳与钢板的碰撞恢复系数测量的方法类似，在钢板表面铺上一层均匀的谷壳，试验时首先通过镊子夹取谷壳放置坐标纸的标记处，然后打开摄像机拍摄谷壳下落与钢板上谷壳碰撞的过程，通过拍摄的图片找出谷壳反弹的最高点并记录，每个颗粒重复测量3遍，取平均值作为该颗粒的碰撞恢复系数e2，共做10组试验，结果如表3，e2取值区间在0.07～0.09之间，均值为0.08。

表3 谷壳与谷壳碰撞恢复系数测量结果Table 3 Measurement result of collision recovery coefficient between chaff and chaff

2 垫料颗粒离散元参数标定

2.1 建立垫料颗粒模型

为了提高仿真的准确性，根据1.2.1 小节中求得谷壳颗粒的三轴尺寸，长度的平均值为9.37mm，宽度的平均值为2.38mm，厚度的平均值为1.60mm。结合EDEM 软件建模的特性，并考虑后期生成数十万颗粒床的计算量极大，采用组合颗粒模型的方法[18]，通过5个半径不同的圆球叠加出类似谷壳模型，以中间球为中心，左右两边对称（图4）。

图4 EDEM软件谷壳建模Figure 4 Modeling of chaff in EDEM Software

2.2 选择颗粒力学接触模型

本研究选用Hertz-Mindlin (no slip)，该模型在力的计算方面高效且精度高。颗粒接触模型如图5，颗粒间的法向力和切向力公式如下。

图5 颗粒接触模型Figure 5 Particle contact model

式中：E*为等效弹性模量；R*为等效粒子半径；α为法向重叠量。

式中：Ft为切向力；St为切向刚度；δ为切向重叠量。

式中：Ftd为切向阻尼力；β为系数；St为切向刚度；m为等效质量；vt为切向相对速度。

2.3 垫料颗粒堆积角物理试验

本研究采用圆筒提升的方法进行试验。试验时，通过导管往圆筒中注入一定量的谷壳，然后启动万能试验机以500mm·mim-1的速度匀速上升，等谷壳完全流出圆筒时，停止万能试验机，利用游标卡尺直接测量堆积谷壳的角度记为α（图6），同时利用直尺测量谷壳堆体的高度并通过以下公式计算出堆积角度γ。

图6 堆积角测量图Figure 6 Measurement diagram of stacking Angle

式中：γ为堆积角度；L为谷壳堆积的高度；D为圆台的直径。

堆积角测量试验重复10次，数据如表4，γ的平均值为40.97°，α的平均值为41.14°，两者的误差在5%之内，说明数据可靠，最终堆积角度取两者的平均值为41.1°。

表4 堆积角度测量数据Table 4 Stacking Angle measurement data

2.4 垫料颗粒堆积角模拟试验

2.4.1 单因数堆积角模拟试验 利用离散元仿真软件EDEM2018 版进行谷壳堆积角的模拟试验，具体仿真设置参数如表5，然后建立三维模型，模型由圆筒和半径为130mm 的圆形底板组成；在圆筒的顶端设置一块虚拟的圆板用于颗粒的生成，每秒生成350个颗粒，并以1m·s-1的速度向下运动，总共生成7000个颗粒，重力加速度设置为-9.81m·s-1；将圆筒的速度设置为25mm·s-1匀速向上提升；时间步长设置为瑞利时间步长的30%，总时间设置为23s，数据保存间隔为0.05s；仿真结束时谷壳的堆积状态如图7。对仿真结果进行处理，利用切片工具对其进行切片，切片后效果如图8，通过测量角度工具对其进行测量。

图7 堆积角堆积后图Figure 7 Stacking Angle

图8 堆积角切片图Figure 8 Stacking Angle section diagram

表5 离散元仿真所需参数Table 5 Parameters required for discrete element simulation

在其他参数不变的前提下，提升速度分别设为10，25，50mm·s-1进行仿真试验，前两者堆积的状态与堆积角度接近，而50mm·s-1堆积角度小于前两者并且颗粒落在底板上的范围更广，速度快的时候中心轮廓出现扰动现象和边缘的颗粒离中心更远。说明圆筒的抬升速度会对堆积结果产生影响，速度越小越能保证颗粒堆积后堆体的稳定性与精度。过小的速度会增加计算量造成仿真时间过长，过大的速度会造成仿真结果的不精确，本研究将圆筒的提升速度设为25mm·s-1。在其他参数不变的前提下，谷壳与谷壳之间的静摩擦系数分别设为0.21，0.31，0.41 进行仿真试验。颗粒堆体的高度依次增加，堆积角逐渐增大，说明静摩擦系数对仿真结果影响明显。在其他参数不变的前提下，谷壳与谷壳之间的滚动摩擦系数分别设为0.05，0.1，0.15 进行仿真试验。颗粒堆积状态的高度依次增加，但增加的速度比改变静摩擦系数来的缓慢，说明改变滚动摩擦系数对仿真结果存在一定的影响。在其他参数不变的前提下，谷壳的剪切模量分别设为3.75e+7，3.75e+8，3.75e+9Pa 进行仿真试验。颗粒的堆积高度随着剪切模量的增大颗粒堆积的高度缓慢变高，但在剪切模量增大的过程中，计算机的运算时间明显的增长，剪切模量越大对计算机的性能要求越高，说明剪切模量对堆积角模拟试验存在一定的影响，但在满足要求的前提下，剪切模量不易设置过大，会影响仿真效率。

2.4.2 正交设计堆积角模拟试验 由于多个因素对仿真模型都有存在影响，本研究通过设计三因素三水平的正交试验方案，来找出最佳组合，使得模拟试验中颗粒的堆积角与现实中颗粒的自然堆积角一样。通过实际测量得到颗粒的堆积角为41.1°，模拟试验结果中的颗粒堆积角与该角度值越接近，说明仿真模型可信度越高。通过单因素试验了解到对颗粒堆积角影响比较明显的参数有静摩擦系数、滚动摩擦系数和剪切模量3个参数，3个参数相互独立，无交互作用，取剪切模量A，静摩擦系数B，滚动摩擦系数C 为试验因素，D 为空集，选择因素水平如表6。

表6 正交试验因素水平Table 6 Orthogonal test factor levels

EDEM 仿真模拟试验结果如表7，m表示模拟试验的堆积角度与颗粒实际堆积角度的百分比，m值越接近百分百，说明该组的参数组合的仿真模型可信度越高。第3次试验结果为98.52%，该组模拟试验的堆积角值与实际颗粒的堆积角值最为接近，说明该组的参数组合仿真模型可信度极高，可以通过微调他们的接触参数来保证该模型与实际中颗粒堆积的模型更加的接近。为分析3 个参数对颗粒堆积角的影响，找到影响堆积角的显著性参数，对正交试验结果进行极差分析，结果如表8。由表8可知剪切模量的极差为7.04，静摩擦系数的极差为24.14，滚动摩擦系数的极差为10.71，静摩擦系数的极差大于其他两个参数，说明静摩擦系数对模拟堆积角影响最大，剪切模量对模拟堆积角影响最小。

表7 试验结果Table 7 Test result

表8 正交设计极差分析结果Table 8 Range analysis result of orthogonal design

因为只通过极差的大小来分析各个因素对试验指标影响的程度没有固定的评价标准，本研究对数据又进行了方差分析，方差分析如表9。由表9 可知，静摩擦系数B 的p=0.04＜0.05，是显著的，说明静摩擦系数B 是影响试验结果的重要因素，其他两个因数相比于静摩擦系数B是影响该试验的次要因素。

表9 正交设计方差分析结果Table 9 Analysis result of variance of orthogonal design

3 讨论与结论

本研究以谷壳作为异位发酵床垫料，对谷壳物理特性进行测定分析，选用Hertz-Mindlin (no slip)模型作为颗粒之间的接触模型，以堆积角作为评价指标，通过单因素试验与三因素三水平正交试验相结合对垫料谷壳颗粒离散元模型参数进行筛选，最终得到颗粒堆积角影响比较明显的参数有谷壳间静摩擦系数、滚动摩擦系数和剪切模量3个参数，因为谷壳的硬度不佳，滚动性不强，谷壳间静摩擦系数对堆积角影响最显著，这与相关研究[19-21]中结果一致。模拟堆积角试验与实际堆积角试验相对照分析，得到的结果更准确。

本研究主要测量垫料谷壳颗粒的物理参数与力学参数，物理参数测量包括垫料颗粒三轴尺寸为长度（L）的平均值为9.37mm、宽度（W）的平均值为2.38mm、厚度（T）的平均值为1.60mm，千粒重为3.76g，堆积密度为103.08kg·m-3，含水率为14.48%。测量力学参数包括垫料谷壳颗粒之间的摩擦系数为0.31、垫料谷壳颗粒与铁板之间的摩擦系数为0.57、垫料谷壳颗粒与钢板之间的碰撞恢复系数为0.19、垫料谷壳颗粒之间的碰撞恢复系数为0.08。对颗粒堆积角影响比较明显的参数有谷壳间的静摩擦系数、滚动摩擦系数和剪切模量3个参数，且静摩擦系数影响最大。当剪切模量为3.75e+8Pa，静摩擦系数为0.41，滚动摩擦系数为0.15时，模拟试验的堆积角值与实际谷壳颗粒的堆积角值最为接近。说明试验标定的参数精准，可为离散元仿真颗粒床建模提基本参数的标定提供理论依据，并为翻抛刀具的设计提供研究基础。