烧结是粉末冶金生产过程中最基本、最重要的环节之一，对最终产品的微观结构和性质起着决定性作用。烧结过程是多种因素(如烧结温度、烧结压力、烧结气氛、粉末粒度、纯度等)影响下的一系列物理、化学变化过程，结果是使颗粒之间发生粘结，烧结体强度增加，密度也相应提高^[1]。因此，烧结的定义^[2]为：烧结是在原子尺度上通过质量扩散的颗粒粘结的热过程，颗粒粘结使结构强度增加，系统能量降低。烧结理论研究揭示了粉末烧结过程的物理本质，Frenkel^[3]将复杂的颗粒系统简化为简单的两颗粒模型，导出了晶体颗粒粘性流动的烧结颈长速率动力学方程。此后Kuczynski^[4-5]、Kingery等^[6]和Coble^[7]运用球-板模型、球-球模型建立了烧结初期体积扩散、表面扩散、晶界扩散、蒸发凝聚等微观物质迁移机制下的颈长动力学方程，奠定了烧结理论的基础。

计算机模拟技术的出现为烧结理论的发展提供了强有力的工具。计算机模拟技术用于烧结研究的目的是对尽可能靠近实际情况的复杂烧结模型进行系统的模拟，从而对烧结现象和过程进行更深入的认识和有效的控制^[1]。

笔者通过总结近年来相场法、蒙特卡罗法、有限元法、分子动力学法、离散单元法等不同数值模拟方法在粉末烧结过程中的应用，介绍各种方法取得的代表性成果，并比较不同方法的优缺点。通过综述使读者对复杂烧结过程的数值模拟方法有较全面的认识。

1 烧结物质迁移机制

颗粒之间形成颈部后，在烧结驱动力的作用下产生颈部长大，颈长的动力是形成晶界附近的空位浓度梯度，它会驱动空位流动。Kuczynski等推导出了不同扩散机制的稳定颈长动力学方程，如体积扩散、表面扩散、蒸发凝聚、晶界扩散机制等，如图 1所示。

Kuczynski将颈长动力学方程的形式归纳为

$ {\left({\frac{x}{a}} \right)^n}=\frac{{F\left(T \right)}}{{{a^m}}}t, $

(1)

式中：指数n和m为烧结颈长指数；F(T)为仅与烧结温度有关的函数，不同的n和m值对应不同的扩散机制，如表 1所示。

烧结过程的影响因素很多，往往产生交互影响，而且烧结材料千变万化，烧结工艺也不断改进和发展。同时，烧结过程的阶段性很强，烧结的扩散机制也比较复杂，很难用某一动力学方程去描述不同材料或工艺的整个烧结全过程^[1]。因此，近年来许多研究者对多种扩散机制耦合作用下的烧结模型进行了研究，程远方等^[8-9]对21世纪的烧结理论进行了总结。张道礼等^[10]考虑表面扩散和体积扩散建立了烧结初期综合作用机制下的动力学方程。Wakai等^[11]对晶界扩散和表面扩散耦合机理进行了研究。Maximenko等^[12]得到烧结第一、二阶段考虑表面、晶界和体积扩散的动力学方程。吴艳青等^[13]引入考虑表面能级数形式的应力函数，建立热压条件下表面扩散与应力场耦合的力学模型，分析表面扩散和弹性变形共同驱动的致密化规律。Djohari等^[14-15]在细观尺度对粘性流动、表面空位扩散和晶界扩散机理进行了研究。Zhang等^[16-17]研究了二维和三维情况下考虑表面扩散和晶界扩散的两颗粒模型，给出了颗粒形状变化和表面曲率变化的规律。Klinger等^[18]研究考虑表面和相间扩散的两非混相颗粒烧结机理。

2 烧结的计算机模拟方法

烧结理论的研究使烧结过程的物理本质逐渐清晰，为使用计算机模拟来预测烧结颈长规律提供了可靠的物理模型。计算机模拟技术在计算机科学、数学、力学等学科的发展和进步，为进行烧结过程模拟仿真研究提供了有利的基础和工具^[19]。Nichols等^[20]以表面扩散物质迁移为主导机制，首次用数字计算机计算了球排、球板烧结颈长的颈部轮廓变化趋势，表明颈部轮廓是连续变化的。在模拟烧结变化过程中涌现了大量的数值方法，如蒙特卡罗(Monte Carlo，MC)模拟、相场法(phase field method，PFM)、有限单元法(finite element method，FEM)、分子动力学(molecular dynamics，MD)和离散单元法(discrete element method，DEM)等。

2.1 蒙特卡罗法

蒙特卡罗方法^[19]也称为随机模拟方法，是一种以概率统计理论作为指导的数值计算方法，通过多次重复抽样来求解确定性问题，当样本足够大时，具有很高的精度。由于其算法的简洁性，对于复杂的烧结模型反而具有更高的计算效率。Qiu等^[21]使用多状态Potts动态模型，以温度、能量最近邻域加权以及界面、表面能之比为参数，对复杂颗粒系统的晶界迁移和蒸发凝聚机制进行了模拟，并研究了无量纲参数$\frac{{{k_{\rm{B}}}T}}{J}$对纳米TiO₂颗粒颈部尺寸的定量变化规律的影响，当$\frac{{{k_{\rm{B}}}T}}{J}$=0时，颈部几乎不发生长大，因为此时没有驱动力来打破颈部平滑的凹表面，而随着$\frac{{{k_{\rm{B}}}T}}{J}$的变大，这种平衡被打破，颈长规律近似为一组斜率为0.22的随时间变化的平行线。笔者取$\frac{{{k_{\rm{B}}}T}}{J}$=0.5分别模拟了烧结过程中两颗粒、三颗粒和多颗粒系统烧结过程中颈部长大过程，如图 2所示。牛玉^[22]通过引入多种烧结机制和动态守恒算法，完善了经典烧结模型，解决了蒙特卡罗算法的质量不守恒问题。刘建元^[23]研究了单个圆形颗粒在液相中的溶解与析出过程及颗粒熟化长大，并对固相烧结中后期晶粒长大、孔洞迁移以及空位湮灭等烧结微观结构演变行为进行了模拟。Luque等^[24]考虑颗粒的几何特征，定义熔化和凝固概率分布函数研究液相烧结时固、液之间的状态改变。Aldazabal等^[25]针对二面角模型的不完整性，提出了考虑溶质的扩散和表面能作为参数的蒙特卡罗模型模拟液相烧结过程。Moria等^[26]采用粘塑性有限元方法计算宏观非均匀收缩，而微观的晶界生长和气孔消失则采用蒙特卡罗方法模拟。Braginsky等^[27]基于动力学蒙特卡洛方法在二维情况下模拟晶粒生长、空位扩散和晶界孔隙湮灭的微观结构演化，获得了烧结压力，并且给出了细观尺度和宏观连续介质模型之间的关系。马非等^[28]基于经典的晶粒生长机理，采用改进的蒙特卡罗方法模拟了不同的烧结温度和保温时间下的晶粒生长演化过程，并用Al₂O₃陶瓷烧结实验进行验证。

2.2 相场法

相场法^[29]是以经典热、动力学为理论基础的一种半唯象方法，通常采用吉布斯-朗道型模型，通过设定相关的场变量，把空间内的微观结构作为一个整体，用偏微分方程来描述微观结构演化过程，从而能够详细了解任意时刻的微结构形貌和空间分布情况。在烧结过程模拟中可以进一步获得颗粒尺寸、粒子分布及随时间变化规律等微观信息，该方法可以实现单相/双相材料的烧结、晶粒长大以及多种材料的烧结过程模拟。Kumar等^[30]研究了非等径颗粒烧结过程中的晶粒生长过程，将烧结和晶粒生长过程分为3个子过程：颈长、缓慢的晶粒粗化过程、晶粒的快速迁移，如图 3所示。分析不同子过程中的热力学驱动力，得到了晶粒生长速度的理论解析式见式(2)，晶粒生长速度强烈地依赖于颗粒的几何尺寸，其中v为晶粒生长速度，m为晶界迁移率，γ_s和γ_gb分别为表面能和晶界能，X，r_s和r_l分别为颈部半径、小颗粒半径和大颗粒半径，β/2为小颗粒与颈部的圆心角。谢辉等^[31]引入了Cahn-Hillard方程控制的相对密度场和Ginzburg-Laudau方程控制的长程取向场参量，模拟了表面扩散和晶界扩散机制共同作用下的四球模型和五球模型中Mo粉末压坯微粒间烧结颈长大以及孔隙的演化过程，并指出表面扩散是主导机制。贾磊等^[32]建立了表面扩散、晶界扩散、体积扩散和蒸发-凝聚机制共同作用下的Mo粉末烧结模型，观察到烧结颈的形成与长大，气孔的封闭、球化和消失以及晶粒的吞并与长大等现象，并指出前期烧结颈长的主导物质传输机制为表面扩散，气孔封闭受体积扩散和界面迁移率的影响。黄创高等^[33]采用一组随时间和空间连续变化的取向场变量和浓度场变量来表征烧结过程微结构特征，用扩散界面相场模型研究陶瓷烧结过程中五球模型内部的气孔扩散与演化过程。刘明治等^[34]采用Cahn-Hilliard方程和Allen-Cahn方程来控制相对密度场和长程取向场的变化，建立了原子扩散系数和晶界扩散、表面扩散和体积扩散的函数关系式，模拟两相多孔材料烧结过程中的微观组织演化。Zhang等^[35]考虑表面扩散和晶界扩散，认为表面区域是相同化学势不同组分的混合相，研究Fe和Cu烧结过程中的微观结构和溶质分布。Asp等^[36]考虑原子的空位扩散机制，研究固态烧结时润湿的影响，认为固体的表面将出现连续变化的特征。Wang^[37]考虑粉末颗粒的刚体平移和旋转、边界迁移、晶粒生长和各种扩散机制建立固相烧结模型，发现刚体运动在粉末压块致密化过程中扮演着重要的角色。

$ v=m\left[{{\gamma _s}\left({\frac{1}{{{\gamma _s}}}-\frac{1}{{{\gamma _l}}}} \right)+{\gamma _{gb}}\frac{{2\left| {\cot \left({\beta/2} \right)} \right|}}{X}} \right]. $

(2)

2.3 有限元法

有限元方法是将所研究的区域进行离散化处理，根据本构关系和相应的初始边值条件，采用多项式插值算法来求解各个单元中状态变量的近似解。该方法已广泛应用于烧结过程的温度场、应力场和微观结构演变方面的模拟。Fukui等^[38]基于有限元方法数值模拟了多物理场作用下氧化锆(YSZ)和钛(Ti)混合颗粒在放电等离子烧结过程中颗粒变形和氧扩散，结果如图 4所示。在外加压力为15 MPa、温升率为160 ℃/min、烧结温度为1 000 ℃，对氧化锆体积分数为12.5%和25%的多颗粒系统进行了研究，塑性变形区集中在钛颗粒接触部，但塑性变形不能使空洞完全消失。当温度达到800 ℃时，颈部空洞会快速收缩，同时，氧会从氧化锆颗粒扩散在钛颗粒，当温度达到1 000 ℃，钛颗粒中氧均匀分布。Djohari等^[14-15]在细观尺度下研究了空位扩散的粘性流动机制、晶界扩散机制和致密化过程。Jeong等^[39]基于塑性理论，将静水压力添加到一般的多孔材料的屈服函数，提出了一个统一的模型来描述压实和烧结，采用有限元分析来预测不均匀密度分布和最终产品的尺寸变化。Kiani等^[40]提出用致密化数据(密度作为时间的函数)代替本构模型来预测烧结变形，但仅限于无压烧结。Leclerc等^[41]基于体积、表面、晶界和气孔扩散的本构模型模拟固相低压力烧结致密化过程。Kraft等^[42]研究了SiC粉末的综合烧结模型，对固态烧结的形状进行了预测。郑涛等^[43]以Fe粉和Al₂O₃粉为原材料建立了磁性磨粒热压烧结过程的传热模型，对模具温度场进行分析，从而优化热压工艺。宋毅等^[44]将粉末材料的热电性能参数假设为温度和密度的双变量函数，提出了改进的放电等离子烧结粉末材料的有限元模型，实现了温度场、电场和位移场的耦合。Bruchon等^[45]采用基于欧拉描述的有限元方法计算表面和晶格扩散。Pan等^[46-49]采用耦合的有限元和有限差分法计算二维情况下的表面和晶界扩散，并得到了气孔消除的动力学方程。Park^[50]采用晶粒生长本构模型来研究W-Ni-Fe合金的液相烧结过程，考虑重力、基底摩擦、表面张力和固体含量的影响。

2.4 分子动力学方法

分子动力学方法是一种确定性模拟方法，通过求解体系内所有粒子的运动方程，从而获得与粒子运动路径相关的基本过程，最终得到相应的组织演变行为^[51-52]。对烧结过程的模拟，其空间尺度大多在微米乃至毫米级。吴茂等^[53]采用嵌入原子势能模型，研究直径为3~10 nm的纳米金颗粒的熔化及烧结过程，在烧结初期，两颗粒迅速形成烧结颈，之后存在烧结颈相对稳定阶段；当达到起始烧结温度后，烧结颈逐渐长大；达到金纳米颗粒表面层原子开始熔化温度，烧结颈快速长大；达到熔点时，烧结颈剧烈增加，此时颗粒完全融化，模拟过程如图 5所示。Ding等^[54]采用分子动力学方法探讨了烧结过程的分子动力学模型和连续扩散模型之间可能的区别。Moitra等^[55]考虑温度、压力和粒子结构的影响，对晶体钨粉末在亚微观层面的颈长和收缩过程进行模拟。Cheng等^[56]认为氩族纳米丝和纳米颗粒混合系统比纳米颗粒和纳米颗粒混合系统更易烧结。Bai等^[57]研究了温度梯度驱动晶界迁移，并且验证二氧化铀颗粒的温度梯度驱动力方程。Tsuruta等^[58]发展了SiC晶粒的n阶并行算法来模拟颈长过程。闫秋实等^[59]对碳化硅纳米颗粒的电火花烧结的初期与中期阶段进行了模拟分析，研究了温度、压力等因素对晶粒增长、晶界形成的影响。Xu等^[60]采用多颗粒模型研究了纳米镍颗粒和镍/氧化锆颗粒的烧结过程。Song等^[61]建立了三级烧结方案，充分考虑分层结构研究烧结过程中扩散系数的影响。

2.5 离散元方法

离散元方法以牛顿第二定律为基础，把离散体看作是由大量独立颗粒单元组成的集合体，每个单元具有单独的物理属性，通过不断循环计算模型中每个颗粒单元的受力和位移变化来对整个系统的运动变化过程进行模拟^[62]。Wonisch等^[63]基于离散元方法研究了氧化铝各向异性烧结过程，并确定了本构参数。Martin等^[64]介绍了一个简单的粗化模型并用它模拟了各向异性收缩的致密化过程。Fu等^[65]在不同的烧结温度下对3个空间方向上粘结剂去除后的收缩行为进行了分析，讨论了各向异性收缩和颗粒微结构之间的关系。Nosewicz^[66]提出一个采用大增量步的粘弹性离散元模型，并模拟了NiAl粉末的烧结过程。Henrich等^[67]研究了细观尺度下的颗粒重排对宏观致密化速率、体积和剪切粘度的影响。Olmos等^[68]对金属和陶瓷(Cu-Al₂O₃)混合粉末烧结现象进行研究并与实验结果进行对比。

2.6 其他方法

对复杂烧结过程的模拟还有其他的方法，如元胞自动机法^[69]，但该算法在时间和空间上受到局部性限制，难以很精确地模拟全局的动力学行为^[22]。此外，Wakai等^[70]用Surface Evolver Program来研究非等径颗粒表面扩散和晶界扩散过程，并讨论了晶界能和表面能比值的影响。Zhu等^[71]采用自编程来模拟考虑表面扩散和晶界扩散的非等径颗粒颈长过程。

2.7 多种方法优缺点比较

烧结是一个多颗粒系统在多物理场作用下发生一系列物理化学变化的复杂过程。不同的数值模拟方法一般都能模拟典型的微观结构演化，如烧结颈部的形成和长大、晶界迁移和晶粒生长以及孔隙消除现象，不同的方法都还存在着自身的优势和劣势^{[19, 22-23]}。蒙特卡罗方法程序结构十分简单，模拟过程中可以直接观察晶粒形貌，计算精度高，可以模拟颗粒系统和复杂几何结构的烧结过程；但是获得模拟参数比较困难，计算时间长，尤其对三维情况计算效率偏低。相场法对界面做弥散化处理，各个参数的物理意义清晰，界面两侧成分和结构连续变化，不需要跟踪界面，计算效率较高；但是用于描述微结构演化的偏微分方程数学推导过程非常复杂, 严格求解复杂几何结构和多种烧结机制的共同作用有很高的难度，可以通过GPU并行计算来提高计算效率。有限元方法需要求解大量方程，计算量非常大，对新的烧结理论本构关系的加入较为困难，需要编写复杂的子程序进行再次研发，这使有限元方法在烧结领域的应用有明显的局限性。分子动力学方法能够分析微观粒子的物质扩散路径，得到颗粒的微观组织演变过程，在模拟烧结过程时计算空间尺度较大，需要极大的计算量，一般只用于分析少数颗粒模型，对复杂几何结构和大规模颗粒系统烧结难以计算。离散元方法是从细观尺度研究颗粒材料的有效方法，在计算烧结过程时由于颗粒数量巨大，求解时间步较短，需要很大的计算量，不适合单独用于烧结过程的分析，但可以与其他方法联合使用，从而更准确地模拟烧结过程。Shinagawa^[72]提出了一个精确的相场法和离散单元法耦合方法，对等径和非等径两球模型以及多球模型晶界迁移和对心收缩进行模拟，采用PFM计算烧结颈表面和内部的烧结力，再将此结果代入到DEM模型中，很好地描述了烧结颈长和收缩行为，并且验证了烧结团聚和晶粒异常长大现象，如图 6所示。

3 结语

烧结理论的研究进展为使用计算机模拟烧结颈长的动态规律提供了可靠的物理模型，计算机技术的快速发展为分析复杂的烧结过程提供了强大的模拟工具。不同的数值模拟方法如相场法、蒙特卡罗法、有限元法、分子动力学法、离散单元法均可以用来模拟复杂的烧结过程；但受算法本身和颗粒烧结过程复杂性的限制，很难完全实现烧结过程模拟，可以通过联合不同方法来模拟烧结过程。此外，需要与实验相结合来验证模拟过程的准确性和可靠性。