基于PFC3D的玉米籽粒离散元参数研究

戴继光

(沈阳农业大学，辽宁 沈阳 110000)

我国农业机械水平不断提高，针对粮食物料特性也有了深入的研究。玉米是我国主要粮食作物之一，玉米堆可视为各向异性的离散体聚集而成的集合体，有限元、流体等分析方式不适用于处理离散物质大量复杂行为信息，源于分子动力学的离散元法是解决散体颗粒的内部运动和受力情况的主要研究方法[1]。目前离散元已经应用于不同领域的多方面。农业工程方面主要应用于筒仓卸料、颗粒与粉体加工等[2]。许多学者将离散元法应用于农业生产、机械等方面用于分析物料在农业生产过程的运动状态[3-7]。物料的离散元参数是研究的一个重点[8-11]，深入研究玉米颗粒之间复杂信息需要得到更贴合实际的离散元参数，基于PFC3D软件的玉米籽粒参数尚不完善，本文就解决这一问题做出研究。

于庆旭等采用EDEM软件对三七种子进行参数标定[12]，模拟堆积试验得出三七种子离散元参数，堆积试验的静止角反应是粮堆摩擦，摩擦又分为内摩擦和外摩擦。内摩擦可以用内摩擦角和静止角衡量，内摩擦角是用以计算仓储设备重要力学参数[13]，静止角试验进行参数标定也是离散元参数研究的一个重要方向[14-16]。直剪试验一直是研究土体抗剪强度、测定内摩擦角和内聚力等参数的有效途径。唐福元等利用直剪仪对不同含水量的玉米的内摩擦角进行了试验测定得到准确的内摩擦角范围[17]。室内真实直剪试验只能获得宏观力学表现，史乃伟等采用PFC3D对沙土进行了直剪试验模拟以及力链分析[18]。离散元软件PFC3D中不同的接触模型下需要标定的参数不同[19]。

本文基于直剪试验用PFC3D程序对三维玉米颗粒进行填充，着重研究了线性接触模型之下在模拟直剪试验中对玉米对内摩擦角影响较大的因子。以内摩擦角为响应值采用中心组合设计试验分析出其回归方程，通过比对真实直剪试验的数据得出较为准确的模拟参数。

1 材料与方法

1.1 直剪试验方法

玉米对的内摩擦角采用直剪仪测定法，此方法是基于Coulomb理论的一种试验方法。直剪仪主要由剪切盒和伺服装置组成，剪切盒分为上下两部分，试验过程上盒固定并且施加一个垂直压力，伺服装置在下盒施加一个水平推力使得试验完全破坏。依据莫尔理论：

τ=c+σtgφ

(1)

式中，τ为剪应力，kPa；σ为正应力，kPa；φ为摩擦角，°；c为内聚力，kPa。

假设粮食不具有内聚力，得出内摩擦角与应力关系式：

φ=arctan(τ/σ)

(2)

1.2 室内真实直剪试验

试验采用玉米品种为“辽单502”，含水量11.02%，密度1197kg·m-3。依据剪切盒与试料颗粒的比例关系选择剪切盒净空尺寸为100mm的正方体剪切盒[20]，剪切速率参考其它非粘性体快剪试验选为0.5mm·s-1,剪切过程压力应根据粮食30m深度的受力载荷为依据，取最大σ为240kPa，载荷等级应该是50kPa、100kPa、150kPa、200kPa。利用土工合成材料综合测定仪在4个不同压力下重复测定5次得到的内摩擦角的平均值见表1。

表1 室内直剪试验结果

1.3 数值直剪试验

1.3.1 种粒模型建立

常见的玉米颗粒可以分为3种，体积较大的类马齿形、体积中等的类锥形和类球形。随机选取实体试验所用的玉米颗粒100粒，玉米粒下底>5mm视为类马齿型，<5mm视为类锥型，体型较小且类似球体的视为类球型，按照上述标准分类并统计尺寸数据如表2。按照尺寸表进行三维建模并导入PFC3D软件进行填充[21]。

表2 玉米籽粒参数

图1 玉米籽粒填充模型

1.3.2 线性接触模型

线性接触模型的结构如图2所示，接触力可分为线性部分F1和阻尼部分Fd。线性部分提供线弹性、摩擦行为，阻尼部分提供粘性行为。线性力通过具有恒定刚度的线性弹簧产生(ks、kn)。线性弹簧不能满足张力，通过摩擦系数μ对剪切力施加库伦准则满足滑移条件，阻尼力由阻尼器产生。

图2 线性接触模型示意图

数值直剪试验剪切盒尺寸为100mm×100mm×100mm,在后续标定试验中使竖直方向压力为200kPa，经循环后达到给定压力状态，试样在稳定压力下剪切设定速率为0.06mm·s-1,剪切位移设定20mm,对整个过程的受力状态进行检测。剪切过程如图3所示。

图3 数值直剪试验过程

1.3.3 试验设计方案

通过PB试验确定玉米籽粒参数中显著因子，根据线性接触后模型仿真需要的相关参数如表3所示。选取8个变量，每个参数选取高、低2个水平，以内摩擦角为响应，通过12组PB试验选择3个显著的参数。

表3 仿真参数设定

经过PB试验选取的3个显著参数变量，多次试验结果选取3个因素的中心点进行三因素五水平的中心组合设计试验，对得到的试验结果进行响应面分析，得出回归方程并绘出直观的三维立体响应面图，通过计算预测出合适的3个参数。

为验证得到参数的有效性，将得到的参数再次进行数值直剪试验，并且与室内真实直剪试验结果比对两者之间的误差。

2 仿真试验结果

2.1 PB试验结果

根据PB试验的分析得出3个显著的因子是有效模量阻尼系数和种间的摩擦系数，因此在后续试验中采用三因子中心组合设计试验。模拟试验结果如表4所示，由试验结果可知，第10组试验结果相较于室内真实直剪试验的结果误差最小，因此取第10组3显著因子的值作为中心点，取玉米有效模量为4.09,阻尼系数为0.72,种间摩擦系数为0.55。后续中心组试验中其它非显著因子也保持第10组PB试验中的数据。

图4 PB试验帕累托图

表4 Placket-Burman试验设计及结果

2.2 中心组合试验结果

依据PB试验得到的中心点进行三因素五水平的中心组合设计试验。试验设计以及试验结果如表5所示。

根据试验结果，进行回归方程拟合，得出关于内摩擦角和3因子之间的二次回归方程：

(3)

利用Design-Expert软件分析该模型的方差，结果如表5，得到该模型P=0.0003，说明该模型自变量与因变量相关性极其显著，方程的拟合度极好，而拟失项P=0.0147，差异显著，必须对该方程进行手动优化。

表5 中心组合试验及结果

依据方差分析在原有二次回归方程上增加高阶项ABC、A2B、A2C后发现，该模型P<0.0001,该方程的拟失项P=0.3202，相关系数R2=0.9867校正决定系数R2Adj=0.9640，相较于优化以前有较大提升，新拟合方程：

(4)

表6 优化后方差分析

新拟合的模型A有效模量和B阻尼系数交互作用明显。且B阻尼系数与C种间摩擦系数交互作用明显。应用Design Expert软件绘制响应曲面，图5a显示有效模量与阻尼系数的效应面相对于阻尼系数有效模量的影响较大，而图5b显示阻尼系数与种间摩擦系数的效应面类似于马鞍面，两者的交互作用明显。

图5 有效模量、阻尼系数与种间摩擦系数效应面

内摩擦角的目标值为23.30°，应用分析软件对新模型预测得到最优组合：玉米剪切模量为385kPa、阻尼系数为0.62、种间摩擦系数为0.051。

图6 不同法向载荷下数值模拟试验的应力变化图

数值模拟试验和真实直剪试验的对比如表7所示，随着法向载荷的变化，室内真实直剪试验和数值模拟试验的结果误差不大，均呈现随着法向载荷增大内摩擦角逐渐减小的趋势，数值模拟直剪试验剪切应力曲线趋势相似。整体而言，试验结果充分证明响应面法标定参数的可行性，也证明了数值模拟试验的正确性。

表7 验证试验分析

3 总结

通过响应面试验筛选出在线性接触模型中对试验内摩擦角影响显著的3个因子分别为玉米弹性模量、粘性阻尼系数、玉米颗粒种间摩擦系数。

通过响应面试验标定得到的3个显著因子参数为玉米剪切模量为385kPa、阻尼系数为0.62、种间摩擦系数为0.051，不显著因子玉米与钢板之间的摩擦系数为对比之间试验的结果采用响应面法标定玉米离散元参数可行，玉米颗粒之间应不存在内聚力，而仿真结果得出存在内聚力的原因是玉米颗粒之间存在相互咬合自锁行为。