随着工业技术水平的蓬勃发展，人们对粉末冶金制品性能的要求越来越高。在粉末压制成形过程中，由于粉末颗粒具有明显的离散性，其压制行为可以从三个尺度来区分，即颗粒内部分子等微观运动、颗粒为特征的细观运动和颗粒体系的宏观运动。不同尺度间的运动具有一定的关联性，粉末颗粒宏观性质更大程度上取决于颗粒相互作用形成的细观结构^[1]。对粉末颗粒压制过程和机理的研究多采用实验的方法，耗费了大量的人力与物力^[2]。近几十年来，研究者们发现有些数值模拟方法可以从颗粒层面将粉末压制过程更好的展现出来，有助于更清楚地了解细观层面颗粒的运动规律。目前，粉末压制成形工艺主流的数值模拟方法有三种，分别为有限单元法、离散单元法和多粒子有限单元法。

有限单元法（finite element method，FEM）是一种连续体力学研究方法，是解决高度非线性问题的强大数值方法。这种方法是将粉末视为各向同性介质，直接通过本构方程描述粉末致密行为，不考虑颗粒之间相互接触的作用力^[3]。该方法很难从更小的尺度去探讨成形致密化过程中粉末颗粒的形状、粒径以及颗粒的分布对成形质量的影响。离散单元法（discrete element method，DEM）是建立在牛顿第二定律上，用来研究非连续性颗粒物质结构和运动规律的一种数值方法。这种方法根据颗粒间重叠量计算接触力，以此更新每个颗粒的速度和位置^[4]。该方法对于颗粒的大变形问题无法准确描述^[5]，但在颗粒材料力链相关研究中应用较广，因为粉末材料的特殊性由其“颗粒性”决定，离散元法正是研究颗粒物质运动规律的一种数值方法。多粒子有限单元法（multi-particle finite element method，MPFEM）是将物体中的每个粒子都视为可变形的连续体，对每个粒子进行有限元网格划分，从而分析粒子之间的接触和变形行为。在细观领域，考虑粉末材料颗粒的离散特点，对于每个颗粒进行建模，通过颗粒间作用力模型，进而推导大量颗粒的压制过程。该方法能够从微观层面研究颗粒的大变形行为，但是需要极高的计算成本，因此，只能对模型的一个子域的粉末结构进行建模，进而推导整个压制过程的细观力学行为。

由于有限单元法将颗粒团视为一个整体，不适合从细观层面分析粉末颗粒参数对粉末压制最终成形质量的影响。虽然离散单元法对于高致密、大变形等问题的模拟精度较低，但是该方法却是研究颗粒物质力链形成、演化的重要手段。多粒子有限单元法可以解决颗粒大变形带来的仿真结果不够精确的问题。因此，本文将从离散单元法在颗粒力链研究和多粒子有限单元法在颗粒大变形研究两个方面来论述数值模拟在粉末压制细观分析中的应用。

1.   离散单元法在颗粒材料压制细观研究中的应用

颗粒材料是离散的、固体的、宏观的颗粒集合体。由于它们的离散性，每个粒子只能与有限数量的邻居接触。在受外部载荷作用时，力从一个边界到另一个边界的传递只能通过这些粒子间的接触发生，力的传递路径称之为力链。力链对外部载荷的变化比较敏感，不断进行断裂重组，以更合适的链条匹配外力变化。力链是颗粒材料稳定性、弹性和流动性等力学性能的决定因素^[6]。

粉末压制所使用的材料属于颗粒材料，外力的变化会导致粉末内部接触力的变化，从而促使力链演化，而力链网络为压制力的传播提供了路径。因此，研究力链问题对粉末压制中力的大小、力的传递以及致密化程度等至关重要。颗粒之间的摩擦、施加压力和加载方式等因素均会影响到力链演化过程。

1.1   摩擦对力链的影响

任何两个表面不光滑的物体接触运动都会产生摩擦，粉末颗粒也不例外。在粉末压制中，颗粒在压力的作用下，相互接触、挤压形成，摩擦会对力链的形成与发展造成影响。Guo^[7]采用离散元法研究了密实颗粒材料中临界力链长度及其屈曲问题，发现颗粒间摩擦对临界力链长度没有直接影响，但会影响颗粒柱屈曲前的曲率。基于同样的方法，孟凡净等^[8]研究了碳钢粉在不同摩擦系数条件下对力链方向性的影响，并将接触力和力链强度联系起来。结果表明，随着摩擦系数的增大，颗粒的力链方向倾向于压应力的方向，且颗粒间接触力越大，力链强度越强。张炜等^[9]研究了铁粉间摩擦系数和铁粉与侧壁间摩擦系数对力链量化特征的影响。结果表明，较低的粉末间摩擦和粉末与侧壁摩擦，力链的数目和方向、整体承载不均匀度和屈曲度均比高摩擦时要好，更有益于致密化过程。虽然摩擦系数小对粉末压制更有利，但是在粉末压制过程中随着压力的增大，摩擦系数不是一成不变的。Zhang等^[10]研究发现，当力链充满颗粒系统的时候，在两种摩擦模型下，平均颗粒间摩擦系数会达到一个稳定值，如图1所示。除了颗粒间摩擦系数之外，颗粒与模具壁之间的摩擦系数和力链也存在关系。Zhang等^[11−12]研究了不同压实速度下侧壁摩擦系数的变化和力链定量特性演化规律，得到了侧壁摩擦系数先增大后趋于稳定的结论。在压实速度变化较小时，压实速度对力链特性和侧壁摩擦特性影响较小。

图  1  颗粒摩擦系数与轴向应变关系^[10]

Figure  1.  Relationship between the particle friction coefficient and axial strain^[10]

下载: 全尺寸图片 幻灯片

综上，摩擦力与颗粒间接触力及摩擦系数有关。随着压制的进行，颗粒之间的接触增多，接触力和摩擦系数也会增大，直到压制过程结束，压坯形成，接触力和摩擦系数趋于稳定。

1.2   压力对力链的影响

压力是粉末压制不可或缺的重要因素。Peters等^[13]采用离散元法研究了压力对力链长度和力链数量的影响，发现力链长度成指数分布的规律，力链数量随着压力增大呈减小的趋势。Bassett等^[14]使用Community Detection技术和Geographical Null模型模拟了无摩擦条件下力链结构与压力的关系，发现力链长度在不同压力下均成指数分布，但力链的几何形状受压力影响较小。Huang和Daniels^[15]在这项技术的基础上，研究了摩擦系数（μ）和压力（P）对力链的影响，得到了摩擦系数和压力的临界值，临界值两侧的力链大小和强度发生了转变。

王飞等^[16]采用离散单元法研究发现，随着压制的进行，力链从一开始的松散状态变成力链强弱区别明显的状态，强力链起支撑作用，弱力链起辅助作用，这时力链分布比较稳定。张炜等^[17]模拟了在双轴压缩条件下，四周压力对铁粉力链的影响。结果表明，随着压力的深入，强力链沿压制方向逐渐增多，然后屈曲、断裂重组，最终达到强力链中包含弱力链的稳定状态，如图2所示。张超等^[18]使用铝粉研究了不同冲击速度对侧壁压力分布的影响，发现随着压制进行，大量弱力链在强力链周围生成，进一步证明了强弱力链互相补充支撑，共同维持粉末整体的稳定。

图  2  轴向应变与力链演变^[17]：（a）轴向应变0%；（b）轴向应变4%；（c）轴向应变7%；（d）轴向应变13%

Figure  2.  Evolution of force chain with the different axial strain^[17]: (a) axial strain 0; (b) axial strain 4%; (b) axial strain 7%; (b) axial strain 13%

下载: 全尺寸图片 幻灯片

1.3   加载条件对力链的影响

在细观分析中，外界的一个小的条件变化可能给力链特性带来不一样的变化。在实际压制过程中，有普通压制和高速压制之分，有单向压制和双向压制之别，不同加载方式对力链的影响也不尽相同。Zhang等^[19]采用离散单元法进行研究，发现准静态双轴条件下，力链数的演化和偏应力的演化是一致的，力链屈曲会带来粒子动能的增加，屈曲发生在孔隙率较大的地方；在动态冲击条件下，力链屈曲会导致局部膨胀，引起能量耗散，影响压制效果。张炜等^[20]采用同种方法研究发现，在高速压制过程中应力传递与力链演化保持较为一致的趋势。该课题组还研究了高速压制过程冲击速度（v）对力链特性的影响^[21]。结果表明，初始速度越高，力链的数目、强度、准直性都会变大，如图3所示。普通压制时，粉末初始平衡状态改变比较小且缓慢，而高速压制使初始平衡态迅速被打破，力链瞬间进行断裂重组，最终重新生成新的力链网络。

图  3  力链特性的变化^[21]：（a）数目；（b）长度；（c）强度；（d）准直性

Figure  3.  Changes of the force chain properties^[21]: (a) number; (b) length; (c) strength; (d) collimation

下载: 全尺寸图片 幻灯片

1.4   粒径比对力链的影响

在实际生产中，粉末尺寸并不一致，粉末粒径也会对力链的形成与发展产生影响。王海陆等^[22]采用离散单元法研究了铁粉粒径分布对力链分布的影响，发现单一粒径粉末强力链数目较多，多粒径粉末弱力链较多，在压实完成后给强力链的支撑起到了辅助作用。两者进行比较，多粒径更容易压实。

除了利用数值模拟研究力链外，还可以通过实验进行研究。Zhang等^[23]实验研究了纯剪切系统和各向同性压缩系统中的力链特性。Gendelman等^[24]利用外力和粒子接触方向解决了粒子之间法向力和横向（摩擦）力的确定问题。Iikawa等^[25]利用受重力垂直敲击的光弹性圆盘系统研究了有序到无序过程中力链取向参数的变化。

对力链的研究有助于揭示力的传播途径，如果能够使力更容易传播到压坯底部，粉末压制会容易许多，也有助于阐明颗粒材料宏观响应的细观机制，为后续对力链的研究打下基础。目前还有很多因素并没有考虑到或研究的还不够深入，比如材料属性不同，力链的演变过程或者量化描述是怎样变化的；颗粒的形状不同，会对力链特性造成哪些影响；不同粒径比粉末力链长度、数目和强度与单一粒径有什么不同等。综上所述，在粉末压制领域，利用离散单元法研究与探索力链已经取得了一定成果，但是想要更清楚的从细观层面揭示颗粒离散性对粉末压制样品质量的影响还不够，尤其是涉及到颗粒大变形时，此时就需要用到多粒子有限单元法。

2.   多粒子有限单元法在粉末压制中的应用

2.1   多粒子有限单元法数值模拟的难点

尽管多粒子有限单元法已用于粉末压制的研究中，但是粉末压制过程是一个非常复杂的高度非线性、非连续性的数值模拟过程，目前还存在许多难题。

首先是颗粒物质随机分布的建模。建立物理模型是数值模拟的第一步也是最关键的步骤。想要从细观角度研究粉末相对密度的影响因素就不能简单把颗粒群作为连续体处理，需要建立颗粒模型。为了使粉末颗粒之间的排列更符合实际情况，需要随机生成，而目前商业化有限元软件较难实现颗粒随机性排列建模。当前主流办法是先使用离散单元法生成随机颗粒，再导入有限元分析软件中模拟。其次是颗粒在压制过程中的大变形问题。由于多粒子有限单元法结合了离散单元法和两者的优点，既可以研究颗粒的位移，也可以研究颗粒的变形问题。而且商业有限元软件都带有自适应网格技术，结合多粒子有限单元法可以很好地解决粒子大变形带来的计算不精确问题。再次，在粉末压制成形过程中，由于颗粒间的摩擦系数受润滑条件、接触压力、颗粒流动行为等因素的影响，且颗粒-壁面摩擦关系复杂，颗粒之间的接触状态一直处于变化当中，所以颗粒间摩擦系数的大小难以确定，其影响因素主要包括施加载荷的速度、压制力、模具的质量和粗糙度等。李达等^[26]根据平面剪切式摩擦测定方法，利用自行设计的模壁摩擦系数测试装置，测定了不同性质粉末以及不用压制力下的摩擦系数变化规律，为后续模拟分析摩擦系数提供参考。

除此之外，粉末颗粒形状的不确定性导致模拟结果与实际压制情况有一定出入，主要表现在颗粒形状的不规则性、颗粒与颗粒之间形成的拱桥以及压实后颗粒的大变形等，并且由于粉末颗粒数量过于庞大，如果每个颗粒都单独划分网格，对计算机的要求非常高。在目前大多数研究中，将颗粒半径放大，建立数量比较少的颗粒模型，将粉末颗粒理想化成圆形或者球形颗粒是一种常态。Peng等^[27]为了方便运算，将颗粒简化成半径均为1 mm的球形，建立了190个球体模型进行仿真模拟。研究表明，由于弹性应变能的释放，卸载后粉末压坯存在回弹性的问题，将会导致压缩坯体尺寸变化，密度和内部应力状态重新分布，并且回弹量与最终尺寸、最终密度之间的定量关系难以明确。Zhou等^[28]考虑了压坯在卸载过程中的回弹问题，发现回弹率与弹性应变能折减率相关，与卸载过程中静摩擦力密切相关，整体回弹率随摩擦系数和单位压实能量的增大而减小。

除以上这些难点之外，从细观角度研究粉末压制过程的理论基础还不完善，还存在颗粒间相互作用固定化、颗粒粘结、颗粒之间形成拱桥等很多现象难以准确描述。

2.2   多粒子有限单元法在粉末压制中的应用

虽然多粒子有限单元法还存在许多难题，但许多学者采用该方法对单一粉末、混合粉末等各种不同粉末进行了细观尺度的仿真模拟，讨论了初始堆积结构、粉末含量、粉末粒径比和烧结温度等众多因素对粉末压实过程中的影响，揭示了粉末细观运动的机理，对提高材料的最终相对密度有很大的指导意义。

2.2.1   颗粒粒径和颗粒含量对压实过程的影响

Han等^[29−30]、Huang等^[31]、Wang等^[32]和Feng等^[33]采用多粒子有限单元法分别研究了Fe/Al、Al/SiC、TiC/316L和NaCl/Al混合粉末粒径比和硬质颗粒含量对压坯相对密度分布的影响，发现粒径比增加，最终相对密度有所下降；硬质颗粒含量增加，致密化所需的压力也就越大，相同压力下相对密度也会有所降低。图4、图5和图6分别为Fe/Al粉末、TiC/316L粉末和NaCl/Al混合粉末相对密度与压力关系曲线，其中R_M/R_N为M粉末与N粉末的粒径比。

图  4  压实压力和Fe/Al粉末粒径比（R_Fe/R_Al）对压坯相对密度的影响^[30]

Figure  4.  Influence of compaction pressure and Fe/Al powder size ratio (R_Fe/R_Al) on the relative density of compact^[30]

下载: 全尺寸图片 幻灯片

图  5  压实压力和TiC/316L粉末粒径比（R_316L/R_TiC）对压坯相对密度的影响^[32]

Figure  5.  Influence of compaction pressure and TiC/316L powder size ratio (R_316L/R_TiC) on the relative density of compact^[32]

下载: 全尺寸图片 幻灯片

图  6  NaCl/Al混合粉末的压实压力和相对密度实验数据和仿真结果^[33]

Figure  6.  Experimental data and simulation results of compaction pressure and relative density for the NaCl/Al mixture powders^[33]

下载: 全尺寸图片 幻灯片

Xu等^[34]分别从2D和3D角度模拟研究了SiC/6061Al混合粉末在热压过程中的致密化行为。结果表明，压力不变，提高温度可增加相对密度；另一方面，随着粒径比的增加，相对密度略微降低，而且随着SiC含量增加，抗压强度上升，导致相对密度难以提升，进一步说明了粉末粒径比和硬质颗粒含量对相对密度的影响。由此可知，为了使压坯质量得到保证，满足实际生产需求，粉末粒径比和硬质颗粒含量存在一个最佳范围。当粒径比超过适当值，小颗粒会被大颗粒“保护”起来，进而造成小颗粒周围的孔隙无法填充，相对密度不高；当硬质颗粒超过适当值，压制完成时的应力释放会使粉末压坯面临碎裂的风险。

在实际生产中还有很多是单一粉末压制的情况。Li等^[35]采用多粒子有限单元法模拟了三种粒径分布（单一粒径分布、二元粒径分布和正态粒径分布）的Ti6Al4V粉体对热等静压过程的影响。研究发现，粒径比越大，小颗粒与大颗粒变形程度的差距越大，大颗粒与小颗粒圆度的差异越大。对于单粒径粉体的填充结构，其致密化机理是颗粒重排破坏桥接结构和大孔隙，从而加速致密化过程；对于二元粉末或正态分布粉末的Ti6Al4V致密体，其机理可以推断为小颗粒的极端变形促进大颗粒的相对刚性运动。

2.2.2   温度对压实过程的影响

Jia等^[36]对不同粒径钨粉的单轴模压和固态烧结过程进行了多粒子有限单元法数值模拟，得到了一组可以获得压制样品质量较好的压制和烧结参数。Zou等^[37]采用同种方法系统地研究了压力和温度对热等静压过程中钨粉致密化行为的影响。结果表明，温度越高，达到相同相对密度所需的压力越小，并且压制结束后颗粒的总应变能越低，如图7和图8所示。由此可知，适当的温度可以使粉末更容易产生变形，进而填充孔隙，提高相对密度。

图  7  不同温度下相对密度与压力关系^[37]

Figure  7.  Relationship between relative density and pressure at the different temperatures^[37]

下载: 全尺寸图片 幻灯片

图  8  总应变能随压力和温度变化关系^[37]

Figure  8.  Relationship of the total strain energy, pressure, and temperature^[37]

下载: 全尺寸图片 幻灯片

2.2.3   堆积结构对压实过程的影响

粉末颗粒的堆积结构在自然状态下是不一样的，给随后的压缩过程带来不一样的变化。Zhang等^[38]采用多粒子有限单元法研究了准静态条件下等粒径铜颗粒在单轴模压成形过程中的致密化过程，图9展示了其随机结构的压缩过程。在有序堆积结构压实中，致密化主要是颗粒的塑性变形，而对于随机和蜂窝堆积结构，致密化是由压缩变形和颗粒之间的相对滑动和滚动共同作用的结果。Zou等^[39]研究了冷等静压过程中各向同性和各向异性的初始堆积结构和烧结温度对钨坯及烧结件性能的影响。结果表明，随机堆积结构压实后内部平均应力比有序堆积结构要大，而且不同结构会导致烧结达到相同相对密度所需的温度不同。不同初始充填结构压实后，压实力链和压实力传递也不同。在力链中局部应力的驱动下，形成“拱”结构的颗粒被迫移动以填满下方的大孔隙，随着压力的进一步增大，颗粒发生塑性变形填满小孔隙。

图  9  随机结构压制过程^[38]

Figure  9.  Random packing compressing process^[38]

下载: 全尺寸图片 幻灯片

2.2.4   其他因素对压实过程的影响

Loidolt等^[40−41]采用多粒子有限单元法研究了粘结接触和应变路径对粉末压实屈服性能的影响。结果发现屈服面和粘结强度成函数关系，如图10所示，其中C_max为粘结强度，P为压力，q为等效应力。屈服面随相对密度的变化规律基本一致，屈服面取向和接触内聚力无关。以上研究考虑到了粉末颗粒之间粘滞性的问题，但作者只研究了铜材料，研究结果对于其他材料不一定适用。

图  10  屈服面和粘结强度函数图^[40]

Figure  10.  Yield surfaces as the function of cohesion strength^[40]

下载: 全尺寸图片 幻灯片

基于同样的方法，Zhou等^[42]研究了卸压过程中等粒径铜球粉末宏观和细观性质的演化，结果如图11和12所示。卸压后，压坯的内应力和位于压坯中心的一些颗粒的内应力完全释放，但相邻颗粒接触区域出现较大的残余应力，最大应力出现在颗粒边缘，而且孔隙填充量越小，压力释放后应力降低越大。Zhou等^[28]研究了卸压后回弹率（δ）的变化及单位质量能量（E_m）和生坯密度的联系，结果如图13和图14所示。如图所示，存在一个最佳单位质量能量值对应相对密度最高值。卸压之后的弹性后效是一个颗粒内部应力释放的过程，由于颗粒的弹塑性，这一过程无法避免，只能尽量减小回弹量，以提高最终相对密度。

图  11  孔隙边缘节点Von Mises应力分布^[42]

Figure  11.  Von Mises stress distribution at the pore edge nodes^[42]

下载: 全尺寸图片 幻灯片

图  12  应力释放前后节点Von Mises应力^[42]

Figure  12.  Nodal Von Mises stress before and after stress release^[42]

下载: 全尺寸图片 幻灯片

图  13  单位质量能量和摩擦系数对回弹率的影响^[28]

Figure  13.  Influence of unit mass energy and friction coefficient on springback^[28]

下载: 全尺寸图片 幻灯片

图  14  卸压后回弹初始模型^[28]

Figure  14.  Initial model of resilience after unload^[28]

下载: 全尺寸图片 幻灯片

数值模拟中有很多现成的本构模型及摩擦模型可以选用。Güner等^[43−45]研究了不同材料模型和摩擦模型对铜球粉末压实过程的影响。结果发现Wanheim-Bay模型和Amonton-Coulomb模型适用于变形较大的颗粒，Von Mises材料模型最适合多粒子有限单元法模拟仿真。除此之外，Peng等^[27]从初始堆积密度角度对190个等粒径W−Cu颗粒进行模压模拟分析和实验验证，结果如图15所示。引入了粉末的平均旋转度和平均等效应变来定量表征颗粒重排和变形，将致密化过程分成三个阶段：颗粒重排、铜颗粒变形、钨颗粒变形，如图16所示。发现提高初始堆积密度和粉体均匀度可提高产品的质量。

图  15  W−Cu颗粒模压模拟和实验验证^[27]

Figure  15.  Molding simulation and experimental verification of W−Cu particles^[27]

下载: 全尺寸图片 幻灯片

图  16  致密化过程中粒子的平均等效应变^[27]

Figure  16.  Average equivalent strain of particles during densification^[27]

下载: 全尺寸图片 幻灯片

总体来说，颗粒的压制都会经历颗粒重排和颗粒变形两个阶段，研究人员利用多粒子有限单元法探讨粉末颗粒细观层面的压制过程，将众多可能影响到粉末压坯质量的因素考虑在内，揭示了颗粒粒径与含量、温度和初始堆积结构等对粉末相对密度的影响。随着研究的越发深入，数值模拟越来越接近实际，逐渐将粉末压制中存在的因素考虑到模拟中，得到的结果对于指导实际生产、提高最终产品性能有重要意义。

3.   结论与展望

（1）利用离散单元法从细观层面对力链特性及力链演化过程进行探讨，揭示了力链数目、长度和准直性等性质对样品最终质量的影响，对颗粒细观运动有更进一步了解。

（2）利用多粒子有限单元法模拟粉末压实并进行了实验对比分析，通过对模型的简化有效解决了粉末颗粒的随机性、颗粒间摩擦、颗粒形状粒径以及颗粒卸载后回弹等问题，获得了有利于指导实际的结论。但对颗粒间相互作用固定化、颗粒粘结、颗粒之间形成的拱桥等难点的研究还不够深入。

（3）无论是对力链的研究还是对颗粒压制过程的研究，外部施加的载荷都是恒定，在动载荷条件下（粉末电磁脉冲压制成型）和恒定动态载荷条件下粉末压制过程中力链演化与颗粒变形致密的关系是未来研究的方向。