基于离散元法的非球形原砂颗粒堆积特性研究

屠熹远 杨伟东 王媛媛 高翔宇

摘要 在基于微滴喷射的砂型喷墨打印工艺中，砂粒铺展过程对制件的打印精度具有重要影响，故对砂堆的堆积特性进行基于离散元法的仿真研究。首先针对砂粒复杂形状和不对称性，采用了一种近似模拟8种典型砂粒形状的多球体建模方法，并给出了相应的力学模型；然后基于EDEM软件进行了不同形状颗粒随机混合生成，模拟了不同圆筒和底盘尺寸的休止角试验装置及颗粒不同形状对休止角的影响。仿真和实验结果表明休止角试验装置的尺寸选择和颗粒形状分布对颗粒堆积特性具有重要影响；当圆筒直径为2.5 mm，底盘直径为3.5 mm时休止角随装置尺寸变化最小；杆状颗粒比片状颗粒流动性好，片状颗粒所占比例增加使砂堆休止角变大。

关 键 词 喷墨打印；离散元；休止角；尺寸选择；形状分布

中图分类号 TB44     文献标志码 A

Accumulation characteristics of non-spherical sand particles based on discrete element method

TU Xiyuan， YANG Weidong， WANG Yuanyuan， GAO Xiangyu

（School of Mechanical Engineering， Hebei University of Technology， Tianjin 300401， China）

Abstract In 3D printing technology based on droplet jet， the spreading process of sand grains has an important influence on the printing accuracy of the workpiece. Firstly， for the complex shape and asymmetry of sand particles， a multi-sphere modeling method is used to approximate eight typical sand particle shapes， and the corresponding mechanical model is given. Then different kinds of particles were randomly mixed and generated based on EDEM software， and the effects of different cylinder and chassis sizes and particle shapes on the repose angle were simulated. The simulation and experimental results show that the size selection of the repose angle device and particle shape distribution have important influence on the particle accumulation characteristics. When the diameter of the cylinder is 2.5 mm and the diameter of the chassis is 3.5 mm， the repose angle of sand piles changes the least with the size of the device. It is also found that the fluidity of rod particles is better than that of flake particles ， and the increase of the proportion of flake particles in the whole ones makes the repose angle of sand pile much larger by the corresponding simulation.

Key words 3D printing; discrete element; repose angle; size selection; shape distribution

0 引言

喷墨打印技术（3DP）是增材制造技术中一项高性价比的快速成型工艺，具有非常广阔的市场前景，针对该技术的研究工作发展迅速。其中，铺粉系统影响粉层的致密度[1]，并最终使制件性能发生变化。本文以原砂颗粒作为研究对象，目前针对原砂颗粒模拟大多采用球形颗粒作为基本模型[2]，但其实际颗粒形状并非球形，颗粒呈现出不规则形状，使得模拟仿真的结果与实际实验的偏差较大。大量物理实验和模拟表明，非球形颗粒的引入对颗粒运动和堆积产生重要影响[3-4]，同时颗粒形状的变化会导致流动性的变化[5]，从而改变了颗粒的分布及平均运动速度。非球形颗粒的空间取向分布直接影响砂层的致密度和孔隙率。因此，砂粒形状的准确描述和分析砂粒形狀对颗粒堆积特性的影响，对3DP工艺具有重要的指导意义。

近些年来，离散元法（DEM）被广泛应用于颗粒力学问题，但在使用过程中通常对颗粒形状和颗粒间接触模型等方面进行简化[6]。在非球形颗粒的离散模拟过程中，颗粒之间的接触判断更为复杂，接触力学模型相对于球形颗粒也要发生变化。目前描述三维非球形颗粒形状的方法主要有[3-4]：椭球、超二次、多面体、多球体组合和复合几何元素组合。其中多球体组合方法采用球体作为基本单元，原则上能以足够精度模拟任何复杂形状，最大优势是可以直接应用为球形DEM开发的接触检测评估程序，是非球形离散元仿真中最通用的方法。Alireza等[7]采用离散元方法模拟漏斗排出颗粒的过程，并用多球体方法模拟真实颗粒形状，实验与仿真的一致性良好。Haeri等[8]在增材制造领域采用离散元方法仿真杆状颗粒的铺粉过程，用多个球体组成杆状颗粒，发现颗粒纵横比越大或棍子的平移速度越大粉床的表面质量越差。Kruggel-Emden等[9]提出了用多球体组合模拟真实非球体颗粒，并发现多球体的宏观碰撞特性与圆球体仿真相比更接近真实的情况。张程林等[10]认为砂土颗粒的形状对砂土的性质具有重要影响，提出根据颗粒加权圆度模拟其真实形状。

纵观国内外研究发现，多球体方法适用于增材制造铺粉过程的颗粒形状模拟。本文在上述研究基础上，用多球体的方法解决对实际颗粒形状的描述，其中涉及到不同球体（来源于砂粒显微镜观察的结果）的混合方式（采用等概率方法），研究分析了数值模拟中试验装置尺寸及颗粒形状分布对砂堆流动性的影响，对砂型打印机的工艺参数的选择具有一定的指导意义。

1 力学模型

本文选用软球模型作为颗粒系统基本模型，这种方法更适用于碰撞持续时间与粒子自由飛行的持续时间相当的系统。颗粒间接触力学模型如图1所示。

Hertz-Mindlin接触模型[11]认为法向力[Fnij]由作用在颗粒i和颗粒[j]上的弹性力和阻尼力组成，分别由弹簧和法向阻尼器产生。对于二维颗粒，弹性力与法向重叠量大小成正比，阻尼力与颗粒相对速度成正比，[Fnij]可表示为

式中：[δ]是法向重叠量；[kn]是颗粒的法向弹性系数；[cn]是法向阻尼系数；[vij]是颗粒[i]相对于颗粒[j]的速度；[n]为颗粒[i]的中心到颗粒[j]中心的单位矢量。

切向力[Ftij]可表示为

式中，[kt]、[ct]、[α]、[vct]分别为切向弹性系数、阻尼系数、切向位移、接触点的滑移速度。

仿真中的滚动摩擦[12]可由一作用在接触面上的力矩表征：

式中：[μ]为滚动摩擦系数；[wi]是颗粒在接触点处的单位角速度向量。

Johnson等人充分考虑了表面能对物体接触的影响，考察了两弹性球体间粘附效应。运用能量计算方法，得到了两球体的接触半径与粘附能、外载荷、球体尺寸之间的关系，提出了著名的JKR模型[13]，如公式（4）（5）所示：

式中：[a、γ、E*、P]分别是两球体接触半径，表面能密度和外载荷；[E*]和[R]分别为等效弹性模量和等效半径，由公式（6）和（7）决定：

式中，[Ei，Ej，vi，vj，Ri，Rj]分别为两球体的弹性模量、泊松比、半径。JKR模型在DEM的仿真中的可行性已经得到了广泛验证。

本文的休止角数值模拟使用EDEM软件，在原砂颗粒间的力学作用采用上述的Hertz-Mindlin with JKR模型。

2 参数设置与方法

2.1 仿真参数设置

砂堆休止角是反映原砂材料流动性的数据，休止角越大说明该颗粒材料流动性越差。为了研究喷墨打印干砂颗粒的堆积特性，进行了球形砂堆和非球形砂堆的仿真模拟，并与原砂砂堆的物理实验进行了对比。仿真参数的设定对仿真结果有重要影响，通过查阅文献[14-15]、手册和实验测定得到所有需要用到的相应参数见表1。仿真网格划分设置为[2rmin]，时间步长设置[16]满足如下公式：

式中，[rmin]，[mmin]分别为仿真中所用到的最小圆球颗粒的半径和质量。

本文采用英国马尔文公司型号为 Mastersizer 2000 的激光粒度分析仪，测试未进行筛选的焙烧砂的粒度，得到该砂样的粒径分布曲线。其平均粒径为0.173 mm，最大粒径0.450 mm，粒径分布如图2所示，符合正态分布。

2.2 多球体颗粒建模

多球体方法的一个优势是颗粒间碰撞的建模，法向接触力不会产生力矩，相对于非球体形状，球体间的接触点的判断较为简单。但形状精度与组成单元的数量息息相关[17]，随着球体数量的增加，模拟得到的非球形颗粒的形状精度越高。但由于铺砂过程中砂粒的数量较多，为了减小计算工作量和提高效率，模拟计算中近似逼近砂粒形状。

为了更好地模拟出砂粒的形状，本文利用光学显微镜对2 000颗不同形状的砂粒形貌进行了观测，根据颗粒形状的分布选取了8种典型砂粒形状，如图3所示。使用Proe软件按照颗粒形状，用多个球体组合形成颗粒，近似模拟得到相应形状的外形，如图4所示。

2.3 多球体颗粒生成混合方法

在EDEM软件中直接生成仿真所用颗粒的粒子工厂不能实现多种形状粒子在同一容器内的随机混合，通过该软件可进行二次开发的接口模块，在立方体容器中设计粒子生成算法，实现过程如图5所示。首先把立方体内部区域划分为20 mm×20 mm×30 mm的网络单元，每个网格单元为边长0.3 mm的立方体。然后选取底层左下角落处网格为第1个颗粒的生成区域，颗粒生成过程完成后循环到X轴方向的下一个网格生成，直到X轴方向的这一列网格全部生成完成颗粒。之后Y轴坐标增加1个边长长度，生成下一列粒子，直到这一层的所有网格内都生成了一类粒子。之后Z轴坐标再增加1个边长长度，生成下一层的粒子，以此类推。直到颗粒的总数量达到要求时，颗粒生成算法结束。

针对非球形颗粒，程序每循环1次生成1个颗粒，每次循环生成的颗粒形状是上述8种类型多球体颗粒的随机一种，采用等概率分配方式，概率均为1/8。每个网格内的粒子三轴坐标有0.1 mm范围的随机浮动。旋转矩阵与欧拉角实现三维空间坐标变换以保证生成颗粒的随机摆放角度。使用Box Muller方法令颗粒粒径分布服从正态分布。生成的颗粒受重力作用沿着一漏斗流入圆筒中，因EDEM软件的颗粒的接触检测算法在圆筒中形成互不交叠、无不合理孔隙的砂堆，之后去掉立方体容器和过渡漏斗。至此，圆筒仿真试验的多球体砂粒的混合准备工作完成。

3 结果与分析

3.1 圆筒试验尺寸选择

在休止角的测量中一般圆筒的直径应大于颗粒最大直径尺寸的4～5倍[18]，为了保证足够的砂粒形成稳定的休止角，圆筒高度设置均大于直径的3倍。

为了避免尺寸效应，取不同尺寸的圆筒和底盘。在仿真过程中，首先采用球形颗粒进行仿真，圆筒堆积着所需颗粒，以5 mm/s的速度提升圆筒，最终形成稳定的圆锥形砂粒堆。休止角形成过程如图6所示，球形颗粒砂堆如图7所示。

本文的休止角测量方法[19]使用Matlab软件对休止角图像进行图像处理和曲线拟合，处理过程如图8所示。首先对剖面图像进行二值化处理，然后使用Zerocross算子对二值图像进行边缘检测，并得到所有砂层表面轮廓上点的像素坐标，运用最小二乘法对这些像素竖直方向的坐标进行线性拟合，最后得到拟合直线和对应的休止角度。不同圆筒直径和底盘直径条件下球形颗粒休止角的模拟仿真结果如表2所示。

由表2可以看出，当圆筒的直径从2 mm增加到4 mm，底盘直径相应从3.5 mm增加到6 mm，此时休止角受装置尺寸影响产生的变化最小，故选圆筒直径为2.5 mm，底盘直径为3.5 mm作为测量砂粒休止角模拟的装置，采用该大小的装置进行数值模拟测得的休止角为31.72°。

将图4中多球体的砂粒通过Proe软件导入EDEM软件中，按照相同概率生成8中不同形状的多球体颗粒，采用圆筒直径为2.5 mm，底盘直径为3.5 mm作为测量砂粒休止角模拟的装置。针对非球形砂粒的仿真结果如图9所示，测量得到的多球体颗粒休止角为34.4°。

将数值模拟得到的装置尺寸按比例扩大搭建试验装置如图10所示，底盘直径为140 mm，圆筒高度为320 mm，直径为100 mm。使用直流电机带动圆筒匀速提升得到砂堆，用同样的休止角测量方法测得的结果为33.6°。实验测量的休止角与球形颗粒仿真相比，误差为5.6%。而与多球体颗粒休止角相比，误差为2.4%，其误差小于球形颗粒时的仿真模拟。

3.2 颗粒形状对流动性影响

为了进一步研究颗粒形状对颗粒流动性的影响，将图4中多球体颗粒的1号、2号、4号、5号、7号划分为杆状颗粒，3号、6号、8号划分为片状颗粒。选取圆球体颗粒、3号片状颗粒及4号杆状颗粒用于数值模拟。

不同形状颗粒混合后休止角的仿真结果如表3所示。表3中仿真结果显示，不同比例的颗粒形状组成砂堆对休止角产生重要影响。其中球形颗粒的流动性最好，两种多球体颗粒相比杆状颗粒的流动性比片状颗粒好。当采用其中两种不同形状的颗粒混合时，圆筒试验形成的砂堆休止角在这两种单一形状颗粒形成的休止角之间变化。在颗粒形状表征精度足够的前提下，多球体方法会让数值模拟的休止角结果更加趋近于实验数据，流动性更接近于真实形成的砂堆。

4 结语

本文基于离散元方法，建立了砂粒间相互作用的力学模型。通过光学显微镜观察砂粒，对8种典型砂粒形状构建了多球体颗粒物理模型，得到了如下结论。

1）基于EDEM软件采用了一种针对不同形状颗粒生成混合方法，完成了多球体颗粒组合的混合过程。

2）为避免尺寸效应，在休止角模拟仿真实验中修改圆筒和底盘的尺寸，最终采用圆筒直径为2.5 mm，底盘直径为3.5 mm作为测量砂粒休止角模拟的装置。并按测量休止角模拟装置比例完成了休止角试验装置的搭建和实验。与砂粒休止角实验数据相对比，多球体颗粒仿真模拟的休止角的误差小于球形颗粒仿真模拟的休止角。砂粒的多球体模型流动性与其球形模型相比更接近实验数据。

3）研究了颗粒形状变化对休止角产生的影响，发现球形颗粒的流动性最好。其中2种多球体颗粒相比，杆状颗粒的流動性比片状颗粒好。片状颗粒所占比例增加使砂堆休止角变大。

参考文献：

[1]    ZIEGELMEIER S，CHRISTOU P，W?LLECKE F，et al. An experimental study into the effects of bulk and flow behaviour of laser sintering polymer powders on resulting part properties[J]. Journal of Materials Processing Technology，2015，215：239-250.

[2]    张帅，单忠德，张杰. 基于离散元方法的型砂流动性仿真研究[J]. 铸造技术，2016，37（2）：288-291.

[3]    LU G，THIRD J R，M?LLER C R. Discrete element models for non-spherical particle systems：from theoretical developments to applications[J]. Chemical Engineering Science，2015，127：425-465.

[4]    COETZEE C J. Review：Calibration of the discrete element method[J]. Powder Technology，2017，310：104-142.

[5]    BERRETTA S，GHITA O，EVANS K E. Morphology of polymeric powders in laser sintering （LS）：from polyamide to new PEEK powders[J]. European Polymer Journal，2014，59：218-229.

[6]    谭援强，郑军辉，张浩，等. 基于离散元法的锥形筒仓中颗粒流体的数学模拟[J]. 过程工程学报，2015，15（6）：916-922.

[7]    KHAZENI A，MANSOURPOUR Z. Influence of non-spherical shape approximation on DEM simulation accuracy by multi-sphere method[J]. Powder Technology，2018，332：265-278.

[8]    HAERI S，WANG Y，GHITA O，et al. Discrete element simulation and experimental study of powder spreading process in additive manufacturing[J]. Powder Technology，2017，306：45-54.

[9]    KRUGGEL-EMDEN H，RICKELT S，WIRTZ S，et al. A study on the validity of the multi-sphere Discrete Element Method[J]. Powder Technology，2008，188（2）：153-165.

[10]  張程林，周小文. 砂土颗粒三维形状模拟离散元算法研究[J]. 岩土工程学报，2015，37（增刊）：115-119.

[11]  王国强，郝万军，王继新. 离散单元法及其在EDEM上的实践[M]. 西安：西北工业大学出版社，2010.

[12]  HAN Q Q，GU H，SETCHI R. Discrete element simulation of powder layer thickness in laser additive manufacturing[J]. Powder Technology，2019，352：91-102.

[13]  JOHNSON K，KENDALL K，ROBERTS A D. Surface energy and the contact of elastic solids[J]. Proceedings of the Royal Society of London A Mathematical and Physical Sciences，1971，324：301-313.

[14]  姚志雄. 三维离散元法解析砂型泥石流细观运动特征研究[J]. 地下空间与工程学报，2017，13（5）：1414-1423.

[15]  顿国强，杨永振，陈海涛，等. 基于EDEM的控位分层施肥开沟器作业过程的仿真与试验[J]. 湖南农业大学学报（自然科学版），2018，44（1）：95-100.

[16]  MEIER C，WEISSBACH R，WEINBERG J，et al. Modeling and characterization of cohesion in fine metal powders with a focus on additive manufacturing process simulations[J]. Powder Technology，2019，343：855-866.

[17]  FERELLEC J F，MCDOWELL G R. A method to model realistic particle shape and inertia in DEM[J]. Granular Matter，2010，12（5）：459-467.

[18]  贾富国，韩燕龙，刘扬，等. 稻谷颗粒物料堆积角模拟预测方法[J]. 农业工程学报，2014，30（11）：254-260.

[19]  韩燕龙，贾富国，唐玉荣，等. 颗粒滚动摩擦系数对堆积特性的影响[J]. 物理学报，2014，63（17）：173-179.