散体颗粒物质广泛存在于自然界中，在采矿、冶金、食品、能源、化工和农业等多个领域对相关散体颗粒物料的研究屡见不鲜。在工、农业机械化作业过程中，散体颗粒物料可能是其工作对象，也可能是其承载介质^[1]。散体颗粒物料与机械装备的接触关系和动态响应特性直接影响到机械装备的工作效率和作业性能，因而对有关散体颗粒的研究得到了工、农业领域研究学者的广泛关注。

随着计算机及光学技术的发展，一种基于三维激光扫描的点云数据获取技术在逆向工程领域得到应用与发展^[2]。李骅等^[3]及刘彩铃等^[4]通过三维激光扫描方法建立稻谷等非规则颗粒的三维实体模型并与实际稻谷尺寸进行对比，满足模型精度要求。Cundall等^[5]于1979年提出了离散单元法来计算稠密颗粒集的运动行为，并被许多学者研究和使用。在离散单元法数值计算过程中大部分学者选择球形颗粒代替颗粒物料。实际问题中，形状简单的球形颗粒无法真实模拟如矿石、玉米、花生、水稻等具有非规则外形颗粒间的碰撞、翻滚等特性^{[1, 6]}。随后Abbaspour-Fard^[7]和Kruggel-Emden等^[8]提出一种多球聚合体方法(multi-sphere method，MSM)对表面较圆滑的农业物料(如：水果、蔬菜等)进行离散单元法颗粒模型建立，但该方法填充球数较少，颗粒模型与实际轮廓相差较大。近年来，许多学者又对非规则形状离散物料的颗粒建模方法进行了研究。于亚军等^[6]设计出基于颗粒聚合体的玉米果穗离散元颗粒建模软件。Markauskas等^[9-10]通过聚合多球体的方法构建了玉米籽粒颗粒模型并通过连续体的方法来获取颗粒簇的自然离散分布情况；同时，又使用分层多球法构建椭球颗粒模型。

随着离散单元理论的完善与颗粒离散物料处理要求的提高，许多学者对颗粒离散元建模精度进行研究。Parafiniuk等^[11]在DEM软件中运用多球法构建了8个不同形状颗粒模型(具有相同的体积，高宽比为1.0~2.5)，研究不同形状对扁长颗粒集单轴压缩效果的影响。Ka Acˇ1 ianauskas等^[12]应用多球近似法和统计评价法来描述不规则准球面颗粒模型，并用多球随机重叠法来构造颗粒模型表面。Ferellec等^[13]在DEM中用多球重叠法建立非规则形状物料颗粒模型，研究球体重叠对非均匀密度颗粒惯量扭矩的影响；并通过修改密度参数来提高数值模拟精度。

综上所述，大部分研究学者基于CAD图形方法采用手动填充方式完成对非规则形状颗粒物料的颗粒建模。虽然该方法可以建立出形状比较复杂的颗粒模型，但是不能保证颗粒模型的精确度且花费时间较长。农业颗粒物料仍然是非规则形状颗粒的代表，本文以‘中早39号’早稻籽粒为研究对象，通过3D激光扫描技术准确重构水稻籽粒三维CAD模型，并对比分析边界包络手动填充和工程离散单元法(engineering discrete element method, EDEM)中API自动填充2种颗粒建模方法的建模效率与颗粒模型精度，通过仿真与试验对比籽粒自然堆积情况来检验颗粒建模方法的可靠性，以期为应用离散单元法研究非规则形状颗粒的动态仿真提供一种高效精确的建模方法。

1 材料与方法 1.1 离散单元法的接触模型

离散单元法把分析对象看成充分多的离散介质，每个颗粒或块体为1个单元，根据牛顿第二定律求出全过程中的每一时刻各颗粒间的相互作用，计算接触力和每个离散单元的运动状态。引入时间步长的概念，递归地求解离散单元的运动过程，交替反复运算，实现颗粒运动情况的预测。根据处理问题的不同，一般有硬球模型和软球模型2种^[14]。本文使用软球模型，将颗粒与颗粒、颗粒与边界的接触采用振动运动方程进行模拟。将接触模型表示成振动模型，振动运动的法向运动和切向运动、颗粒与颗粒间、颗粒与边界之间的滑动如图 1所示。

将颗粒接触过程中的振动运动进行法向和切向分解，法向振动运动方程为：

(1)

颗粒接触过程的切向振动运动表现为切向滑动与颗粒的滚动：

(2)

(3)

式中：m_1，2为颗粒i，j的等效质量；I_1，2为颗粒i，j的等效转动惯量；s为旋转半径；μ_n、μ_s分别为颗粒的法向和切向相对位移；θ为颗粒自身的旋转角度；F_n、F_s分别为颗粒所受外力的法向分量和切向分量；M为颗粒所受的外力矩；K_n、K_s为接触模型中的法向及切向弹性系数；c_n、c_s为接触模型中的法向及切向阻尼系数。

1.2 水稻籽粒颗粒模型重构 1.2.1 水稻籽粒3D模型重构

本文采用我国南方‘中早39号’水稻籽粒为研究对象。通过非接触式激光三维扫描仪^[2]得到谷物的三维点云图数据(图 2)，对点云图进行格式转换、消除噪声、平滑、对齐、合并、插值补点等一系列数据处理^{[2, 15]}，然后导入点云图数据处理软件进行谷物三维实体模型重构，得到如图 3所示稻谷三维实体模型。

1.2.2 水稻籽粒3D模型与实际籽粒尺寸对比分析

对比籽粒3D模型与实际籽粒得到长、宽、厚和体积误差均在5%以内(表 1)，符合使用要求。

1.3 水稻籽粒颗粒建模 1.3.1 基于边界包络手动填充法颗粒建模

多球聚合法构造水稻籽粒颗粒模型就是通过设置若干个粒径不一的球体三维坐标及物理半径与接触半径使其交叠在一起形成一定形状的过程，形成的固定空间有重叠也有间隙(图 4)，通过互补的方法使其和实际模型具有近似的体积和质量。以‘中早39号’早稻籽粒为样品，通过获取的三维点云图数据及处理软件(Geomagic Studio)进行三维实体建模形成包络线并得到样品坐标信息，设置球体三维坐标、接触半径及物理半径使填充颗粒模型三轴尺寸及体积接近统计结果(球体填充示意图如图 5所示)，通过手动填充操作得到图 6所示水稻籽粒颗粒模型，籽粒填充球数为164个，耗时约25 h。

1.3.2 基于EDEM中API自动填充颗粒建模

1) 网格单元划分及网格坐标信息获取：网格单元划分是实现DEM快速填充建模的一个重要环节，网格质量的好坏直接影响颗粒模型的球体填充数量和仿真计算规模。为建立正确、合理的有限元模型，划分网格时应考虑网格的数量、疏密、阶次、质量、分界点、布局、编号等。由于水稻籽粒模型形状较为复杂且曲率变化较大，在ICEM CFD中应用其自动划分网格功能来生成自由网格^[16]。考虑到颗粒填充得到的模型与真实稻谷的相似程度和仿真中所需要的计算时间，将水稻籽粒3D模型划分7种不同数量网格进行对比分析。同时利用Fluent计算流体动力学仿真软件中用户自定义函数(uesr-ddfined function，UDF)对FLUENT进行个性化设置，通过C语言编写可获得网格坐标的“udf.h”文件，得到网格坐标信息。

2) 水稻籽粒颗粒建模：API(application programing interface)是EDEM的应用程序编程接口的简称，利用EDEM所提供的API接口，可以通过编写插件来自定义接触模型、外部耦合和声场颗粒模型^[14]。在EDEM的仿真过程中插件的调用顺序如图 7所示，每个计算流程的主函数运行前会与相应的插件连接。在本次颗粒建模试验中通过VS2010编程软件创建相应的颗粒工厂动态链接库文件(BlockFactory_x64.dll)，并调用Fluent中生成的网格坐标文件(Block_Factory_Data.txt)进行EDEM中颗粒模型的自动填充。获得的颗粒模板信息中包括：球体物理半径系数、球体接触半径系数、球体ID和球体三维坐标等。EDEM多球体填充过程中实际填充球体物理半径与接触半径等于所设置值乘以相应系数求得。多球聚合法颗粒建模过程如图 4、图 5所示，构建不同数量球的填充颗粒模型如图 8，单个水稻籽粒颗粒模型花费时间约1 h。为进一步研究自动填充颗粒建模精度问题，建立7种不同球数颗粒模型，每种模型中软球接触半径保持默认为1 mm，软球物理半径分别取1、1.1、1.2、1.3、1.4、1.5、1.6 mm等，共生成49个颗粒模型进行对比分析。

1.4 试验设计 1.4.1 自然堆积角仿真试验设计

自然堆积角是表征颗粒物料流动、摩擦等特性的宏观参数，该值测定物料的流动性级别及计算物料的内摩擦角^[17-19]。本文通过无底圆筒法在EDEM软件中模拟水稻籽粒自然堆积过程，在EDEM中API颗粒快速填充后得到稻谷颗粒模型，水稻籽粒颗粒模型长轴约7.3 mm，圆筒直径应大于稻谷长轴4~5倍，高度与直径之比为3 : 1。故取圆筒直径35 mm、长度105 mm进行数值模拟。圆筒上升速度^[20]选用0.05 m · s^-1。颗粒工厂生产颗粒约2 000颗，时间步长为3×10^-6 s。根据国家标准^[21]对水稻籽粒物理特性参数进行测定，具体模拟参数如表 2所示。

1.4.2 水稻籽粒自然堆积角验证试验设计

为验证水稻籽粒颗粒模型的精度，以真实水稻籽粒、不锈钢无底圆筒进行自然休止角验证试验。将收获后‘中早39号’稻谷随机选取2 000颗(约57.14 g)，人工去芒去杂。为最大限度符合模拟条件，将选取的稻谷含水率调至3%~6%，不锈钢圆筒规格按模拟尺寸确定，速度与仿真条件一致为0.05 m · s^-1速度牵引无底圆筒缓慢匀速上升以降低圆筒抖动、偏心而造成的误差。自然休止角测量工具为电子数显游标卡尺，精度0.01°。试验重复7次，并进行自然休止角度测量。

2 结果与分析 2.1 水稻籽粒颗粒模型尺寸分析

2种颗粒建模方法得到的籽粒颗粒模型与实际籽粒得到的长、宽、厚和体积对比结果如表 3。可以看出：手动填充颗粒模型尺寸误差均在5%以上；自动填充颗粒模型除厚度误差为5.86%外，其余尺寸误差均在5%以下，自动填充得到的籽粒颗粒模型符合离散单元仿真计算要求。

2.2 籽粒颗粒模型堆积密度分析

在多球聚合颗粒建模研究领域中，堆积结构是学者们研究的重点。它决定了仿真材料或系统的宏观性能，有效控制球体堆积结构成为颗粒建模的关键要点。堆积结构主要通过堆积对象球体和堆积形成的孔隙来描述，描述参数主要有堆积密度、经向分布函数、球体的配位数和孔隙大小分布等^[22]。水稻籽粒为非规则形状颗粒，采用不同直径球体堆积而成，故通过堆积密度来评定颗粒模型的优劣。堆积密度是堆积系统中所有球体的体积与堆积空间体积之比。堆积密度与堆积系统的孔隙率密切相关，根据Torquato的理论，可通过下面公式得出：

(4)

式中：φ为堆积密度；V(S_i)表示球S_i的体积，其值为为堆积空间Ω的体积。

将得到的49个水稻籽粒颗粒模型通过堆积密度指标进行对比分析，进行高斯拟合得到图 9和公式(5)。堆积密度越接近1，颗粒模型的精确度越高，考虑到填充球数越多数值仿真时间越长与颗粒模型堆积密度等因素，选用填充球数203、填充球体物理半径1.2 mm的颗粒模型进行水稻籽粒自然堆积角的数值模拟。

(5)

式中：ρ为颗粒模型堆积密度；erf为高斯误差函数，；x为填充球数(个)；y为填充球体物理半径(mm)。

2.3 自然堆积角仿真与试验对比分析

水稻籽粒数值模拟自然堆积角的测量采用图像处理技术，x、y表示实际像素点，不具有实际量纲：1)提取堆积图像的一半，利用MATLAB读取稻谷堆积单侧图像并进行灰度处理；2)将籽粒堆积单侧图像二值化处理并搜索图像边界；3)选取图像坡度边界，并对选取的边界进行线性拟合，得到方程斜率k，如图 10。最终自然堆积角测定公式如下：

(6)

式中：θ为籽粒数值测量堆积角(°)；k为斜率。

试验中自然堆积角的测量方法为：对堆积后形成椎体的底面直径(d)和高度(h)以不同体位进行多测量并取平均值(图 11)，水稻籽粒试验堆积角(θ′)由以下公式计算得出：

(7)

式中：θ′为测量的水稻籽粒堆积角；h为稻谷堆积高度(mm)；d为稻谷堆底直径(mm)。

对仿真和实际水稻籽粒自然堆积角试验分别重复7次，得到如表 4所示结果。数值仿真较试验结果平均值误差为0.76%，证明离散单元法快速多球填充颗粒建模方法有效可行。

3 结论

本文在对我国南方水稻籽粒物理特性进行分析研究的基础上，通过2种方法进行水稻籽粒离散单元法颗粒建模，并进行自然堆积角试验，得到以下结论：1)对我国南方‘中早39号’早稻籽粒物理特性参数进行测定，通过非接触式三维激光扫描方法重构籽粒三维实体模型，与实际籽粒对比尺寸误差均在5%以内。2)通过手动填充及自动填充2种方式对水稻籽粒进行颗粒建模得到：手动填充法得到颗粒模型填充球数164个，耗时25 h，与实际籽粒尺寸对比误差均在5%以上；自动填充法得到颗粒模型填充球数203个，耗时1 h，与实际籽粒尺寸对比除厚度误差较大外其余误差均在5%以内，自动填充法颗粒建模具有较高效率及精度。3)对于自动填充建模方法，通过高斯拟合法对不同填充参数(填充球数、球体物理半径等)与颗粒模型堆积密度之间的对应关系进行非线性拟合，并得到填充球数203个，物理半径1.2 mm时颗粒模型最优。4)水稻籽粒自然堆积角仿真结果与试验结果平均值误差为0.76%。证明基于EDEM中API多球聚合自动填充方法重构水稻籽粒颗粒模型有效、可行。该方法减少了颗粒建模的工作时间，提高了颗粒模型精度，对于非规则形状散体物料的离散单元法研究具有一定价值。