随着增材制造技术的发展，砂型3D 打印技术得到了越来越多铸造技术人员的关注。目前常用的砂型3D 打印技术主要分为基于激光烧结原理和基于微滴喷射原理的砂型3D 打印技术。前者通常是利用激光器融化砂粒表面的覆膜将砂粒黏结，后者通常是采用微滴喷头喷射黏结剂到预混了固化剂的砂粒表面，黏结剂和固化剂发生交联反应将砂粒黏结[1-3]。这些主流的砂型3D 打印技术存在的共性技术是基于粉床的砂型快速制造技术。砂型制造过程中涉及了粉末颗粒的存储、供给、铺设、黏结以及回收等工序。砂粒铺设过程中砂层颗粒的密实度和表面平整度是衡量砂粒铺设质量的重要参数[4-6]。铺设砂层的密实度越高，砂型3D 打印制造的砂型强度就会越高。砂层表面平整度越好，砂型整体均匀性越好。因而，铺砂装置的稳定运行对砂型制造过程的连续性和可靠性具有重要的作用。目前存在的铺砂装置主要分为铺粉辊式、刮刀式、移动料斗式以及他们的组合。众多学者在塑料、尼龙粉末等高分子材料以及金属粉末的增材制造方面关于铺粉工艺和装置的研究较多[7-10]，但在砂型3D 打印方面，铺砂工艺基础试验研究较少。
近年来，离散单元法作为一种新颖的数值计算方法，被广泛用于求解与分析复杂离散系统的运动规律与力学特性。离散单元法将系统看作是由一个个离散的个体组成，个体之间存在着相互运动、接触以及能量的联系。Cundall 等[11]最初建立离散元单元模型是用来求解岩石力学领域的问题。现在，越来越多的研究人员将离散单元法应用到土壤、食品、物料加工及运输等散体颗粒的研究方面[12-15]。随着研究问题的深入，为了描述黏性颗粒间接触问题，Derjaguin 等[16-17]基于Hertz 接触模型，提出了一种考虑黏性功的接触模型。后来，国内外学者也发展了其他的接触模型，其中JKR(Johnson-Kendall-Roberts) 模型通过引入了表面能对颗粒接触的法向接触力项进行了修正，并用来模拟具有较强黏性力的颗粒间的接触问题[18-21]。
目前砂粒铺设工艺中由于砂粒预混了固化剂，其流动性变差，铺设砂层的过程中存在砂层致密度差以及破坏已成形砂型等问题，严重影响砂型3D打印制作的砂型质量。本文采用离散单元法，基于JKR 接触模型对采用刮刀式铺砂装置的铺砂过程进行了仿真研究，并开展了砂粒铺设密实度实验研究。通过铺砂实验，研究了砂粒铺设过程中固化剂添加量以及铺砂速度对成形砂层密实度的影响规律并与仿真结果进行对比，修正了仿真模型参数。在此基础上借助仿真模型研究了砂粒铺设过程中固化剂添加量、铺砂速度以及刮刀倾角对已成形砂型的影响。
1 离散单元力学模型
1.1 EDEM 动力学模型
EDEM(Engineering Discrete Element Method)软件是英国DEM-Solution 公司开发的通用CAE 软件，可以用来准确描述颗粒系统的运动以及颗粒间的相互作用力。在EDEM 软件中，颗粒运动的力学模型基于经典的牛顿第二定律。模型中将两个颗粒碰撞过程用弹簧-阻尼模型进行描述，颗粒的运动方程表示为：
式中，Mij 代表颗粒的质量；Iij 为颗粒的转动惯量；Kn，Kτ 分别为颗粒的法向和切向刚度系数；Cn，Cτ 分别为颗粒的法向和切向阻尼系数；ui-n，ui-τ 分别是颗粒速度的法向和切向分量；θ 为颗粒接触过程中转动的角度。颗粒间滑动的力学模型可以表示为：
式中，μ 代表了颗粒间的滑动摩擦系数；
代表了摩擦力的方向。此外，EDEM 在每个时间子步的计算过程中将颗粒的当前所受外力F 和外力矩T设定为常数，采用差分方法得到下一时刻新的颗粒速度值(
)N+δ 和(
i)N+δ，通过对式(3)积分得到新的颗粒位移值，通过不断循环更新颗粒的速度和位移得出颗粒在任意时刻的运动状态。
1.2 JKR 接触模型
在EDEM 软件中颗粒所受外力主要就是颗粒间接触作用力。为了描述不同属性颗粒之间的接触行为，EDEM 软件定义了硬颗粒接触和软颗粒接触两种不同的接触方式。JKR 接触模型属于软颗粒接触方式，即允许颗粒接触时存在一定量的穿透。Hertz-Mindlin with JKR Cohesion 接触模型可以有效地模拟接触区内范德华力以及具有较强黏性力的颗粒系统，其接触力学模型表达式为[22]：
式中，Fij-n,Fij-τ 分别为颗粒之间接触力的法向和切向分量；F
为法向接触力的弹性力项，JKR 接触模型对该弹性力项进行了修正，引入了表征颗粒间黏性力的
项，其中γ 为颗粒的表面能，a 为接触半径；
为法向接触力的阻尼力项，它与颗粒间法向相对运动速度
有关。
为切向接触力的弹性力项，是采用了Mindlin 切向力模型；
为切向接触力的阻尼力项。E*，G*，R*，M*分别为颗粒接触模型的等效弹性模量，等效剪切模量，等效半径和等效质量，其表达式如下式所示。
δn 和δτ 分别是颗粒接触的法向重叠位移和切向重叠位移，其表达式如下所示。式中，
为颗粒圆心位移向量；
为颗粒圆心连线方向的单位向量
为颗粒的切向相对运动速度。
1.3 Bonding 模型
砂型3D 打印铺砂工序中对已成形砂型采用Bonding 模型来模拟砂粒之间的接触行为。Bonding模型可以用来描述砂粒之间的一种固结状态[23-24]。在EDEM 软件中，采用Hertz-Middlin with Bonding模型，即一种不连续系统的三维模型计算程序(PFC3D)，对砂粒黏结接触行为进行模拟，接触力由法向接触力和切向接触力组成。法向接触力在颗粒动力学模型基础上增加了法向黏结力项
以及接触判断条件
A。切向接触力增加了切向黏结力项，
δτ，以及Mohr-Coulomb准则，
，μ=min(μi,μj)。
PFC3D 模型描述了两个球形颗粒接触且在接触区内存在微小的重叠，形成了一个圆形的接触面，接触面半径为aB。在一个时间子步内，球形颗粒在接触面上产生的接触力可以表示为：
式中，颗粒接触面面积A 和转动惯量J 可以表示为：
当黏结接触的颗粒所受的外力或外力矩超过设定的临界值时，颗粒之间会发生黏结破裂。颗粒间黏结的临界状态可表示为：
2 砂型打印过程的砂粒铺设仿真模型
砂粒铺设过程的仿真研究是在宏观空间里研究微观颗粒的运动，涉及到空间尺度的转变。为了描述更加真实的颗粒运动，本模型根据砂粒的真实几何尺度进行仿真研究。砂粒铺设过程的仿真模型主要包括了砂粒和铺砂装置的几何模型与物理模型，本文在EDEM 软件中建立仿真模型并进行计算分析。
2.1 砂粒仿真模型
针对铸造常用的石英砂颗粒开展砂粒铺设仿真试验研究。采用光学显微镜(OLYMPUS BX51M)观察了石英砂颗粒的颗粒形貌，如图1a 所示。为了使模型更加准确，考虑石英砂颗粒多为类似三角形和椭球形颗粒，通过颗粒长径比将颗粒分为近似三角形颗粒和椭球形颗粒，并且在软件中建立了两种不同形状颗粒的集合来模拟石英砂颗粒，如图1b所示。
图1 砂粒几何模型的建立：(a)砂粒形貌，(b)砂粒模型
Fig.1 The silica sand particles:(a)microscope figure,(b)particle model in EDEM
复杂形貌颗粒的建立方法是采用多球模型，即用多个球组合形成复杂形貌颗粒的近似模型。多球模型的各个球之间是刚性接触，内部不会产生作用力。此外，采用筛分法测量了石英砂颗粒的粒径分布，发现其粒径尺寸在70 目左右，并且近似正态分布。因而在EDEM 软件中建立了具有相同正态分布的颗粒集合，如图2 所示。
图2 砂粒粒径分布模型
Fig.2 Size distribution of the particle model
砂粒铺设过程中颗粒之间会发生碰撞、挤压以及摩擦。由于砂粒要预先添加固化剂，因而采用JKR 模型模拟颗粒之间的黏性力，试验研究了不同固化剂添加量下的砂粒铺设过程，早期的砂粒剪切试验表明不同固化剂添加量使得砂粒之间摩擦力不同[25]，砂粒剪切试验采用剪切测试仪(EDJ-I)测量砂粒的有效屈服轨迹，得到不同固化剂添加量下的砂粒间摩擦系数μ1。本模型通过定义不同的砂粒间摩擦系数μ1 来模拟固化剂添加量对砂粒铺设过程的影响，具体的砂粒模型参数如表1 所示。
表1 不同固化剂含量的砂粒物理模型参数
Tab.1 Sand particle model parameters of different contents of curing agent
试验样固化剂含量本 #1 #2 #3 #4 #5，w/% 0.00 0.25 0.50 0.75 1.00 μ1 0.708 0.684 0.662 0.635 0.622
2.2 砂粒铺设过程仿真模型
铺砂装置是采用刮刀式设计，材料为不锈钢。砂粒铺设过程中刮刀水平移动推动砂粒进行铺砂。铺砂过程是砂粒与刮刀、砂粒与砂粒之间相互运动的过程。试验装置模型如图3 所示，图中砂型代表已成形砂型，通过此模型可以获得铺砂过程中各工艺参数对已成形砂型的影响，为进一步优化工艺参数提供参考。此外，为了研究刮刀形状对砂粒铺设过程的影响，试验分别设计了0°、+30°以及-30°等不同倾角的刮刀。
图3 砂粒铺设过程的数值模型
Fig.3 Schematic of the coating process
铺砂工艺中为了提高砂型打印效率，砂粒铺设的速度是一个重要的工艺参数。为了研究铺砂速度对砂粒铺设过程的影响，试验设计了不同的刮刀移动速度，如表2 所示。
表2 砂粒铺设过程中铺砂速度
Tab.2 The speed of sand-paving in the simulation test
试验样本 #1 #2 #3 #4 #5铺砂速度/(mm·s-1) 50 120 150 20 80
2.3 仿真模型控制参数的设置
仿真模型控制参数的设置主要包括时间步长的设置和网格尺寸的设置。时间步长是指模型计算过程中两个子步之间的时间差值。DEM 程序在每一个时间子步内更新颗粒的运动状态，时间子步过大会造成计算结果失真，时间子步过小会消耗大量的计算时间。许多研究人员提出了不同的临界时间子步模型[26]，如式(11)所示。
式中，α 为常数；m 为颗粒质量；K 为颗粒等效刚度。研究发现，α 与颗粒的接触状态有关，颗粒堆积越密集α 越小，临界时间子步越小。由于颗粒碰撞过程中大部分能量是通过瑞利波消耗的，因而可以通过瑞利波沿颗粒表面传播速度来确定临界时间子步。瑞利波速度表达式如(12)所示。
临界时间步长应小于瑞利波传递颗粒半球面所需时间，即
。本模型的Δtcrit 为1.983×10-6 s，因而时间子步Δt 为0.793×10-6 s。颗粒之间的接触判断是基于网格单元法进行的，即接触只可能发生在同一网格内部或相邻网格之间。网格内部颗粒数量过多会造成接触判断的时间过长，网格内部颗粒数量越少，接触判断越简单，但需要的网格数量就会增大会占用大量的计算内存。本模型中颗粒的最小尺寸为79.6 μm，网格尺寸设定为200 μm，网格数量为1938176。
3 砂层密实度实验研究及仿真模型校核
3.1 砂层密实度实验研究
砂层密实度是表征铺砂质量的一个关键指标，砂层密实度越高，通过3D 打印技术获得的砂型强度才会越高。借助中国机械科学研究总院集团有限公司开发的数字化砂型3D 打印精密成形机(CAMTC-SMP600)进行砂层密实度试验研究。砂层的密实度D 可用砂粒的堆积密度ρ0 与真实密度ρ的比值来表征，如式(13)所示。试验过程中选取同一位置处的砂层进行取样称量，计算其堆积密度。
固化剂添加量会直接影响砂粒的流动性，研究表明颗粒含湿量的增加会使颗粒间内聚力增大，造成其流动性严重下降[27-29]。笔者研究了石英砂铺砂过程中密实度随砂粒中固化剂添加量的变化规律，如图4a 所示。从图中可以看出，固化剂添加量的增大会造成砂层密实度严重下降，不同铺砂速度下试验得到的固化剂添加量对砂层密实度的影响规律曲线具有一致的变化趋势。铺砂速度是影响成形效率的关键工艺参数，因此试验研究了铺砂速度对石英砂砂层密实度的影响规律，如图4b 所示。从图中可知，铺砂速度对石英砂砂层密实度的影响较小，随着铺砂速度的提高，砂层密实度有减小的趋势。
图4 砂粒铺设工艺试验得到的砂层密实度变化规律：(a)固化剂含量，(b)铺砂速度
Fig.4 The change in sand layer compactness obtained by the sand laying process test:(a)content of curing agent,(b)coating velocity
3.2 砂粒铺设过程仿真模型校核
砂粒铺设过程的仿真模型中影响颗粒之间相互运动的关键参数包括了JKR 接触模型中的表面能、颗粒之间的摩擦系数μ1 以及颗粒与刮刀之间的摩擦系数μ2。图4 表明砂层密实度受固化剂添加量的影响较大，砂粒中添加不同量的固化剂使得颗粒的表面能、颗粒之间的摩擦系数μ1 以及颗粒与刮刀之间的摩擦系数μ2 都不同。根据砂粒直剪试验测得不同固化剂添加量的石英砂颗粒之间的摩擦系数μ1 以及颗粒与刮刀之间的摩擦系数μ2，如表3 所示，但是由于仿真模型和真实模型之间存在一定的偏差[25]，因而对颗粒JKR 接触模型中的表面能参数进行修正。
表3 仿真模型校核后的砂粒物理模型
Tab.3 The modified sand particle model
试验样固化剂含量本 #1 #2 #3 #4 #5，w/% 0 0.250 0.500 0.750 1.000摩擦系数1 0.708 0.684 0.662 0.635 0.622摩擦系数2 0.211 0.253 0.292 0.306 0.358表面能/(J·m-1) - 0.043 0.076 0.140 0.160
设计了不同固化剂添加量下的砂粒铺设仿真试验，铺砂速度为120 mm/s，仿真结果如图5 所示，通过修正不同固化剂添加量下的砂粒表面能参数使得仿真结果与试验结果相接近，使仿真模型更加准确。
图5 试验结果与数值计算结果对比
Fig.5 Comparison of test results and simulation results
4 砂型打印过程砂粒铺设动力学研究
砂粒的铺设过程实际上是砂粒在刮刀的作用下发生相互挤压和剪切的过程。砂粒的物理性质以及铺砂工艺参数会使砂粒在铺设过程中对砂型3D 打印前一层已成形砂型产生作用力，严重情况下会破坏已成形部分的砂型，主要的破坏形式是成形部分的砂型受力破坏，以及在刮刀的作用下发生错位等。目前，砂层间作用力很难通过试验手段测量，因此采用修正后的EDEM 仿真模型研究铺砂工艺参数对砂粒铺设过程中已成形砂型的影响。图6 为固化剂添加量S=0.25%，铺砂速度V=80 mm/s 时铺砂过程中已成形砂型的受力云图。从图中可知，铺砂过程中，已成形砂型表面会受到外力作用，当刮刀经过已成型砂型表面时，已成形砂型受力最大。提取已成形砂型表层砂粒所受法向力|FZ|和切向力|FY|随时间的演变规律，获得不同仿真试验条件下的已成形砂型所受法向力|FZ|和切向力|FY|，如图7 所示。从图中可知，在铺砂过程中，刮刀推动着新铺设的砂粒经过已成形砂型表面时，会对已成形砂型表面产生法向力|FZ|和切向力|FY|，不同的固化剂添加量S 对应的已成形砂型表面最大受力不同。
图6 S=0.25%,V=80 mm/s 时，铺砂过程中已成形砂型受力云图
Fig.6 The load snapshot of the formed sand mold during the coating process with S=0.25%,V=80 mm/s
图7 不同固化剂含量下成形砂型的受力演变规律，V=80 mm/s:(a)|FY|,(b)|FZ|
Fig.7 The external force evolution rule of the formed sand mold with different curing agent contents when V equals 80 mm/s:(a)|FY|,(b)|FZ|
4.1 固化剂添加量对成形砂型的影响
固化剂添加量是影响砂型3D 打印成形工艺的关键参数之一。砂粒预混了不同含量的固化剂后会不同程度的改变砂粒的流动特性，砂粒中固化剂含量越高，砂粒的流动性越差，砂粒之间的内聚力增强。图8 为两种不同铺砂速度下，已成形砂型最大受力随砂粒中固化剂含量的变化规律。从图中可知，随着固化剂含量的提高，已成形砂型最大受力也明显增大。当固化剂添加量S=0 时，成形砂层所受剪切力小于法向力，可知干砂颗粒之间接触摩擦力起主要作用，当砂粒中加入固化剂后，已成形砂型所受剪切力大于法向力，可知砂粒间黏性力起主要作用。不同的铺砂速度下，铺砂过程中已成形砂型最大受力随固化剂添加量变化的趋势相同。因此，为了减小砂粒铺设过程对已成形砂型的影响，应尽量减少固化剂添加量。结合砂型3D 打印成形工艺需求，砂粒中固化剂添加量S 应控制在0.5%以下。
图8 不同固化剂含量对成形砂型最大受力的影响
Fig.8 Influence of the curing agent content on the maximum external force of the formed sand mold
4.2 铺砂速度对成形砂型的影响
铺砂速度是影响砂型3D 打印效率的一个至关重要的参数，同时铺砂速度也严重影响铺设砂层的质量。所谓铺砂速度是指砂层铺设过程中刮刀移动的速度。图9 为没有添加固化剂的干砂颗粒和固化剂添加量S=0.5%的湿砂颗粒在刮刀倾角为θ=0°情况下，铺砂速度对已成型砂型最大受力的影响规律。从图中可知，砂粒铺设过程中刮刀移动速度越快，已成形砂型所受外力最大值越大。此外，干砂颗粒在铺设过程中铺砂速度对已成形砂型最大受力影响相比湿砂颗粒的影响小。
图9 不同的铺砂速度对已成形砂型最大受力的影响
Fig.9 The influence of different coating velocities on the maximum external force of the formed sand mold
4.3 刮刀形状对成形砂型的影响
铺砂过程中砂粒是在刮刀的作用下发生相互挤压与剪切，刮刀与砂粒接触的倾角变化会使得砂粒在铺设过程中所受外力发生改变，因而也会在砂粒铺设过程中对已成形砂型造成不同程度的影响。图10a～b 为砂粒铺设过程中，刮刀倾角a 分别为-30°、0、+30°时砂粒中固化剂含量对已成形砂型切向和法向最大受力的影响规律。图10a 为砂粒铺设对已成形砂型切向最大受力的影响曲线，从图中可以看出刮刀倾角a 为-30°时，刮刀对已成形砂型的最大切向力影响最小，刮刀倾角θ 为+30°时，刮刀对成形砂型的最大切向力影响最大。图10b 为砂粒铺设对已成形砂型法向最大受力的影响曲线，从图中可以看出刮刀倾角a 为-30°时，刮刀对已成形砂型的最大法向力影响最小，刮刀倾角a 为+30°时，刮刀对成形砂型的最大法向力影响最大。研究表明，-30°、0、+30°倾角的刮刀在铺设砂粒时分别对应了3 种不同的砂粒剪切状态，即剥离剪切、纯剪切以及挤压剪切。剥离剪切状态使得砂粒在剪切过程中会产生垂直砂层表面向上的剥离力，在砂粒铺设过程中这种剥离力会在一定程度上减小外力对成形砂型的法向作用力，因而对成形砂型产生较小的切向摩擦力；相反挤压剪切状态使得砂粒在剪切过程中会产生垂直砂层表面向下的挤压力，在砂粒铺设过程中这种挤压力会在一定程度上增大外力对成形砂型的法向作用力，因而对成形砂型产生较大的切向摩擦力。纯剪切状态介于剥离剪切状态和挤压剪切状态之间的临界状态，从图10 中可以看出在纯剪切状态下和在剥离剪切状态下砂粒铺设过程对成形砂型最大受力的影响相差不大，这进一步表明纯剪切状态下砂粒之间相互作用更多的是处于剥离剪切状态。
图10 不同刮刀倾角下固化剂含量对成形砂型最大受力的影响，V=150 mm/s:(a)|FY|max,(b)|FZ|max
Fig.10 The influence of the curing agent content on the maximum external force of the formed sand mold at different dip angles when V equals 150 mm/s:(a)|FY|max,(b)|FZ|max
5 结论
(1)砂粒铺设工艺试验结果表明，随着固化剂添加量的增加，铺设砂层的密实度大幅度降低。铺砂速度对砂层密实度的影响较小。
(2)JKR 接触模型的表面能参数可以用来模拟砂粒之间的黏性力行为，通过将砂层密实度试验结果与仿真结果进行对比，修正了砂型3D 打印砂粒铺设过程的数学模型。
(3)砂粒中固化剂添加量以及铺砂速度等工艺参数都会对已成形砂型产生重要影响，固化剂添加量的影响作用更大，结合砂型3D 打印成形工艺的要求，砂粒中固化剂添加量S 应小于0.5%，此外，具有-30°倾角结构的刮刀在砂粒铺设过程中对已成形砂型的影响最小。
[1] DIMITROV D, SCHREVE K, DE BEER N.Advances in three dimensional printing-state of the art and future perspectives[J].Rapid Prototyping Journal,2006,12(3):136-147.
[2] GÜNTHER D, HEYMEL B, GÜNTHER J F, et al.Continuous 3D-printing for additive manufacturing[J].Rapid Prototyping Journal,2014,20(4):320-327.
[3] 朱佩兰，徐志锋，余欢，等.无模精密砂型快速铸造技术研究进展[J].特种铸造及有色合金，2013,33(2):136-140.ZHU P L,XU Z F,YU H,et al.Process of rapid technology without mold in precision sand casting[J].Special Casting&Nonferrous Alloys,2013,33(2):136-140.
[4] SHI Y, LI Z, HUANG S, et al.Effect of the polymer materials on the quality of selective laser sintering parts[J].Proceedings of the Institution of Mechanical Engineers, Part L: Journal of Materials:Design and Applications,2004,218(3):247-252.
[5] KRUTH J P,LENY G,SCHNIDEL R,et al.Consolidation of polymer powders by selective sintering[A].3rd International Conference PMI2008[C].Ghent:University College Ghent.1-16.
[6] VAN HOOREWEDER B,DE CONINCK F,MOENS D,et al.Microstructural characterisation of SLS-PA12 specimens under dynamic tension/compression excitation[J].Polymer Testing, 2010,29(3):319-326.
[7] 晏耐生，林峰，齐海波，等.电子束选区熔化技术中可控振动落粉铺粉系统的研究[J].中国机械工程，2010,21(19):2379-2382,2389.YAN N S, LIN F, QI H B, et al.Study on controllable vibration powder spreading system in electron beam selective melting[J].China Mechanical Engineering,2010,21(19):2379-2382,2389.
[8] MINDT H W, MEGAHED M, LAVERY N P, et al.Powder bed layer characteristics: The overseen first-order process input[J].Metallurgical and Materials Transactions A,2016,47:3811-3822.
[9] YUU S,UMEKAGE T,KAWASAKI Y.Numerical simulation for the friction mechanism of a powder bed using two-dimensional distinct element method[J].Advanced Powder Technology, 2006,17(3):345-359.
[10] PARTELI E J R, PÖSCHEL T.Particle-based simulation of powder application in additive manufacturing [J].Powder Technology,2016,288:96-102.
[11] CUNDALL P A, STRACK O D L.A distinct element model for granular assemblies[J].Géotechnique,1979,29(1):47-65.
[12] OBERMAYR M, VRETTOS C, EBERHARD P, et al.A discrete element model and its experimental validation for the prediction of draft force in cohesive soil[J].Journal of Terramechanics, 2014,53:93-104.
[13] JIANG M J, YU H S, HARRIS D.Bond rolling resistance and its effect on yielding of bonded granulates by DEM analyses[J].International Journal for Numerical and Analytical Methods in Geomechanics,2006,30(8):723-761.
[14] OWEN P J, CLEARY P W.Screw conveyor performance: Comparison of discrete element modelling with laboratory experiments[J].Progress in Computational Fluid Dynamics, an International Journal,2009,10(5-6):327-333.
[15] BOAC J M, CASADA M E, MAGHIRANG R G, et al.Material and interaction properties of selected grains and oilseeds for modeling discrete particles[J].Transactions of the ASABE, 2010, 53(4):1201-1216.
[16] DERJAGUIN B V.Untersuchungen über die reibung und adhäsion,IV[J].Kolloid-Zeitschrift,1934,69:155-164.
[17] DERJAGUIN B V,MULLER V M,TOPOROV Y P.Effect of contact deformations on the adhesion of particles[J].Journal of Colloid and Interface Science,1975,53(2):314-326.
[18] MORENO-ATANASIO R,ANTONY S J,WILLIAMS R A.Influence of interparticle interactions on the kinetics of self-assembly and mechanical strength of nanoparticulate aggregates[J].Particuology,2009,7(2):106-113.
[19] BARBER J R,CIAVARELLA M.JKR solution for an anisotropic half space[J].Journal of the Mechanics and Physics of Solids,2014,64:367-376.
[20] BORODICH F M, GALANOV B A, KEER L M, et al.The JKR-type adhesive contact problems for transversely isotropic elastic solids[J].Mechanics of Materials,2014,75:34-44.
[21] BRILLIANTOV N V, ALBERS N, SPAHN F, et al.Collision dynamics of granular particles with adhesion[J].Physical Review E,2007,76:051302.
[22] JOHNSON K L, KENDALL K, ROBERTS A D.Surface energy and the contact of elastic solids[J].Proceedings of the Royal Society A:Mathematical Physical and Engineering Sciences,1971,324(1558):301-313.
[23] POTYONDY D O, CUNDALL P A.A bonded-particle model for rock[J].International Journal of Rock Mechanics & Mining Science,2004,41:1329-1364.
[24] JIANG M H, YU H S, LEROUEIL S.A simple and efficient approach to capturing bonding effect in naturally microstructured sands by discrete element method[J].International Journal for Numerical Method in Engineering,2007,69(6):1158-1193.
[25] GUO Z,SHAN Z D,DU D,et al.Experimental investigation on the flow properties of sand granules in the process of sand mold printing[J].Rapid Prototyping Journal,2018,24(9):1599-1608.
[26] O’SULLIVAN C, BRAY J D.Selecting a suitable time step for discrete element simulations that use the central difference time integration scheme[J].Engineering Computations, 2004,21(2/3/4):278-303.
[27] CHEN P, YUAN Z L, SHEN X L, et al.Flow properties of three fuel powders[J].Particuology,2012,10(4):438-443.
[28] GANESAN V, MUTHUKUMARAPPAN K, ROSENTRATER K A.Flow properties of DDGS with varying soluble and moisture contents using jenike shear testing[J].Powder Technology,2008,187(2):130-137.
[29] EMERY E, OLIVER J, PUGSLEY T, et al.Flowability of moist pharmaceutical powders[J].Powder Technology, 2009, 189(3):409-415.