金属增材制造数值模拟研究进展

王 岩 ，刘雨萌 ，刘江伟 ，魏瑛康 ，张亮亮 ，王建勇 ，尚巍巍 ，刘世锋 ✉

1) 西安建筑科技大学冶金工程学院,西安 710055

2) 山东大学能源与动力工程学院,济南 250061

3) 陕钢集团产业创新研究院有限责任公司,汉中 723000

20 世纪80 年代出现的快速成形技术（rapid prototyping,RP）发展为今天的增材制造技术（additive manufacturing,AM），在此几十年中，新材料、新工艺、新设备的迅速更新变革使得增材制造技术成为机械加工和锻造等传统金属材料加工的重要补充，形成了减材、等材和增材的鼎足之势。根据施加热源和材料形式的不同，可以分为选择性激光烧结技术（selective laser sintering，SLS）、选区激光熔化技术（selective laser melting，SLM）、电子束选区熔化技术（electron beam selective melting，EBSM）、激光立体成形技术（laser solid forming，LSF）、电子束熔丝沉积技术（electron beam fusion，EBF）、电弧增材制造（wire and arc additive manufacturing，WAAM）以及直接金属粉末激光烧结（direct metal laser sintering，DMLS）等。

相较于传统成形工艺，增材制造技术具有设计自由、材料利用率高、个性化定制和快速完成复杂构件等优势，被应用于航空航天[1-4]、汽车工业[5-6]、医疗[7-8]和核电[9]等领域。航空航天领域，欧洲空中客车集团创新中心（Airbus Group Innovations）采用粉末床熔合（powder bed fusion，PBF）增材制造技术成功制备出经过拓扑优化的空客320 和380 飞机的钢制舱门托架和发动机舱门铰链，使其减重60%，在承受同样外部载荷的情况下，最大应力减小49%[3]。西北工业大学黄卫东教授和林鑫教授团队采用激光立体成形技术（laser solid forming，LSF）制造出的C919 飞机Ti6Al4V 合金翼肋缘条，经过中国商飞的性能测试，完全满足设计要求[4]。汽车工业领域，Juechter 等[6]采用粉末床熔合技术生产出性能优异的钛铝合金涡轮增压器轮毂。医疗领域，增材制造技术主要应用于与骨科、牙科等相关的临床领域[10]，通过对孔结构的设计，促进与天然骨力学性能相匹配的骨替代植入物尽可能快地被宿主组织吸收、降解和替换[11]。核电领域，中国广核集团应用选区激光熔化技术成功生产出核电站复杂流道仪表阀阀体，经过基础性能测试，该阀体的材料化学成分和基础力学性能满足国际核电标准RCC-M 的要求[12]。

在增材制造过程中，金属在电子束/激光辐照下历经一系列复杂的物理过程，包括能量吸收、传热/传质、熔化/凝固等。除材料本身特性和成形工艺特点外，掌握和调控增材制造过程中低层物理机制，可实现对增材制造构件织构[13]、组织和性能[14]的调控和优化，同时可实现对裂纹、孔隙和球化等缺陷的优化和消除[15]。数值模拟是理解金属增材制造过程中发生的复杂物理过程并为工艺条件优化提供指导的有力工具[16]。数值模拟分析可以针对实验技术存在的稳定性不足、可重复性差、分辨率受限、可观测区域限制及设备成本昂贵等问题，帮助理解和分析增材制造过程中物理状态的变化，指导优化工艺过程。此外，非接触式空气耦合的超声无损检测技术[17]、电荷耦合器件（charge-coupled device，CCD）高速相机技术[18]和高速高能X 射线原位成像技术[19]的发展，为认识和理解增材制造过程中复杂的物理过程提供了技术支撑，保障了增材制造构件的成形精度和质量。

1 宏观尺度

宏观尺度中用的数值模拟方法主要为有限元法。在几乎所有的工程学科中，有限元法已被证明是解决边界、初始和特征值问题的基本工具。在有限元法中，连续的物质（区域）被离散成有限个单元，形成一个网格，从而将连续的问题简化。因此，在涉及多种物理过程的增材制造过程中，有限元法能避免进行大量物理实验，成为预测增材制造工艺产生的残余应力和变形的最常用技术。但是其不足之处是无法满足复杂的边界条件，这是制约其发展的重要因素。宏观尺度中，主要从温度场和残余应力两个方面展开论述。

1.1 温度场

温度及温度梯度是金属增材制造过程中重要的物理量，与材料密度、表面张力、热导率、热容量或热扩散率等密切相关，产生不同的热力学、动力学及机械效应，从而决定了构件的最终质量。因为热源峰值温度非常高（高达金属沸点），温度梯度和加热/冷却速率极大（分别有106K·m-1和106K·s-1），所以实验中温度场的精确测量比较困难[20]，通常采用数值模拟的方法分析增材制造过程中的温度场。因此，在增材制造数值模拟过程中，光束热源模型至关重要。增材制造中热源模型主要分为二维面热源模型[21]和三维体热源模型[22]。表面热源模型认为热源作用于粉末颗粒的表面，热量通过材料热传导向粉末内部和基板传递。二维表面热源通常假设为高斯分布[23]，对于基本横向电磁模式（TEM00），在任意点（x，y）的光束辐照度可以用式（1）表示。

式中：I为光束辐照度，（x，y）为热源相对于扫描速度的坐标，t为时间，β为激光束的吸收率，P为固定激光源的功率，rl为激光束的半径，v为扫描速度。

考虑到激光扫描可以穿透金属粉末床[24]，在二维面热源模型的基础之上建立了三维体热源模型。例如对于常见的激光粉床熔融（laser powder bed fusion，LPBF）增材制造过程，激光能量不仅沉积在粉末床的表面，而且与粉末床内部颗粒相互作用。因此，必须考虑体积热源来描述激光穿透粉末的情况[25]。Zhang 等[26]采用图1 所示的圆柱形、半球形、半椭球形、锥形、辐射传递法、射线追踪法、线性衰减法和辐射衰减法八种热源模型对激光粉床熔融进行数值模拟。通过实验对仿真结果进行了验证，与实验相比，所有八个热源模型都导致熔池深度浅40%以上。为了提高模型的性能，研究者提出了各向异性增强导热系数和吸收率的模型，并将其应用于指数衰减热源的传热模拟，使得熔池宽度和深度的平均误差分别为2.9%和7.3%，提高了热源模型的准确性。

图1 热源模型示意图[26]：（a）圆柱形；（b）半球形；（c）半椭球形；（d）锥形；（e）辐射传递法；（f）射线追踪法；（g）线性衰减法；（h）辐射衰减法Fig.1 Schematic of the heat source models[26]:(a)cylindrical shape;(b)semi-spherical shape;(c)semi-ellipsoidal shape;(d)conical shape;(e)radiation transfer method;(f)ray-tracing method;(g)linearly decaying method;(h)exponentially decaying method

增材制造过程中的工艺参数和工艺过程对温度梯度的分布有显著影响。在研究中发现，在一定功率下，熔池底部的最大温度变化率随扫描速度的增加而增加，但在达到某一临界速度后又减小。在低扫描速度下，不同激光功率下的最大温度变化率之间的差异很小，随着扫描速度的增加，差异变大，激光功率和最大温度变化率之间的关系变得复杂[27]，且最高温度对激光功率的敏感度比扫描速度高[28]。分形扫描方式得到的温度场比S 型扫描方式更加对称[29]。扫描道次对温度场也有一定的影响。考虑到熔化潜热和与温度相关的热物理特性，Tang 等[30]采用高斯体热源分别对衬底-粉末模型、全粉末模型和部分固体模型下的熔池温度场和形貌进行了单轨道和多轨道扫描分析。图2 为273 K 下部分固体模型的多道温度场轮廓图，研究发现单道扫描时，温度峰值和最高冷却速率随预热温度的升高而升高，而最高加热速率则呈现相反的趋势；在多道扫描策略中，每道的温度最大值随着模型中固体体积的减少而增加，且随着扫描轨道数量的增加而增加。此外，预热温度可导致最高温度提高[31]。

图2 273 K 下部分固体模型多道温度场轮廓图[30]：（a）t=5.040 ms；（b）t=15.022 ms；（c）t=25.026 ms；（d）t=35.132 ms；（e）t=45.182 msFig.2 Contour of the multi-channel temperature field at 273 K in the partial solid model[30]:(a)t=5.040 ms;(b)t=15.022 ms;(c)t=25.026 ms;(d)t=35.132 ms;(e)t=45.182 ms

1.2 残余应力

残余应力是当材料与周围环境达到平衡时，留在材料内部的应力[32]。由工艺引起的残余应力所导致的零件变形和分层等缺陷仍然是阻碍增材制造发展的重要挑战之一。目前，主要通过改变工艺参数、调控热输入分布来降低残余应力[33]。

影响残余应力的因素有扫描矢量的长度[34-35]，扫描策略[32]，扫描速度以及粉末床的预热温度。Liu 等[36]发现扫描矢量长度从42 mm 减少到18 mm时，残余应力减少了50%以上。Cheng 等[33]发现45°扫描方式在纵向和法线方向上的残余应力值最低。Dunbar 等[37]将（全长）光栅逐层旋转67°，变形幅度降低了37.5%，增加扫描速度可以有效减小残余应力。Vastola 等[38]对Ti6Al4V 的单道电子束熔化进行了残余应力形成的系统计算研究。研究者先将模型与已发布的实验数据匹配，研究了标称工艺参数的残余应力，显示扫描区域周围存在明确的热影响区；随后，从标称参数开始进行计算，研究光束尺寸、光束功率密度、扫描速度和床预热温度对改变应力场的影响。Farahmand 和Kovacevic[39]研究了扫描速度、扫描功率和粉末床预热温度对残余应力的影响，发现粉末床预热温度影响最大，粉末床温度越高，残余应力越小；在较低的扫描速度下，横向和纵向残余应力更高。Li 等[40]用有限元模型（finite element method，FEM）来预测选择性激光熔化过程中的应力和变形，分别采用了图3 所示的四种不同的扫描策略并且分析对应模式下的残余应力（图4）。结果表明，水平连续扫描模式比垂直和岛状扫描模式更能减小基板变形。因此，可以通过缩短扫描矢量，提高扫描速度，调整扫描方式以及提高粉末床的预热温度减小残余应力。

图3 四种扫描策略示意图[40]：（a）水平扫描；（b）垂直扫描；（c）连续扫描；（d）“最小”热影响Fig.3 Four scanning strategies[40]:(a)horizontal sequential;(b)vertical sequential;(c)successive;(d)“least”heat influence

图4 四种扫描策略的残余应力示意图[40]：（a）水平扫描；（b）垂直扫描；（c）连续扫描；（d）“最小”热影响Fig.4 Residual stress contours for the four scanning strategies[40]:(a)horizontal sequential;(b)vertical sequential;(c)successive;(d)“least”heat influence

在粉末床熔合期间积累的热应力会导致零件变形，目前的计算技术需要大量的时间来预测粉末床熔合增材制造构建过程中产生的零件畸变，新模型的提出可以在保证正确性的前提下有效缩短计算时间。Peng 等[41]用一个热回路网络（thermal circuit network，TCN）模型来估计粉末床熔合期间的部件温度历史，与传统的有限元模型热建模相比，采用热回路网络模型预测金属粉末床熔合增材制造零件的温度历史，在牺牲了不到15%的精度下，计算速度快了两个数量级以上。Li 和Anand[42]提出了一种新颖的后向插值（backward interpolation，BI）模型来快速估计零件变形，预测的变形结果与实际构件的变形结果吻合较好，使用台式计算机在几分钟内获得了零件失真值，证实了使用这种方法显着减少了计算时间。

2 介观尺度

宏观尺度主要针对的是构件尺寸级别，该尺度模拟有助于理解和优化增材制造中的温度场和残余应力，扮演着举足轻重的角色。但对于粉末颗粒，宏观尺度的数值模拟无法满足对熔池特征以及熔池动力学的研究[43]，介观尺度的模拟可以有效解决这个问题，本结主要从粉末颗粒和熔池行为进行综述。在介观尺度中，在考虑熔池流体动力学的情况下，使用的数值方法包括格子玻尔兹曼（lattice Boltzmann，LB）和有限体积法（finite volume model，FVM）。格子玻尔兹曼法的基本思路是划分网格，与有限体积法等传统方法相比，具有以下优点：能处理具有复杂边界条件的过程，编程、前后处理简单、易实现，不需要对计算网格进行转换，可直接模拟有复杂几何边界的连通域流场。有限体积法的优点是可以在质量、动量、能量守恒方程的基础之上将浮力驱动流动、表面张力、Marangoni 对流等复杂的物理现象包括在内。因此，如果研究对象有复杂的集合边界可以采用格子玻尔兹曼法，若研究的重点在各种物理现象对增材制造的影响可以选用有限体积法。当情况复杂时，可以用多种方法进行耦合，满足所研究的内容。

2.1 粉末颗粒

铺粉会显著影响后续操作过程以及最终产品的质量。现有的主要铺粉方式有辊式和刮刀式两种。Nan 等[44]考虑了气体雾化金属粉末通过辊子铺粉的动力学，并通过自由落体离散元（distinct element method,DEM）分析了间隙高度和辊子转速对演化粒子轨迹和铺展层均匀性的影响。结果表明，由颗粒对流/循环引起的狭窄间隙中的瞬时堵塞和铺展堆中的尺寸分离对铺展层的均匀性产生不利影响。偏析程度随着间隙高度的增加或辊子转速的降低而减小。Yao 等[45]发现在316L 不锈钢粉末铺展中，刮刀移动速度越小，粉床质量越高，但工作效率越低。当刮刀移动速度为0.1 m·s-1，增加刮刀间隙高度或减小粒度（即D=30 μm）将增加平均相对堆积密度和结构均匀性。Chen 等[46]对刮刀式铺粉过程进行了实验和计算模型相结合的研究，确定了三种主导粉末扩散过程的沉积机制：内聚效应、壁效应和渗透效应。导致颗粒团聚并因此降低堆积密度的内聚效应随着颗粒尺寸的减小而变得更强；导致粉末层中空位更多的壁效应随着层厚度的减小或粒度的增加而变得更强；双峰粉末颗粒存在渗流效应，导致粉末层内颗粒偏析，从而降低堆积密度。三种沉积机理在散粉过程中相互竞争，对粉层的堆积密度产生综合影响。Zhang 等[47]对辊式和刮刀式的尼龙粉末铺展情况进行模拟，发现在较厚的粉末层中，由于压实作用，通过辊式铺粉可以获得较高密度的粉末层。随着粉末层厚度的增加，辊式铺粉比刮刀式铺粉对偏析指数更敏感。对粉床的数值模拟研究可以发现，不同铺粉方式的最终目的是获得更均匀、更致密的粉床，更高的铺粉效率。

但是，在增材制造数值模拟过程中，往往会忽略颗粒的大小分布以及相对密度变化的随机分布，采用大小均匀的有规律粉床填充简化粉床模型进行数值模拟[48]，但这与实际粉床的分布状态严重不符。因此，一些学者进一步考虑粉末床中颗粒的随机分布建立模型，Körner 等[49]评估了随机粉末床的重要性，并应用降雨模型来生成二维模拟的粉末床，该模型还可扩展到三维[50]。降雨模型的建立过程如图5 所示，通过计算第一接触点的垂直位置，将每个粒子分别放置在前一层上，然后向下旋转，直到达到稳定状态。这种方法的自然相对密度在二维或三维中约为74%或60%。为了调整到粉末床颗粒的合理密度（一般为45%～60%），去除部分粉末颗粒，直到达到所需密度[49]。对颗粒相互作用力的研究[51]揭示了粘结力对粉末床最终相对密度的关键作用，并将其纳入自由落体离散元方法中以得到更真实的粉末分布。对复杂分布过程的简化方法是采用自由落体离散元法[52]代替滚动和倾斜过程。每个粒子都能在受边界墙限制的连续空间中运动，除重力外，法向力和切向力相互作用于质点，改变质点运动。将有限体积法和自由落体离散元法结合可以模拟出粉末床熔合过程中粉末颗粒的沉积过程[53]以及粉末粒度对粉末流动行为的影响[54]。

图5 随机粉床带旋转的随机包装降雨模型（a），降雨模型产生的粉床（b），去除部分颗粒并调整相对密度（c），真实粉末床（钛合金）截面（d）[49]Fig.5 Schematic of the rain model for random packing with rotations(a),powder bed produced by the rain model(b),adjusting the relative density by removing some of the particles(c),cross section of a real powder bed (titanium alloy)(d)[49]

2.2 熔池行为

在增材制造加工过程中，粉末颗粒吸收激光束或电子束的能量，熔化即可形成熔池。熔池中的对流依赖于黏度，并由不同的外力（如重力、浮力、表面张力、毛细作用、Marangoni 效应或蒸发压力）驱动[48]，如图6 所示。与其他力相比，重力的作用较小[55]，可以忽略。材料的热膨胀引起浮力并产生热应力。金属的高表面张力与润湿能力相结合，使得熔池表面光滑而稳定。相反，不稳定的熔池易分裂，表面张力会导致单个熔球的形成，即球化效应[56-57]。Marangoni 力可以让流体远离熔池中心的温度峰值，并增加热量传输。在增材制造加工过程中，因为熔池具有较高的温度，材料蒸发产生的反冲压力能驱动流体运动。这种现象在选区激光熔化过程中尤为明显，反冲压力使熔池的形状为锁孔状而非碗状[58]。在选区激光熔化过程中，激光束穿透到材料的某一层厚度时会发生蒸汽毛细现象[59]，具有挥发性的元素发生选择性蒸发，这会改变局部或是整个构件的物质组成[60]。

图6 部分熔融粉末床中熔化过程的主要物理现象[48]Fig.6 Dominant physical phenomena during melting illustrated in the partially molten powder bed[48]

熔池几何形状受激光强度参数[61]、表面张力[62]、Marangoni[63]、激光功率[64]和扫描速度[65]的影响。随着激光强度参数值的增加，熔池宽度增加，熔池深度变浅。然而，对于较大的强度参数值，激光功率强度的径向分布导致在远离对称轴一定距离处形成具有较大深度的熔池[61]。表面张力使熔池的纵向形态沿凝固轨迹纵向变化，且影响熔池深度的周期性变化。Marangoni 和反冲压力使液态金属表面呈现凹形，还能增大熔池的宽度。激光功率可以增大熔池宽度和深度[64]。Shi 等[66]研究了圆形高斯（circular Gaussian，CG）和椭圆（transverse elliptical，TE和longitudinal elliptical，LE）两种形状的激光束对熔池几何形状的影响，图7 为激光功率350 W，扫描速度0.75 m·s−1下圆形高斯、横向椭圆和纵向椭圆型的熔池深度（沿构建方向），分别为54、36、80 μm，其中横向椭圆光束的熔池最宽（沿横向方向），而纵向椭圆光束的熔池最深（沿构建方向）。Zhang 等[67]为了研究熔池热流场及其对Inconel718合金选区激光熔化过程的影响，建立了一种与传热和流体流动耦合的物理模型。模拟结果表明，熔池表面温度梯度引起的表面张力驱动Marangoni 对流，使选区激光熔化过程中的流体流动状态主要是向外的对流。Marangoni 对流包括对流和传导热通量，两者都对熔池形状有影响，但对流热通量的影响占主导地位，因为其幅度比传导热通量大一个数量级。对流热流加快了金属液的流动速度，有利于散热。对流热流使熔池更宽，而传导热流使熔池更深更宽。此外，多次扫描引起的热积累增加了对流和传导热通量，导致熔池宽度和深度增加，但对流热通量对熔池形状的主导作用没有改变。

图7 熔池几何形状与激光束形状的关系[66]Fig.7 Melt pool geometry as the function of the laser beam shape[66]

熔池动力学主要由毛细管力、Marangoni、蒸发压力以及对粉末颗粒和前一层的润湿能力驱动。Zhao 等[68]发现熔融金属的Marangoni 效应主要决定流体行为，是影响CoCrMo 合金选择性电子束熔化（selective electron beam melting，SEBM）熔池不稳定性的关键因素。凝固速率对凝固前沿的局部流体流动很敏感。凝固前沿的主动流体对流增加了新晶粒形成的概率，从而抑制了柱状晶粒的外延生长。图8 上部显示熔体轨迹的3D 视图，下部显示2D 中心纵向横截面，带有速度值的箭头是流体速度矢量的尺度，电子束流功率（P）为600W，速度（V）为300 mm·s-1。Scharowsky 等[55]用高速摄影机对电子束熔化过程进行观察，并分析熔池寿命、大小和振荡。Körner 等[49]发展了一个2D 格子玻尔兹曼方法模型来研究在选择性光束熔化过程中，随机堆积粉末床在高斯光束照射下的熔化和再凝固。这种方法使得许多在标准连续图像中无法描述的物理现象变得容易理解，例如对粉末密度的影响、随机填充粉末床的随机效应、毛细现象和润湿现象。比较不同束流功率和扫描速度下Ti6Al4V 的熔池深度和宽度，表明数值模拟结果与实验结果非常吻合，应将熔池温度限制在低于4000 K 的实际值内，由于空间分辨率从5 μm 降低到1 μm，蒸发反冲压力计算量增加[48]。

图8 具有多种物理场的模拟熔池[68]：（a）Marangoni 效应、电子反压力和蒸汽反冲压力；（b）没有电子反压力；（c）没有蒸汽反冲压力；（d）没有Marangoni 对流Fig.8 Simulated molten pool with the multiple physics[68]:(a)Marangoni effect,electron counter pressure,and vapor recoil pressure;(b)without electron counter pressure;(c)without vapor recoil pressure;(d)without Marangoni convection

3 微观尺度

近十年来，微结构模拟仿真的兴起是增材制造技术和模拟仿真方法快速发展的结果。与铸造相比，增材制造通过在更小的区域内控制工艺参数，为冷却条件的可控性提供了更多的可能。模拟仿真在小区域微观模拟仿真方面的重要性越来越明显，有助于实现对增材制造加工组织进行定量控制的目标[69]。目前在增材制造凝固组织模拟方法中，使用最多的是元胞自动机（cellular automata，CA）法。此外，相场（phase field，PF）法在物理模型表示方面具有较高的精度和分辨率，可以捕获凝固组织中的亚晶粒特征，但相场模型通常需要较高的计算成本。蒙特卡罗（Monte Carlo，MC）法在模拟晶粒生长行为，特别是再结晶行为方面简单高效，但只能提供大量的“试错”后的结果。图9 显示了每个模型的基本思想[70]。元胞自动机模型的基本思想在于通过许多小的几何对象（2D 中的正方形，3D 中的八面体）的叠加来再现生长的树枝状晶粒的包络，每个对象都由元胞自动机的一个单元单独控制。相场的基本原理是一个平滑变化的函数，称为“相场”，描述了两个流动相之间的界面。蒙特卡洛波茨模型不考虑凝固过程中的枝晶、偏析和过冷，图9（c）显示了在凝固过程中模拟的晶粒结构演变。熔池内的所有细胞都有不同的状态（“自旋”）。在凝固过程中，由于表面能效应，结构开始粗化，粗化仅在热影响区内进行。在较低温度下，即在热影响区之外，粗化停止并且晶粒结构稳定[70]。

图9 不同建模方法示意图[70]：（a）元胞自动机；（b）相场法；（c）蒙特卡洛波茨模型Fig.9 Schematic of the different modeling approaches[70]:(a)cellular automaton;(b)phase-field method;(c)Monte Carlo Potts model

3.1 元胞自动机法

元胞自动机法是一种基于变换规则反映单元集合状态的算法，典型的元胞自动机模型由元胞、元胞状态、元胞空间、元胞邻居、规则函数和时间函数组成。在过去的二十多年里，元胞自动机已经成功地用于微观结构模拟，如动态再结晶[71]和晶粒生长行为[72]。随着凝固模拟技术的快速发展，元胞自动机法被广泛应用于凝固组织模拟，研究不同参数对所得微结构的影响[55]。Zhang 和Zhang[73]将计算机流体动力学（computer fluid dynamics，CFD）和元胞自动机法结合预测粉末床熔合316L 不锈钢凝固组织演变，如图10 所示。在凝固过程中，界面处液体到界面处固体的转变采用“修正偏心方形”方法进行建模，较高的扫描速度会形成更细的晶粒尺寸、更大的激光晶粒角和更高的柱状晶粒含量。Ao 等[74]使用元胞自动机建立了预测合金微观结构与传热和熔池运动耦合的模型，并对合金的微观形态和形成机理进行了研究，发现凝固过程中冷却速率约为105～106K·s-1，等轴晶随着预热温度的升高或扫描速度的降低而增加，随着扫描间隔的增加，外延柱状晶粒变得越来越长且越来越窄。

图10 计算机流体动力学模型（a）、元胞自动机模型（b）和316L 元胞自动机模拟晶粒结构的俯视图（c）[73]Fig.10 CFD model (a),CA model (b),and the top view of the CA simulated grain structure of 316 L stainless steel (c)[73]

3.2 相场法

相场法的理论基础是Landau 提出的朗道连续相变理论[75]，是解决界面问题的一种强有力的计算方法[69]，它提供了一种通过模拟来更好地理解相变机理的途径。相场法模拟主要适用于增材制造凝固过程中的晶粒长大和枝晶形貌。

增材制造过程中晶粒的相场模拟主要有两种方式。方式一，在相场模型中将温度设为与温度梯度和速度有关的项。柱状晶沿热流方向生长，温度梯度越高，其生长速度越快[76]，还与热梯度和冷却速率[77]有关。Sahoo 和Chou[78]通过相场模型（phase field model，PFM）研究温度梯度和扫描速度对微结构的影响。对于实际合金，大多数相场模型只取两种主要的合金元素进行计算，对多合金元素的处理需要探究。方式二，先用相场模型方法进行温度场模拟，将得到的相关数据导入相场模型中用来研究凝固组织生长情况。Karayagiz 等[79]在相场中提取温度梯度（G）和生长速率（R），并将两者导入相场模型中，发现在高的线能量下，冷却速率越高，柱状晶越来越细。Liu 等[80]将相场中得到实际的温度分布导入相场模型中，这比将温度场的温度梯度导入相场模型更接近实际的增材制造过程。增材制造过程中的枝晶形貌相场模拟与枝晶生长有关，Acharya 等[81]研究了增材制造过程中沿扫描方向和垂直扫描方向的枝晶生长情况。Yang 等[82]采用相场法研究了镍基高温合金凝固过程中的三维组织形成，为了加速大规模相场模拟，还提出了一种利用图形处理器（graphics processing unit,GPU）实现并行计算的方法。在等温凝固模拟中，Yang 等[82]详细研究了不同过冷程度下单枝晶和多枝晶组织的演变规律，记录并分析了枝晶尖端的生长速度和固相分数；在定向凝固模拟中，他们还研究了不同温度梯度下初生枝晶臂的形成，在温度梯度为50 K·mm-1时，DD6 单晶高温合金定向凝固一次的枝晶形貌如图11 所示，由图可知模拟得到的显微组织与实验观察吻合较好。

图11 （a）温度梯度为50 K·mm-1 时模拟枝晶形貌俯视图；（b）DD6 单晶高温合金定向凝固一次臂阵列的显微形貌[82]Fig.11 Top view of the simulated dendrite morphology at the temperature gradient of 50 K/mm (a)and the microstructure of the directional solidified primary arm array of DD6 single crystal superalloy[82]

3.3 蒙特卡洛法

蒙特卡洛模拟涉及使用随机数和概率来解决基于试错的复杂问题。在晶粒生长模拟领域，蒙特卡洛法的离散计算域为相同且规则的模式。模拟域中的每个网格都分配有代表不同晶体的整数。那些具有相同晶体取向的相邻网格属于一个晶粒。基于系统自由能的最小化，随机选择的节点将改变其方向，并可视化为晶粒生长行为。为了实现增材制造过程的凝固微观结构模拟，Rodgers 等[83]采用随机蒙特卡罗方法（stochastic Monte Carlo）和波茨蒙特卡罗方法（Potts Monte Carlo），建立模拟凝固和固态演化的三维模型，该模型结合了基于有限差分的热传导求解器，以创建完全集成的微观结构预测工具，揭示了之前未知的重熔循环次数和柱状晶粒尺寸之间的相关性。Sunny 等[84]根据特定的热历史，采用动态动力学蒙特卡罗（kinetic Monte Carlo，KMC）框架预测晶粒形貌的显著变化，为研究人员提供了评估增材制造零件不同区域力学性能变化的新能力，如图12 所示，采用直接能量沉积（direct energy deposition，DED）技术制备的304L 不锈钢的实验结果和动力学蒙特卡罗模拟结果吻合。

图12 直接能量沉积制备的304L 试样实验与KMC 模拟电子背向散射衍射图像比较[84]Fig.12 Comparison of the electron back-scattered diffraction images between the 304L specimen made by DED and the KMC simulation[84]

现有的凝固组织模拟方法各有优点，但应用在增材制造中也存在一些缺陷。元胞自动机法比相场法计算成本低，相比动力学蒙特卡方法具有更坚实的物理基础，是目前最适合增材制造凝固组织模拟的方法。相场方法是建立在物理模型的基础上，再现率高，它还能够描述高阶枝晶、微偏析、固相转变等亚晶粒特征。但计算量大只能在超级计算机上实现，限制了它的应用推广。元胞自动机方法研究部分和整体之间的关系，认为系统的物理状态是一个聚合的本地状态改变；相场方法试图通过连续的场参数来描述系统的状态，从而避免了界面尖锐所带来的问题；对于蒙特卡洛方法，它使用大量的采样来逐步逼近实际。虽然这三种方法之间存在着显著的差异，但它们在增材制造微观结构模拟领域都发展迅速。

在微观尺度的数值模拟，充分应用元胞自动机法、相场法和蒙特卡洛法可以详细了解微观组织的影响因素，预测晶粒的形貌变化和工艺−微观组织−性能关系，有助于提高对缺陷的认知和对工艺的优化。

4 多尺度

由于材料在不同尺度上存在结构差异，如传热/传质，应力应变，微观结构和宏观结构，导致了对应模拟方法和模型的差异。单个尺度的研究较为普遍并且取得了一定的进展，可以用来预测大多数的实验现象。但在研究温度场、熔化/凝固等物理现象时，宏观尺度受微观尺度影响，单一尺度的模拟必然会造成其它尺度信息的缺失，因此多尺度方法模拟可以实现在精确、有效和真实描述之间的平衡[85]。金属增材制造跨尺度数值模拟是力求从多个角度，尽可能全面地对增材制造过程进行研究，宏观尺度可模拟温度和应力，介观尺度可揭示粉末床上所发生的热演变过程，微观尺度可观察机理。

4.1 宏观-介观尺度

电子束熔化（electron beam melting，EBM）为功能梯度材料（functionally gradient materials，FGM）的制造开辟了新途径，但是温度场变化较快使得测量精度不够，Yan 等[20]开发了一个新的多尺度传热建模框架以研究制备功能梯度材料，如图13 所示。热源模型源自微尺度电子−材料相互作用模拟，能够结合材料特性、电子束特性和入射角，从力学角度描述加热现象，计算材料性能和电子束性能。将热源模型应用于介观尺度下单个粉末颗粒的热演化模型中，阐明了不同材料、不同尺寸的混合粉末颗粒的熔融和凝聚过程。另外，研究者还进行了另一种介观尺度模拟，以评估原始粉体床的有效热导率，建立了一个宏观传热模型，该模型可跟踪粉床的聚结状态，以确定粉体床的有效材料性能，并将对熔池尺寸的预测与已发表的实验结果进行比较和验证。

图13 多尺度模型框架示意图[20]Fig.13 Schematic diagram of the proposed multi-scale modeling framework[20]

4.2 微观-宏观尺度

增材制造生产的各种合金具有许多潜在的应用，对于产品质量和安全性要求较高，需要精确的数值模拟来设计结构和部件，应全面考虑合金的特定微观结构以获得足够的模拟精度。然而，适用于增材制造的合金晶体塑性模型和模拟尚未完全建立。Zhang 和Andrä[86]基于三维代表性体积元（representative volume element，RVE），采用快速傅里叶变换（fast-Fourier transform，FFT）方法对激光粉床熔融制造的AlSi10Mg 合金的晶体塑性模型进行求解，利用原位同步X 射线衍射实验结果对模型进行了校准，系统定量地分析了Al 相的硬化参数、Si 颗粒的有效长径比和孔隙率对宏观和微观力学性能的影响，为解释增材制造AlSi10Mg 合金的微观结构-性能关系提供了一种有效的工具，由此产生的宏观材料特性可用于优化宏观组件和结构。图14（b）和图14（c）为随着施加应变大小的变化，真实应力和晶格应变的曲线图，可以看出预测的曲线和文献中的测量数据非常吻合。电弧送丝增材制造技术（wire and arc additive manufacturing，WAAM）可以制备出具有独特力学行为的结构，Watanabe 等[88]从多尺度角度研究电弧送丝增材制造多层钢结构的力学行为。多层结构的拉伸试验可用来评估非均质性对材料宏观性能的影响，通过与拉伸试验的相场模拟对比，研究了非均质多层钢结构的变形行为。结果表明，多层钢结构的非均质性比其设计更为复杂，其宏观力学性能超出了微观力学估计的上下限。

图14 激光粉床熔融制造AlSi10Mg 合金的代表性体积元图（a）、真实应力图（b）和晶格应变图（c）[86-87]Fig.14 RVE (a),true strain (b),and lattice strains (c)of the AlSi10Mg alloys prepared by LPBF[86-87]

4.3 微-介-宏观尺度

Gu 等[89]为包括铝基合金（AlSi10Mg）、镍基高温合金（Inconel718）、陶瓷颗粒增强铝基和镍基复合材料在内的金属材料选区激光熔化提供了多尺度建模和协调控制。在微观尺度上研究了氮化铝颗粒在选区激光熔化Al 基纳米复合材料中的迁移和分布机理，以及WC/Inconel718 复合材料中增强层与基体之间梯度界面的原位形成。在介观尺度上揭示了AlSi10Mg 和Inconel718 合金粉末的激光吸收和熔化/致密化行为。最后，在宏观尺度上，提出了激光局部扫描过程中的应力发展规律和选区激光熔化加工复合材料变形参数依赖控制方法。多尺度数值模拟和实验验证有助于监测激光加工过程中复杂的粉末−激光相互作用、热传质行为、微观组织和力学性能的发展。

5 结论与展望

本文综述了宏观、介观、微观及多尺度下的增材制造过程数值模拟，包括宏观尺度下的温度场和残余应力，介观尺度下的粉末颗粒和熔池行为，微观尺度下的凝固现象以及跨尺度下所涉及的多物理现象。与实验试错相比，金属增材制造技术数值模拟节约了时间和成本，得到了飞速发展，有助于优化工艺参数并提高最终打印部件的机械和几何质量，但仍有几个方面需要解决：首先，计算成本高，运算时间长。怎样在保证精度的前提下降低计算成本是需要迫切解决的问题，并行技术和云计算技术的应用和发展是解决这一问题的思路。其次，金属增材制造过程中涉及多种互相影响的复杂物理现象，要提高模拟结果的精确度需要用到复杂的物理模型，对边界条件、材料特性和模型所做的简化可能在一定程度上引入了的差异，无法通过实验验证。数据挖掘能够在因大量数据导致模拟过程耗时耗力的问题方面发挥重要作用。第三，模拟结果对不同案例研究的可移植性是有限的。在该领域的不同研究人员之间很难找到材料特性建模的共同点，根据不同的材料特性建立物理模型库，便于数据的提取和选用。第四，金属增材制造模拟的专用软件局限于产品设计解决方案，而工艺模拟领域正处于开发的最初阶段，关于工艺模拟专业软件的开发应该提上日程。第五，宏观尺度的仿真模拟会被更广泛的应用于生产实践中，为金属增材制造过程中产生的变形、表面粗糙度大、孔隙和裂纹等问题提供切实可行的解决方案；介观和微观分析的重心将逐步从小范围的科研转向工业化的生产中。金属增材制造过程的模拟会更加系统化，材料的选择、粉末颗粒的分布、构件和工艺的设计、工艺参数的设置、性能的分析以及后处理等因素也将引入仿真模拟中，以确保高效地制备出具有优异性能的构件。总而言之，金属增材制造多尺度及跨尺度的数值模拟是实验和模拟的有机结合，不同尺度之间的耦合会是金属增材制造中尚未充分探索的一个研究领域，这可能是未来几年的潜在研究方向。