选区熔化工艺的细观力学模拟研究现状

钱 亚，闫文韬，林 峰

（清华大学机械工程系，生物制造与快速成形技术北京市重点实验室，北京100084）

选区熔化技术作为金属增材制造的重要方法之一，可较高精度成形复杂结构，定制材料特性，近年来在追求轻质高强的航空航天领域得到了发展与应用[1]。其基本成形原理是利用一束聚焦热源按特定几何路径扫描预铺粉床，粉末吸收热量升温熔化并凝固成致密的结构；依次进行铺粉扫描，直至形成三维元件（图1）[2]。按照热源不同，可分为电子束选区熔化（EBSM）和激光选区熔化（SLM）。为使成形件达到足够的致密度与力学性能，选区熔化的能量输入要确保粉末充分熔化；另一方面，由于粉床几何的复杂性和能量空间分布不均，熔池内部分区域可能过热伴随蒸发而另外部分热量不足产生未熔合，导致层间缺陷的产生[3]。因而对成形工艺参数的合理选择将决定成形件的力学性能与质量。由于成形束斑的直径一般在几十到上百微米之间，粉末的升温熔化过程在微秒级别，通过实验动态观察选区熔化过程存在极大困难，借助数值模拟可重现成形过程并分析多参数对成形质量的影响。

早期对选区熔化的模拟工作集中在使用有限元方法求解粉床的热传导过程，根据温度信息可获得熔池的形状大小及热弹性关系求解应力分布[4]。这类工作的一大特征是将空间离散的粉床等效成连续均一化的实体，颗粒间的热导率、比热等参数也需综合粉床堆积密度、形状等进行相应的转化。尽管合理的简化可得到与实际较接近的计算结果，但这类模型存在先天的不足，使得难以深入捕捉选区熔化的过程：①粉末颗粒的形状大小分布在均一化模型中无法体现，这将可能导致球化产生和表面粗糙度的变化；②颗粒熔化后的自由流动、表面张力、蒸发反冲等因素共同作用，也将对熔池表面形貌产生重要影响。

图1 选区熔化过程示意图[2]

近年来，建立发展的分辨颗粒熔化流动的细观力学模型弥补了这些不足[5-6]。首先需解决的是对粉床的描述，以离散元方法（DEM）为基础，结合实际工艺中的材料、工艺参数，可得到具有物理意义的粉床分布[7]。利用离散元得到的粉床进行选区熔化的计算，除了热传导（能量）方程，还需考虑流体的动量，固液、气液边界，表面张力与反冲压力等因素[8]，一定程度上增加了模型的复杂性与计算量。

细观模型能精确捕捉成形过程中的孔隙与表面形貌的发展，为理解机理优化工艺提供了有力的工具。本文首先对细观模型计算研究现状进行总结，结合自身工作对部分技术细节、难点深入剖析，最后应用实例显示该模型对工艺指导的价值。

1 选区熔化细观模拟研究进展

1.1 粉末堆积与铺粉过程

作为选取扫描熔化过程的初始条件，粉床形成质量对后续工艺中粉末间的热传导、流动、团聚等现象有重要影响。一般为了获得密实的结构，粉床堆积以高相对密度和低表面粗糙度值为工艺目标。计算中，粉床的产生从空间内生成随机排布的颗粒开始，考虑重力作用及颗粒间接触逐步计算演化到稳定的堆叠模式。

Korner等采用的雨滴模型[5，9]和Zhou等采用的顺序堆叠算法[10]均是考虑了颗粒在重力作用下下落以达到能量最低状态，颗粒间存在滚动滑移但未显示计算相互之间作用力。这类方法易于产生大量的粉末颗粒的堆积分布，但由于并未充分考虑粉末颗粒的受力而缺乏物理可靠性，且应用仅限于堆积，无法模拟铺粉过程。

Xiang等采用离散元方法充分考虑了粉床堆积过程中的物理作用，包括赫兹接触的径向切向力、摩擦、范德华力，提高了模型的真实度，并研究了层厚及粉末颗粒大小分布对堆积效果的影响[11]。Parteli等借助DEM开源软件LIGGGHTS模拟了完整的铺粉过程，粉刷采用滚轮，对不同的铺粉速度、粉末颗粒的复杂形状进行了研究，且模型详细考量了各种物理模型和材料参数的选取，具有极高的仿真度和参考价值[7]。Haeri等在此基础上研究了粉刷形状、几何特征对柱状颗粒铺粉质量的影响[12]。

综合以上，对粉床堆积过程的模拟从初始的自由堆积到粉刷铺展，过程及物理作用逐步完善趋向于实际过程，但对工艺的优化探索，如粉刷形状、振动压实等影响有待研究。

1.2 选区扫描熔化过程

以离散化的粉末颗粒为对象，计算其受热熔化流动，这是细观模型的主要思想。但由于粉末几何边界的引入，以及对流体运动变化的描述（表面张力、蒸发反冲等），模型计算的实现增加了复杂度。

首先需解决的是如何准确添加热源描述电子束或激光束与粉床的作用关系。Yan采用Monte-Carlo方法计算电子束流与材料的相互作用，得到了电子束热源的能量空间分布[13]。Klassen根据实验和计算结果拟合出了电子束能量分布函数，可考虑加速电压、材料、入射角度等因素的影响[14]。为了处理电子束跟复杂表面的接触，可将束斑内的面能量源分解成若干子射线，每束射线对应特定空间坐标的几何信息[15]。电子束与材料的作用可穿透一定深度，相比之下，激光束与材料的作用穿透能力有限，但会在材料表面发生多次散射，为了准确计算激光束与粉床作用时的能量吸收，一般采用光线跟踪法（ray-tracing）来描述光线在粉床内部的散射吸收过程[16-17]。Boley等对激光束在粉床表面的吸收率进行了研究，证明了粉末颗粒的形状分布及光线在粉床内部的多次散射对吸收率的影响[18]。使用光线跟踪添加能量的输入能反映实际的作用过程，但由于粉床的熔化流动边界变化及光源移动，需要每个时间步计算造成了极大的计算量。

计算粉床的热传导过程，当温度高于材料的固态点（合金）或熔点时，需考虑熔化潜热，而当材料完全熔化还需求解流体的运动方程。多种计算流体力学（CFD）方法被用于粉末熔化流动过程的计算，如有限体积法（FVM）[19]、格子玻尔兹曼法（Lattice-Boltzmann Method）[20-21]和任意拉格朗日欧拉法（Arbitrarily Lagrange-Eulerian Method）[22]。 这些方法需具有同时处理固液问题的适应性，并对粉末流动时液面边界的自由运动具有很好的追踪能力，一般结合体积分数法（Volume of Fluid）来实现对边界的跟踪[23]。其他一些界面现象（如表面张力、反冲力等）则可通过在边界单元上施加指定大小的外力来等效[24-25]。除了束流能量输入、粉末间接触或流动传热，其他热流路径还有蒸发、辐射及传导对流，目前均在模型内得到相应的处理，但可靠性有待验证。仍待突破的是SLM工艺过程中成形室内气相的引入，实际工艺中产生的粉末飞溅及熔池的运动均可能受此影响，但同时处理固液气三相的难度较大。

德国学者Korner等最早开始细观粉末熔化过程的模拟[5]，使用自主研发的格子玻尔兹曼方法计算了二维电子束选区熔化过程，阐释了液体浸润及粉末随机分布造成的球化机理，并再现了不同能量密度下的单道扫描成形特征，与实验结果基本一致。随后，Korner等又在单层扫描的基础上研究了二维多层扫描工艺，对不同层厚、速度、能量对沿高度方向的影响进行分析，结果表明表面力和颗粒随机性是造成表面粗糙度和缺陷的原因[26]。

美国学者Khairallah等采用内部的ALE3D软件模拟了激光选区熔化的单道扫描过程，模型包含了表面张力、Marangoni流动、蒸发反冲等主要物理现象，并深入研究了剥蚀、飞溅、孔隙形成的微观机理（图2），提出了成形过程中表面力与反冲力的竞争作用关系[6，27]。

除了使用自行研发或内部软件，借助开源或商业计算流软件也实现了对细观选区熔化过程的模拟。Gurtler等用开源软件OpenFOAM计算了激光选区熔化，并揭示了能量不足时缺陷的形成，但模型未考虑颗粒的随机分布、蒸发反冲等重要因素[28]。Qiu等同样用OpenFOAM研究了扫描速度、粉层厚度对表面形貌及孔隙的影响，结果显示在高功率、低层厚下更易形成致密的结构[29]。Lee等用商业软件flow3d研究了激光选区熔化过程中的球化现象，提出了输入能量密度及粉床相对密度对球化的影响[19]。

借助熔化模拟得到的熔池温度信息，近年来开始尝试对凝固和组织生长过程进行模拟，且已取得相应的成果[30-31]，借助耦合的选区熔化-凝固模拟工具，后期可深入研究工艺与组织的关系，优化制件的性能。

图2 SLM工艺过程中的缺陷形成过程[27]

以离散化粉床熔化模拟为特色的细观力学方法能在粉末颗粒层面较真实地反映实际过程，相比于宏观均一化的方法可揭示表面形貌及缺陷的形成机理。但由于涉及到的物理过程复杂，且空间时间分辨率极高，在模型构建和计算效率上均有很大挑战。目前的工作大多局限在单道扫描过程，计算时间在几百甚至上千小时，距离面向实际工艺过程的CAE目标仍有很大差距。

2 多物理场模型构建

针对前文所述的多物理场模型，结合自身工作介绍了构建的DEM-CFD耦合计算电子束选区熔化（EBSM）全过程的模型框架，并结合实例深入探讨了多物理场对工艺过程的影响。

2.1 铺粉过程和影响因素

在用离散元方法计算铺粉的过程中，除了粉末接触挤压产生的径向切向力外，由于粉末颗粒尺寸较小且数量较多，颗粒间的滚动摩擦和粘附也会对粉床的流动性造成影响。模型中使用恒定方向扭矩模式施加滚动摩擦的作用[32]，其表达式为：

式中：ωrel为两接触颗粒相对角速度；μr为滚动摩擦系数；Rr为颗粒半径；Fn为压力。颗粒间的粘附力用简化的JKR模型描述，对相互接触的两颗粒在法向上添加吸附力：

式中：k为吸附能量密度；A为颗粒间接触面积，其计算式为：

式中：Δr为两颗粒中心间距；Ri、Rj分别为颗粒 i、j的半径。

为进一步说明摩擦与粘附对粉床堆积的影响，考虑如图3所示的铺粉过程，分别计算无摩擦无粘附、仅摩擦、仅粘附三种情况下粉床的相对密度，结果见表1。其中，粉末平均粒径为20 μm，粉床层厚为50 μm。

表1 滚动摩擦、粘附效应下粉床的相对密度

除物理因素外，几何条件也会对铺粉质量造成影响。借助离散元方法可比较不同粉刷形状与刮粉速度条件下得到的粉床密实度。从图3可看出，在未考虑滚动摩擦和粘附效应的情况下，速度增大，粉末飞溅角度倾向水平，粉刷与颗粒接触面越光滑，铺粉过程越稳定。

图3 不同情况下的粉床状态

2.2 电子束选区熔化过程

借助CAD软件建模，可将颗粒排布信息导入选区熔化计算模型。使用商业软件flow3d模拟粉床颗粒的熔化，其自带的VOF方法及FAVOR网格映射便于处理复杂几何边界[33]。电子束热源采用文献[13]中的分布函数添加到粉床上，通过将束斑离散成若干子束实现对粉床复杂表面的捕捉，以便计算入射角和实际受热区域（图4）。在选区熔化过程中，蒸发将减少熔池内表面材料的质量，同时带走一部分内能，而离开表面的蒸汽会对熔池表面产生反冲力作用[34]，因此蒸发过程对熔池的温度与形貌有重要影响。液体的蒸发量与温度的关系为：

式中：φ为材料蒸发系数；mA为原子质量；kB为玻尔兹曼常量；TS为熔池表面温度。

熔池表面的反冲压力可表示为：

式中：Lv为材料蒸发潜热；Tb为材料沸点。

图4 热源离散施加示意图[34]

根据上述蒸发和反冲压力公式，在每一个时间步内更新表面单元内的材料体积分数、内能与受力状态，以达到对蒸发过程的模拟。图5是蒸发反冲对熔池状态的影响，由于蒸发带走的热量极大降低了熔池的体积和熔道堆积高度，使成形面趋于平坦稳定。

图5 蒸发反冲对熔池状态的影响

2.3 铺粉-熔化一体化模型

图6是结合铺粉、选区熔化为一体的EBSM全过程模拟框架。以每一层为循环单元，先由DEM计算的铺粉作为熔化模拟的粉床的初始条件，再将CFD计算得到的熔道几何信息导入铺粉模型作为下一层铺粉的边界条件，如此重复即可将计算的规模拓展到多层[35]。

图6 EBSM全过程模拟框架[36]

3 工艺参数研究

依靠细观选区熔化模型，研究不同工艺参数对成形结果的影响，主要集中在熔池表面形貌的评估，部分结果如下：

（1）能量密度

定义能量密度为输入功率与扫描速度之比，其值直接关系到单位距离内粉床吸收的热量。图7是熔池在能量密度不足和过大时的两种特征形貌，粉末材料为Ti-6Al-4V。可看出，在能量密度不足的情况下熔道发生断裂、聚球，这是由表面张力在颗粒随机排布的粉床上作用的结果[5]；而能量密度过大时，熔道表面呈波浪形，这与蒸发反冲、Marangoni流动及Rayleigh不稳定性有关。计算结果与实验现象基本吻合[36]。

图7 不同能量密度下的熔池形貌

（2）速度

除了能量密度值决定粉床吸收的热量外，扫描速度大小对热量的耗散有决定作用。图8是在相同能量密度下，不同扫描速度对熔池形貌的影响。在相同的扫描距离和能量密度下，可认为粉床吸收了相同的能量。当扫描速度低时，热量有更多的时间向周围粉床传播，故有较大体积的粉床耗散这部分能量；当扫描速度高时，热量传播的时间和距离均有限，热作用区域体积较小（图8所示虚线内）。考虑极限情况，当速度趋于极小时，完成相同长度的扫描时间趋于无限长，颗粒吸收的热量被耗散而无法达到熔化的效果；当速度极大时，热量基本没有传导而锁定在束斑作用区域，粉床局部熔化乃至蒸发，但熔池体积有限，颗粒间及颗粒与基板或上一层难以充分熔合，这都不是理想的工艺过程。

图8 不同扫描速度下的熔池形貌

扫描速度过快伴随的另一个问题是熔池长宽比过大而导致的Rayleigh不稳定[37]，这将造成熔道宽度的波动（图8c），对相邻熔道间的成形控制有不利影响。

（3）颗粒分布

结合离散元方法可生成具有特定尺寸分布的粉床，继而研究粉床颗粒分布对成形效果的影响。如图9所示，分别对平均直径为70 μm和40 μm的粉床进行单道扫描计算，结果显示两种粒径分布的颗粒在相同工艺参数下熔道大小接近，颗粒小的粉床其熔池更细长，且熔道边缘相对平整。从图9可看出，颗粒较大的粉床堆积相对松散，内部有明显空隙，而颗粒较小的粉床排布更密实。可以认为，实际的熔化受双重因素的影响：粉床密度与颗粒大小。对颗粒较大的粉床其相对密度不足，导致颗粒间传热率降低，但大颗粒内部比小颗粒的导热效率更高，因而两者共同作用下表现出的不同粉床的熔化效果差距并不明显。

图9 不同粉末颗粒分布下的熔池形貌

4 结束语

运用细观模型研究选区熔化工艺为理解复杂过程提供了有力的工具。其特点是以粉末为基本单元，从铺粉和扫描两个环节分别用离散元法计算粉末颗粒的碰撞运动，以及用计算流体力学法计算粉末的传热熔化过程。模拟中涉及到的多物理参数和现象均对结果有不可忽略的影响，需谨慎处理。在该模型基础上，结合实际工艺分析熔池形貌的演化，从结果可看出，由能量密度决定的粉床吸收热量直接影响熔池的大小和形貌，扫描速度和粉末尺寸与热量的传导耗散有关，并间接影响相邻熔道与相邻层间的结合。细观模型在处理几何随机性、解释形貌和缺陷机理上有独特的价值，但巨大的计算量也制约了其发挥更广泛的作用，如何借助目前有限的计算规模去涵盖复杂的工艺变量和流程，需对理论机制或模型算法进一步完善。

[1] 巩水利，锁红波，李怀学.金属增材制造技术在航空领域的发展与应用[J].航空制造技术，2013（13）：66-71.

[2] KÖRNER C. Additive manufacturing of metallic components by selective electron beam melting-a review[J].International Materials Reviews，2016，61 （5）：361-377.

[3] JUECHTER V，SCHAROWSKY T，SINGER R F，et al.Processing window and evaporation phenomena for Ti-6Al-4V produced by selective electron beam melting[J].Acta Materialia，2014，76：252-258.

[4] SCHOINOCHORITISB，CHANTZISD，SALONITISK.Simulation of metallic powder bed additive manufacturing processes with the finite element method:A critical review [J].Proceedings of the Institution of Mechanical Engineers，Part B:Journal of Engineering Manufacture，2015，231（1）：96-117.

[5] KÖRNER C，ATTAR E，HEINL P.Mesoscopic simulation of selective beam melting processes [J].Journal of Materials Processing Technology，2011，211（6）：978-987.

[6] KHAIRALLAH S A，ANDERSON A. Mesoscopic simulation model of selective laser melting of stainless steel powder [J]. Journal of Materials Processing Technology，2014，214（11）：2627-2636.

[7] PARTELI E J R，PÖSCHEL T,Particle-based simulation ofpowderapplication in additive manufacturing[J].Powder Technology，2016，288：96-102.

[8] MARKL M，KORNER C.Multiscale modeling of powder bed-based additive manufacturing[J].Annual Review of Materials Research，2016，46：93-123.

[9] MEAKIN P，JULLIEN R.Restructuring effects in the rain model for random deposition[J].Journal de Physique，1987，48（10）：1651-1662.

[10]ZHOU J，ZHANG Y，CHEN J K.Numerical simulation of random packing of spherical particles for powder-based additive manufacturing [J].Journal of Manufacturing Science and Engineering，2009，131（3）：031004.

[11]XIANG Zhaowei，YIN Ming，DENG Zhenbo，et al.Simulation of forming process of powder bed for additive manufacturing[J].Journal of Manufacturing Science and Engineering，2016，138（8）：081002.

[12]HAERI S，WANG Y，GHITA O，et al.Discrete element simulation and experimental study of powder spreading process in additive manufacturing[J].Powder Technology，2017，306：45-54.

[13]YAN Wentao，SMITH J，GE Wenjun，et al.Multiscale modeling of electron beam and substrate interaction：a new heat source model[J].Computational Mechanics，2015，56（2）：265-276.

[14]KLASSEN A， BAUEREIβ A， KÖRNER C.Modelling of electron beam absorption in complex geometries[J].Journal of Physics D：Applied Physics，2014，47 （6）：065307.

[15]ZOHDI T I.Rapid simulation of laser processing of discrete particulate materials [J]. Archives of Computational Methods in Engineering，2013，20 （4）：309-325.

[16]CHO D W，CHO W I，NA S J.Modeling and simulation of arc:laser and hybrid welding process[J].Journal of Manufacturing Processes，2014，16（1）：26-55.

[17]AHN J，NA S J.Three-dimensional thermal simulation of nanosecond laser ablation for semitransparent material[J].Applied Surface Science，2013，283：115-127.

[18]BOLEY C D，KHAIRALLAH S A，RUBENCHIK A M.Calculation of laser absorption by metal powders in additive manufacturing[J].Applied Optics，2015，54（9）：2477-2482.

[19]LEE Y S，ZHANG W.Mesoscopic simulation of heat transfer and fluid flow in laser powder bed additive manufacturing[J].Solid Freeform Fabrication Symposium，2015：1154-1165.

[20]KÖRNER C，THIESM，HOFMANN T，etal.Lattice boltzmann model for free surface flow for modeling foaming[J].Journal of Statistical Physics，2005，121 （1-2）：179-196.

[21]ATTAR E，KÖRNER C.Lattice boltzmann model for thermalfree surface flowswith liquid-solid phase transition [J].International Journal of Heat and Fluid Flow，2011，32（1）：156-163.

[22]MCCALLEN C.ALE3D:arbitrary lagrangeeulerian three-and two dimensional modeling and simulation capability [R].Livermore：Lawrence Livemore National Laboratory，2012.

[23]HIRT C W，NICHOLS B D.Volume of fluid(VOF)method for the dynamics of free boundaries[J].Journal of Computational Physics，1981，39（1）：201-225.

[24]KLASSEN A，SCHAROWSKY T，KÖRNER C.Evaporation model for beam based additive manufacturing using free surface lattice Boltzmann methods[J].Journal of Physics D：Applied Physics，2014，47（27）：275303.

[25]ATTAR E，KÖRNER C.Lattice boltzmann method for dynamic wetting problems[J].Journal of Colloid and Interface Science，2009，335（1）：84-93.

[26]KÖRNER C，BAUEREIβ A，ATTAR E.Fundamental consolidation mechanisms during selective beam melting of powders[J].Modelling and Simulation in Materials Science and Engineering，2013，21（8）：085011.

[27]KHAIRALLAH S A，ANDRSON A T，RUBENCHIK A，et al.Laser powder-bed fusion additive manufacturing:Physics of complex melt flow and formation mechanisms ofpores,spatter,and denudation zones [J].Acta Materialia，2016，108：36-45.

[28]GÜRTLER F J，KARG M，LEITZ K H，et al.Simulation of laser beam melting of steel powders using the threedimensional volume of fluid method[J].Physics Procedia，2013，41：881-886.

[29]QIU C，PANWISAWASET C，WARD M，et al.On the role of melt flow into the surface structure and porosity development during selective laser melting [J].Acta Materialia，2015，96：72-79.

[30]RAI A，MARKL M，KÖRNER C.A coupled cellular automaton-lattice boltzmann model for grain structure simulation during additive manufacturing [J].Computational Materials Science，2016，124：37-48.

[31]PANWISAWAS C，QIU Chunlei，ANDERSON M J，et al.Mesoscale modelling of selective laser melting:thermal fluid dynamics and microstructural evolution [J].Computational Materials Science，2017，126：479-490.

[32]AI Jun，CHEN Jianfei，ROTTER J M，et al.Assessment of rolling resistance models in discrete element simulations[J].Powder Technology，2011，206（3）：269-282.

[33]HIRT C W，NICHOLS B D.Flow-3D User's Manual[M].Leland：Flow Science Inc.，1988.

[34]ALLMEN M V.Laser-beam interactions with materials:physical principles and applications [J]. Nuclear Instruments and Methods in Physics Research，1997，124（4）：647-648.

[35]闫文韬，钱亚，林峰.选区熔化过程多尺度多物理场建模研究进展[J].航空制造技术，2017（10）：50-58.

[36]GUO Chao，GE Wenjun，LIN Feng.Effects of scanning parameters on material deposition during electron beam selective melting of Ti-6Al-4V powder[J].Journal of Materials Processing Technology，2015，217：148-157.

[37]GUSAROV A，YADROITSEV I，BERTRAND P，et al.Heat transfer modelling and stability analysis of selective laser melting[J].Applied Surface Science，2007，254（4）：975-979.