粉末压制成形作为粉末冶金工艺中的基本工序，一般情况下指的是模压成形，即将一定量的金属粉末装入到模具中，通过加压的方式使粉末颗粒在模具中成形为具有一定几何形状的压坯^[1]，其广泛应用于机械零部件制备等场合。随着计算机的发展以及计算力的提高，数值仿真技术被广泛应用于粉末压制领域^[2]。目前应用数值仿真技术研究粉末压制成形一般采用连续介质分析和非连续介质分析方法两种，分别是有限元法（finite element method，FEM）和离散元法（discrete element method，DEM）。有限元法主要用于对连续体的仿真，即将粉体视为连续介质考虑。吴斌等^[3]使用有限元软件研究了高速压制下模壁摩擦因素对铜粉和铁粉的影响。离散元法将颗粒看成单独的个体，更有利于研究粉末颗粒之间的接触并且与实际情况更为吻合，但目前将离散元法应用于岩土力学等其它领域较为广泛，而应用于粉末压制中较为有限。张璐栋等^[4]利用三维离散元软件对高速压制条件下铝粉颗粒的动态响应进行了数值模拟。莫品强等^[5]通过离散元法模拟了不同围压下的三轴排水剪切试验，采用边界去除法对不同初始各向异性的试样进行堆积试验，揭示了组构演化规律，证实了应力拱的存在，分析了应力拱对堆积体宏观特征差异的影响。

通过离散元法来研究粉末压制过程，可从微观角度对压制过程中颗粒的流动、变形等特性进行分析并对宏观现象做出相应解释，弥补有限元等连续体方法的不足。对此需要引入力链的概念：散体材料中，力经由颗粒间接触从一个颗粒传递至另一个颗粒，这一传力特征在可视化后表现为链状结构，直观地称为力链^[6]。张超等^[7]使用离散元软件模拟冲头加载过程，研究了不同冲头速度撞击金属粉末颗粒所产生的透射波以及对侧壁的压力分布，且发现了颗粒间力链的传播和分布规律。王海陆等^[8]通过建立不同粒径分布的粉末颗粒仿真模型，得到压制过程中相对密度变化规律与力链分布情况。Zhang等^[9]通过分析双轴压缩下不同应变阶段下的力链，发现力链的结构方向随着压制过程渐渐指向最大主应力方向，力链的各向异性程度也受到初始围压的影响。Zhu等^[10]通过模拟二维干颗粒双轴压缩，研究了力链在介观尺度上的几何形状以及介观结构对力链行为的影响，并通过屈曲现象指出力链的稳定性伴随着介观结构的变化而变化，而后者随着系统总体积变化。

相关力链研究多集中于二维颗粒体系状态，在力链演化的同时，颗粒系统也在逐渐开始堆积，从类流态转变为类固态，这就是阻塞现象，即粉末压制过程中可能发生的局部结构固结现象，在一定程度上将对粉体致密化过程造成影响。在研究阻塞时，首先必须确定阻塞现象发生的条件，当粉体颗粒系统达到称为“阻塞点”的特定时刻时，就会发生阻塞现象。关于阻塞点的确立，相关学者开展了广泛研究。Göncü等^[11]通过配位数在无摩擦三维颗粒中找到了阻塞点，Radhakrishnan等^[12]通过配位数发现了悬浮液中干颗粒的阻塞点，他们都在相似的无放射性条件下分析了颗粒系统中的颗粒。不同的颗粒形状和颗粒空间结构也会对“阻塞点”产生一定的影响。当阻塞发生时，颗粒的剪切力从原来的均匀分布转变为集中起来。Ward等^[13]通过收集铜粉压实过程中的相对密度、声呐管功耗和温度的现场测量数据，发现了阻塞现象，并总结了压制过程中三个不同的致密化阶段，从宏观角度分析粉末压实过程中的阻塞现象。在相关颗粒系统中，经常通过颗粒体系配位数来确定阻塞点，而对阻塞点附近力学行为的关注却较为有限。粉体颗粒体系中的力链对阻塞可能具有很大影响，力链网络可以支持颗粒系统的稳定性，从而引起颗粒阻塞。

目前，关于粉末单轴压制过程中阻塞现象的研究较少，亦缺乏对阻塞现象状态下介观力学行为分析，尤其针对三维粉体颗粒状态下的介观力学结构力链分析。本文在其他学者研究的基础上，从微观角度分析了粉末阻塞的现象。基于颗粒物质力学理论，对三维粉体颗粒压制过程进行模拟，通过总结粉体颗粒的物理运动规律，找到粉末压制过程中的阻塞点。同时进一步分析了力链随压制过程的演变和阻塞现象的变化，对力链进行了拓扑分析，发现了造成阻塞现象的介观原因，讨论了压制坯体的密度变化以及压制成形后的密度不均匀现象，为后续研究优化压坯的致密化提供了理论依据和研究方向。

1.   数值模拟

1.1   仿真模型

通过离散元软件按照级配分布生成了不同粒径的铁粉，粒径分布如表1所示^[14]。仿真模型为圆柱体模具，直径为20 mm，高度为30 mm，圆柱体模具由三面墙体组成，圆柱面墙体作为阴模模壁，顶部墙体为可做竖直运动的上模冲，底部墙体作为底面固定不动。生成颗粒体积分数为51.9%。赵婷婷和冯云田^[15]利用量纲分析方法研究得出放大定律适用于离散元接触模型，并通过筒仓侧壁压力和休止角两个算例验证了放大定律的有效性。Chen等^[16]以及Grima和Wypych^[17]使用了粒径放大进行仿真，得出的仿真结果与实验结果非常匹配，故将颗粒粒径进行放大对仿真结果不会造成影响。因此为了提高计算效率，节约计算时长，对颗粒粒径进行放大，放大倍数为10倍，粒径按照高斯分布共生成49039个三维铁粉末颗粒，颗粒位置随机生成。颗粒参数如表2所示。

表  1  粒径级配占比（体积分数）

Table  1.  Proportion of the particle size gradation

粒径 / μm	23	61	113	150
体积分数 / %	35.1	33.5	25.3	6.1

下载: 导出CSV 

| 显示表格

表  2  颗粒仿真参数

Table  2.  Simulation parameters of the particles

物理量	数值
墙体法向刚度，kn / (N·m−1)	2×1012
墙体切向刚度，ks / (N·m−1)	2×1012
墙体摩擦系数，μw	0.25
压制速度，V / (m·s−1)	1
阻尼系数，η / (N·s·m−1)	0.7
初始体积分数 / %	51.9
颗粒摩擦系数，μp	0.25
颗粒密度，ρ / (kg·m−3)	7800
颗粒泊松比，ν	0.25
颗粒弹性模量，E / GPa	209

下载: 导出CSV 

| 显示表格

1.2   Hertz接触模型

Hertz接触模型是基于Mindlin和Deresiewicz理论建立的近似计算模型^[18−19]。采用Hertz接触模型对摩擦接触中光滑弹性球体变形进行分析，产生法向力和剪切力。当两个球形颗粒发生弹性接触时，根据Hertz接触模型的力−位移定律，颗粒间的法向力（F_n）如式（1）～式（2）所示，颗粒间的切向力（F_s）如式（3）所示。

式中：$\left\langle R \right\rangle $为颗粒等效半径，$\left\langle E \right\rangle $为等效弹性模量，δ_n为法向相对位移，R₁、R₂分别为颗粒1和2的半径；$F_{\text{s}}^0$为初始切向力，k_s为位移增量的切向刚度，Δδ_s为切向相对位移增量。其中k_s由初始切向刚度（k_s0）和历史切向加载状态决定，具体计算如式（4）所示。

式中：$ {G_1} $和$ {G_2} $分别为颗粒1和2的剪切模量，$ {v_1} $和$ {v_2} $分别为颗粒1和2的泊松比。

在颗粒生成后通过重力沉积模拟实验装粉过程（重力加速度为9.8 m·s⁻¹，方向沿Z轴向下），中止条件为不平衡力比率降低到0.00001%，此时颗粒系统达到平衡状态，沉积后的结果如图1所示。设置上模冲速度1 m·s⁻¹、方向沿Z轴向下的指令模拟压制过程，压制过程终止条件为上模冲压强达到600 MPa，此时样件基本成形且性能良好。

图  1  经过重力沉积后的离散元模型

Figure  1.  Discrete element model after gravity deposition

下载: 全尺寸图片 幻灯片

1.3   模型验证

为确保仿真数据的可靠性，需要对模型进行验证。图2所示为压制过程中的力−位移曲线，分别对仿真结果和实验结果进行多项式拟合，图2（a）力−位移曲线拟合方程为Y=9.83712−64.29832X+145.4765X²，拟合结果相关系数R²=0.99992；图2（b）力−位移曲线拟合方程为Y=30.11595−30.18838X+11.73733X ²，拟合结果相关系数R²=0.99087。仿真和实验拟合结果相关系数均接近于1，且两条曲线的趋势相近，说明模型具有较高的可信度。

图  2  力−位移曲线模型验证：（a）仿真曲线；（b）实验曲线

Figure  2.  Force-displacement curve model verification: (a) simulation curves; (b) experimental curves

下载: 全尺寸图片 幻灯片

2.   力链分析

2.1   力链的定义与筛选

在颗粒系统中，力通过颗粒之间的接触从一个颗粒传递到另一个颗粒，这一特征在被可视化之后的一个链状结构被称为力链。在金属粉末压制过程中，力链起着传递外部压制力的关键作用。目前，关于力链的定义有两种说法，一种是按照接触力进行筛选，另一种是按照应力进行筛选。Peters等^[6]首先提出通过主应力寻找力链，并发现主应力的大小既表征了力链的方向，也表征了力链的大小，并基于此提出了构成一条力链所需要的3个条件。为了更好的研究力链演化对阻塞现象的影响，本文采用主应力法对颗粒进行了分析。

主应力法的力链要满足以下三个条件^[20−21]：（1）颗粒b主应力（σ₁）大于颗粒系统的平均主应力，并将这些颗粒定为高应力颗粒；（2）颗粒a球心和颗粒b球心的连线与颗粒a的主应力（σ₁）方向的空间夹角（α）小于45°；（3）满足要求1和要求2的一条链状结构中的颗粒数目≥3。图3为一条检索出的力链，力链中的四个颗粒皆为高应力颗粒，且球心连线与主应力（σ₁）的方向夹角都小于45°。

图  3  颗粒系统中一条力链（黑色线条为接触力）

Figure  3.  Force chain in the particle system (the black line is the contact force)

下载: 全尺寸图片 幻灯片

高应力颗粒的计算公式如式（5）和式（6）所示。

式中：$\sigma _1^k$、$ \sigma _1^t $为颗粒的主应力，N为颗粒数，应力$ \sigma _{ij}^t $是X、Y和X−Y方向的应力，计算如式（7）所示。

式中：N^c为颗粒的总接触数，$V$为颗粒的体积，$ f_i^c $为作用在接触处的力的第i个分量，$ {{{\boldsymbol{r}}}}_j^c $为连接粒子中心和接触处的向量分量。

在力链识别的过程中，如果有两个颗粒满足三种条件。此时，选择条件2中角度最小的颗粒作为力链的延伸方向，如图4所示，其中α及β均为满足条件2的接触夹角，则选取较小夹角α对应颗粒为力链延伸方向。

图  4  力链颗粒延伸方向的筛选

Figure  4.  Selection of force chain particles in the extension direction

下载: 全尺寸图片 幻灯片

2.2   力链颗粒数目占比

通过力链定义中的第2条可知，力链颗粒是基于高应力颗粒筛选的，所以研究力链演化时，对高应力颗粒的分析也至关重要。由图5可知，随着压制过程的进行，高应力颗粒数目占比在初期迅速上升至40.0%以上，接着逐渐变得稳定。因为在压制刚开始的阶段，应力主要是由重力堆积产生，高应力颗粒主要集中在压坯下层，而随着压制的进行，上模冲的压力逐渐替代重力成为主要因素，使得高应力颗粒的分布逐渐变得平均、稳定。从图6可以看出，力链颗粒数目占比在上模冲压强达到90 MPa左右时达到峰值，随后开始逐渐变少，最终稳定在12.5%以上。

图  5  压制过程中高应力颗粒数目占比

Figure  5.  Number proportion of high stress particles in compacting process

下载: 全尺寸图片 幻灯片

图  6  压制过程中力链颗粒数目占比

Figure  6.  Number proportion of force chain particles in compacting process

下载: 全尺寸图片 幻灯片

2.3   稳定力链颗粒数目占比

为了研究力链的演化规律，对力链颗粒进行追踪，通过对比两个相邻数据点的力链颗粒信息，分析稳定力链颗粒数目占比，其中稳定力链颗粒是在相邻的两个数据点中都为力链的颗粒，规律如图7所示。在压制刚开始阶段，能稳定成为力链颗粒的数目占比只有20%左右，说明此时颗粒系统极不稳定，处于一种类流态。随着压制过程的进行，稳定力链颗粒在所有力链颗粒中的数目占比在初期迅速上升，随后趋于稳定，与高应力颗粒的数目占比趋势类似。最终有接近90%的力链保持不变，说明此时颗粒系统已经从类流态转变为类固态。

图  7  压制过程中稳定力链颗粒数目占比

Figure  7.  Number proportion of stabilizing force chain particles in compacting process

下载: 全尺寸图片 幻灯片

2.4   力链拓扑

为了深入研究力链演化的规律，本文引入了Cycle的概念^[22−23]，Cycle指沿一定形状的封闭循环，即在同一颗粒开始和结束的循环。通常，Cycle是对颗粒系统总体的接触力进行分析，本文在力链中进行筛选，找出力链中符合Cycle规则的个体，对它进行了量化分析。对比图8（a）与图8（b）可以发现，在粉末压制过程中，力链条数在刚开始压制时急剧增加，随着上模冲压强超过60 MPa，力链条数开始逐渐变少。当上模冲压强大于90 MPa时，力链条数依然在逐渐下降，与之相对的Cycle力链数量在上模冲压强超过90 MPa后，在保持稳定的前提下有逐步上升的趋势。Cycle力链与普通力链最大的区别就是力链判定的第2个条件，即角度。说明随着压制过程的进行，力链逐渐开始弯曲，导致有些力链不再满足条件2从而使力链条数下降。力链弯曲并出现环状是导致阻塞现象出现的重要因素之一。

图  8  压制过程中的力链条数：（a）普通力链；（b）Cycle力链

Figure  8.  Number of force chains during compacting process: (a) ordinary force chain; (b) Cycle force chain

下载: 全尺寸图片 幻灯片

3.   阻塞现象

当颗粒系统进入阻塞态时，颗粒相互交叠形成具有力学特性的刚体，最为明显的特点就是此时颗粒处于一种相对稳定的状态，对此我们通过分析不同区域下的颗粒物理现象去寻找颗粒系统在压制过程中的阻塞点。将颗粒系统轴向均匀地分为五等分，每个区域的高度随着压制过程而变化，确保不会出现测量范围的空洞。每个区域都放置五个测量球，作为测量区域，如图9所示。从上到下依次定为区域一到区域五。测量球的半径为压坯高度的十分之一，随压制过程的进行而变化。每一层测量球的球心位置在X、Y轴线上，与原点距离为10 mm，位置分布如图10所示。

图  9  测量球示意图

Figure  9.  Schematic diagram of the measuring balls

下载: 全尺寸图片 幻灯片

图  10  测量球位置示意图：（a）俯视图；（b）主视图

Figure  10.  Schematic diagram of measuring ball position: (a) top view; (b) front view

下载: 全尺寸图片 幻灯片

3.1   不同区域的颗粒速度

对颗粒系统中的五个区域颗粒速度进行分析，将颗粒分为力链颗粒和普通颗粒，求出每个区域中的力链平均速度与所有颗粒平均速度。为了研究压制过程中的力链演化与阻塞现象，重点分析轴向速度，结果如图11所示。从图中可以看出，当上模冲压强在0～30 MPa时，颗粒速度都有明显提高，这是由于突然对颗粒加压产生的冲击导致的。当上模冲压强在30～60 MPa时，颗粒速度保持稳定，在60～90 MPa时开始下降，当压强大于90 MPa时，颗粒速度会在一定范围内波动。对比图11可以发现，第一区域由于接触上模冲，所以可位移的空间比其余四个区域都大，总体速度也是最大的，与之对应的第五区域靠近下底面，所以速度也是最小的。对比每张图的力链颗粒与所有颗粒可以发现，力链颗粒的速度大于所有颗粒的平均速度，且区域越深越明显。随着区域逐渐深入，速度的波动范围逐渐变小但波动幅度越来越大，且力链颗粒比普通颗粒更加明显。综上所述，暂时将上模冲压强达到60 MPa时定为阻塞点。

图  11  压制过程中颗粒平均速度：（a）第一区域；（b）第二区域；（c）第三区域；（d）第四区域；（e）第五区域

Figure  11.  Average particle velocity during compacting process: (a) first region; (b) second region; (c) third region; (d) fourth region; (e) fifth region

下载: 全尺寸图片 幻灯片

3.2   不同区域的配位数

阻塞是颗粒系统从类流态转变为类固态的一种现象，而颗粒系统中表达结构趋于稳定最关键的一个特征就是配位数，配位数为单个颗粒的平均接触数目，配位数越大，压坯致密化程度越好，反之，压坯致密化越差^[24−25]。对不同区域中的配位数进行了分析，结果如图12所示。当上模冲压强为0～30 MPa时，第一区域颗粒的配位数变化趋势最大，因为第一区域的颗粒受重力沉积的影响最小，是五个区域中致密化最低的，初始平均配位数只有2.6。由于距离上模冲距离最近，所以配位数上升速度最快，且在压强大于60 MPa之后迅速趋于稳定。相对的，第五区域一开始受到重力沉积的影响最大，所以初始平均配位数就达到4.4，但处于最底层的第五区域受上模冲影响较小，所以整体的平均配位数上升趋势缓慢，最终上模冲达到600 MPa时，平均配位数也低于其余的四个区域。五个区域总体上来说平均配位数在0～30 MPa时都是一个急剧上升趋势，并在大于60 MPa之后逐渐趋于稳定。所以将阻塞点定为上模冲压强达到60 MPa时。

图  12  不同区域的配位数

Figure  12.  Coordination number in different regions

下载: 全尺寸图片 幻灯片

当上模冲压强达到60 MPa后，随着颗粒间的相互挤压，颗粒的配位数经过迅速增长期后变化趋势逐渐平稳，颗粒的应力变化逐渐稳定，力链开始弯曲导致阻塞现象的出现，而阻塞现象使得压坯的密度不均匀，压坯的上层比下层紧密。

4.   结论

（1）在单轴压制过程中，当上模冲压强在0～60 MPa时，力链颗粒数目与高应力颗粒数目占比都有明显的上升趋势。当压强大于60 MPa时，力链颗粒数目占比开始趋于稳定，且稳定力链颗粒数目占比增长幅度也逐渐变小。

（2）当上模冲压强大于60 MPa之后，力链条数减少，但Cycle力链数量增加，力链开始逐渐弯曲甚至出现环状结构，而环状结构是导致阻塞现象出现的主要原因之一。此时颗粒系统开始趋于稳定，颗粒间开始出现阻塞现象。

（3）将颗粒系统划分为五个区域，并对五个区域内的颗粒速度、配位数进行对比分析，发现最上层的区域速度最大，速度波动幅度最小且配位数变化最明显；最下层区域的速度最小，速度波动幅度最大而配位数变化最小。由于应力大，力链颗粒平均速度要高于普通颗粒，速度的变化幅度也比普通颗粒明显。在压制的最后阶段，五个区域中相对较为松散的区域是最下层的第五区域。