激光选区熔化铺粉过程的数值模拟及粉层表征1)

刘海林 易 敏 王建祥 易 新†,2)

* (北京大学工学院力学与工程科学系,湍流与复杂系统国家重点实验室,北京 100871)

† (中国兵器工业导航与控制技术研究所,北京 100089)

** (南京航空航天大学航空学院,机械结构力学及控制国家重点实验室,南京 210016)

†† (北京大学南昌创新研究院,南昌 330096)

引言

激光选区熔化技术逐层制造的成形方式和快速熔化凝固的热循环特点使其可以调控零件微观组织演化,生产性能优异的复杂构型零件,在制备高熵合金、高温合金和难熔金属等构件方面具有广泛应用[1-5].优异的成形零件质量是激光选区熔化技术得以广泛应用的基础.成形零件质量受粉层质量、熔池流动、微观组织演化和热应力分布的影响[6].目前研究人员分别通过模拟铺粉过程[7-8]、细观尺度传热传质现象[9-10]、微观尺度晶粒生长行为[11-12]和宏观尺度零件成形过程[13]来对上述影响因素进行研究.

良好的金属粉体颗粒流动性是粉层金属颗粒紧密均匀排布的前提,也是激光选区熔化成形性能优异零件的基础条件之一.颗粒流动性常用休止角来反映,流动性越好,休止角越小.金属激光选区熔化粉体一般为尺寸几十μm 的金属颗粒.此时颗粒间黏附力与颗粒表面能密度呈正比,并对颗粒流动性具有重要影响.离散元模拟指出,金属颗粒间黏附力越大,休止角越大,颗粒流动性越差[14].此外,金属颗粒的形貌和尺寸分布等也对颗粒流动性有着重要影响.由于摩擦效应和互锁效应,形貌不规则、表面粗糙的金属颗粒往往有着较差的流动性[15].对于同时包含小粒径和大粒径金属颗粒的粉层,小粒径金属颗粒在小体积分数情况下,会倾向于填充大粒径颗粒间的缝隙,而当小粒径颗粒体积分数较大时,其容易团聚成块.该团块会显著阻碍颗粒流动,降低混合金属颗粒的流动性,增加粉层中空白区域的数量[16].

除颗粒流动性外,粉层厚度和铺粉速率也对粉层质量有重要影响.Yao 等[17]以粉层密度和粉层密度标准差作为衡量粉层质量的指标,使用离散元方法模拟了316L 不锈钢粉末的铺粉过程,研究了粉层厚度和铺粉速率对粉层质量的影响.研究发现,高铺粉速率会增加颗粒流动的不稳定性,低粉层厚度会增加壁面效应和颗粒堵塞概率,高铺粉速率和低粉层厚度都会降低粉层质量[17].Chen 等[18]以粉层密度和粉层表面粗糙度为指标衡量粉层质量,使用数值和实验方法研究了铺粉速度对粉层质量的影响.研究发现,增大铺粉速度会降低粉层质量.粉层质量的含义是多方面的,如粉层密度、配位数分布、粉层表面粗糙度和粉层均匀性等.确定合适的粉层厚度和铺粉速率需要对粉层质量进行全面的表征.然而,上述研究仅选择其中若干个指标研究了粉层厚度和铺粉速率对粉层质量的影响.

多层铺粉过程中,成形区的颗粒流动行为不仅受颗粒流动性影响,还受成形区表面形貌影响.Parteli等[19]使用组合球方法构建具有一定粗糙度的成形区表面,研究了多层铺粉过程中成形区表面承受载荷的时空演化,发现颗粒内部的力链演化会导致成形区表面载荷的大幅波动.孙远远等[7]以单道为基本单元建立粗糙表面模型,研究了铺粉方向、单道夹角、单道搭接率和粉层厚度对粉层质量的影响,研究发现成形区的粗糙表面有利于提高粉层质量,适当地提高粉层厚度与优化铺粉角度可以增加成形区粗糙表面上的粉层均匀性.Nan 等[20]使用搭接率为50%的平行圆柱表示成形区粗糙表面,并在此基础上研究了颗粒-保护气相互作用对粉堆内部颗粒流动行为和粉层均匀性的影响,发现保护气体的拖拽力减缓了粉堆内部颗粒的对流循环速率,增大了粉堆休止角,恶化了粉层质量.Yao 等[21]使用矩形方块表示成形零件,研究了颗粒-成形零件间协同作用对颗粒行为的影响,并分析了影响机理.然而,上述研究使用简单几何体或组合几何体表示的成形区粗糙表面十分具有规律性,是对真实成形区粗糙表面的简化,与实际情况有较大偏离.

目前,在粉层厚度和铺粉速率影响粉层质量的研究中,对粉层质量的表征还不全面.此外,构建成形区粗糙表面的现有方法不能准确刻画单道熔迹凹凸变化的真实粗糙表面,也无法表示部分熔融金属颗粒.本文基于离散元模拟方法,获得了金属钨颗粒的滚动摩擦系数和表面能密度,分析了它们对颗粒流动性的影响,定量研究了铺粉速率和粉层厚度对粉层质量的影响,并构建了高效多层铺粉离散元计算新模型.

1 粉层质量表征指标

粉层质量包含粉层密度、配位数分布、粉层表面粗糙度和粉层均匀性等.然而,目前尚无评定表征粉层质量的国内和国际标准,相关文献中的粉层质量表征指标也不统一.表1 总结了文献中提到的粉层质量评定指标及其物理含义.

2 离散元模拟的基本框架

本文将建立获取金属颗粒性质的离散元计算模型、单层铺粉离散元模型和多层铺粉离散元模型,这些离散元模型具有相似的基本框架.本节将集中介绍这些离散元模型的基本框架.

2.1 颗粒接触模型与黏附模型

在离散元模拟中,颗粒的运动受到重力、颗粒间的接触力、接触力矩和范德华作用导致的黏附力共同影响.颗粒间的接触力可用Hertz-Mindlin 模型描述[42-44].考察相互接触的颗粒i与颗粒j(图1),颗粒半径分别为Ri和Rj,质量分别为mi和mj,角速率分别为ωi和ωj.

图1 颗粒i 与颗粒j 接触示意图Fig.1 Schematic of contact between particle i and particle j

此时,颗粒j作用于颗粒i的法向接触力Fn和切向接触力Ft分别为[32,45]

式中,下标n 和t 分别表示法向和切向的物理量,en(t)为法向 (切向) 单位向量,δn=Ri+Rj-Rij为法向重叠量,δt为切向相对位移,为δn(t)的时间导数,Kn(t)为法向 (切向) 刚度系数,γn(t)为法向 (切向) 阻尼系数.

法向刚度系数Kn的表达式为[32]

其中,E为颗粒杨氏模量,v为颗粒泊松比,有效颗粒半径R*满足1/R*=1/Ri+1/Rj.

切向刚度系数Kt的表达式为[32]

其中,有效剪切模量G*=E/[4(2+v)(1-v)].

法向阻尼系数γn的表达式为[32]

切向阻尼系数γt的表达式为[32]

切向接触力满足|Ft|≤μs|Fn|,μs为静摩擦系数,有时也称库伦摩擦系数[32].

此外,颗粒j作用于颗粒i的接触力矩为[35]

其中,μr为滚动摩擦系数.

由范德华作用导致的颗粒间黏附由经典JKR黏附模型描述,黏附力FJKR表达式为[47]

其中,γ是颗粒表面能密度,a是颗粒间接触区域的半径.a和δn满足如下方程[19]

在离散元模型中,设置基板、刮刀和颗粒的材料与物性参数相同.当颗粒与基板或刮刀发生相互作用时,视基板或刮刀为半径和质量都无穷大的颗粒.

2.2 利用瑞利波法确定时间步长

实验发现,当金属颗粒接触碰撞时,瑞利波消耗的能量占总能耗的70%[48].因此,可以通过沿球形金属颗粒表面传播的瑞利波确定离散元模型的临界时间步长.

瑞利波方程为[48]

瑞利波方程根βR的近似表达式为[48]

金属颗粒表面的瑞利波波速vR为[48]

颗粒间的接触作用只影响发生碰撞的两个金属颗粒,而不应该通过瑞利波的传播影响其他金属颗粒.因此,为了保证计算的稳定性与精度,离散元模拟的时间步长Δt需小于瑞利波传递半个金属颗粒球面所需要的时间TR,即

其中,Rmin取为发生碰撞的两个金属颗粒的半径最小值.这里设置时间步长为Δt=30%×(TR)min,时间积分格式选择 velocity Verlet 格式[14,30].

2.3 放大离散元模拟的时间步长

图2(a) 展示了在后文相关实验中使用的球形钨颗粒形貌,图2(b) 为相应的钨颗粒粒径分布.

图2 球形钨粉颗粒的形貌与粒径分布Fig.2 Morphology and size distribution of tungsten powder

如图2(b) 所示,钨金属颗粒的尺寸为5～40 μm,即Rmin=2.5 μm.钨的密度ρ=1.935×104kg/m3,杨氏模量E=340 GPa,泊松比v=0.28[49-50].此时离散元模型的时间步长Δt=9.75×10-10s 过小,计算量过大,必须考虑适当地放大时间步长.

从式(13)可知,放大时间步长的方法包括增大最小金属颗粒半径Rmin,增大钨密度ρ和减小钨杨氏模量E.通过增大钨密度ρ来放大时间步长的方法适用于重力可以忽略的准静态模拟.然而,铺粉是一个动态过程且重力不可忽略.这里通过忽略部分小体积钨颗粒与减小钨杨氏模量来适当放大时间步长.

为了减小最小金属颗粒半径Rmin对时间步长的限制,在离散元模拟中忽略半径为2.5～6 μm 的钨颗粒.此时,离散元模拟中的钨颗粒粒径分布如图3所示.

图3 离散元模拟中,钨颗粒质量百分数和累计质量分数随颗粒直径的变化Fig.3 Mass percent and cumulative mass percent versus tungsten powder diameter

大量的数值研究表明,在离散元模拟中减小杨氏模量对颗粒行为与模拟结果影响不大[14,30,51].因此,在离散元模拟中将钨的杨氏模量减小两个数量级至E′=3.4 GPa,此时Δt=2.38×10-8s.

当颗粒杨氏模量减小时,颗粒刚度会随之减小.当颗粒刚度减小,而颗粒表面能密度保持不变时,颗粒将表现得更加具有“黏附性”,颗粒整体行为将发生改变[52].因此,为了使刚度减小后的颗粒整体行为不变,按照Hærvig 等[53]与Washino 等[54]提出的准则缩减钨的表面能密度,即

式中,γ为钨的表面能密度,γ′为钨的缩减表面能密度.

3 确定颗粒静摩擦系数、滚动摩擦系数与表面能密度

目前,离散元模拟中已知的物性参数包括钨的密度ρ,缩减杨氏模量E′和泊松比v.尚待确定的物性参数包括颗粒静摩擦系数μs、滚动摩擦系数μr和表面能密度γ.通过斜面仪来测定静摩擦系数μs,对比模拟与实验得到的休止角和致密度来确定滚动摩擦系数μr和表面能密度γ.需要特别指出,本节得到的滚动摩擦系数和表面能密度可以使离散元模型中的钨颗粒整体表现与真实颗粒整体表现相同,但并不与真实钨金属颗粒的滚动摩擦系数和表面能密度相同.这是因为颗粒形貌和尺寸分布对颗粒整体表现也有重要影响,而离散元模拟中颗粒尺寸分布与真实颗粒尺寸分布存在一定差异,且真实的钨颗粒形貌不是模拟所采用的完美球形.

3.1 测量钨颗粒间静摩擦系数

采用图4 所示的自制斜面仪来确定钨颗粒间静摩擦系数.将底部用双面胶沾满球形钨颗粒的滑块被放置于斜置钨板上,缓慢增大斜置钨板的倾斜角度,记录滑块发生滑动时的钨板倾斜角度αc.实验测得临界倾角αc=15°,得到静摩擦系数μs为

图4 静摩擦系数μs 测量实验示意图Fig.4 Schematic of measuring coefficient μs of static friction

3.2 测量休止角和致密度

采用漏斗法测量钨颗粒的休止角.图5 展示了漏斗法测量休止角的实验示意图,其中,D1=50.75 mm,D2=2.68 mm,d=30.00 mm,H1=34.96 mm,H2=9.96 mm,H3=46.14 mm.

图5 漏斗法测休止角α 的实验示意图Fig.5 Schematic of measuring repose angle α using the funnel method

实验结果如图6 所示,其中D=32.48 mm,H为粉堆高度,α=arctan(2H/D)为钨粉堆休止角.这里,D包含金属圆筒壁面厚度,而图5 中的d不包含金属圆筒壁面厚度.

图6 漏斗法测休止角α 的实验结果Fig.6 Experimental result of measuring repose angle α using the funnel method

图7 漏斗法测量休止角α 的离散元计算模型Fig.7 Schematic of discrete element modeling of measuring repose angle α using the funnel method

表2 列出了利用漏斗法测量的休止角结果,根据这些数据,可估算出钨粉堆休止角为

表2 钨粉堆休止角α 测量结果Table 2 Measured repose angles of tungsten powder piles

休止角测量结束后,将超出圆筒顶端的钨粉抹去,测量钨粉的松装密度.测得空桶质量为 67.0 g,空桶与剩余钨粉总质量为 340.5 g,空桶体积为πd2H1/4=24 710.9 mm3,得到钨粉松装密度为11.068 g/cm3,致密度ρr为

3.3 离散元模拟确定滚动摩擦系数与表面能密度

在休止角和致密度测量中使用的钨颗粒数量巨大 (约10 亿颗),离散元模拟很难计算如此数量的颗粒运动,缩小离散元模型的几何尺寸与颗粒数量是减小离散元模拟计算量的常用方法[16-17,51].此外,根据 Meier 等[14]的研究结果,缩小离散元模拟的粉堆几何尺寸对得到的参数无明显影响.因此,本小节通过缩小离散元模拟的粉堆几何尺寸与颗粒数量来提升计算效率.

图8 为测量钨颗粒致密度ρr的离散元模型,其中,长方体尺寸为a=0.5 mm,h=2 mm,长方体水平方向设置为周期边界条件.在离散元模拟的初始时刻,在长方体中生成球形钨颗粒约3.2 万个.待钨颗粒在重力的作用下完成沉积静置时,取粉堆中间部分测量粉堆致密度.

图8 测粉堆致密度ρr 的离散元计算模型Fig.8 Schematic of measuring relative density ρr of the powder pile using discrete element modeling

在测休止角和致密度的离散元模拟中,设置滚动摩擦系数μr的取值分别为

设置缩减表面能密度γ′的取值分别为

图9 绘制了休止角α随滚动摩擦系数μr和缩减表面能密度γ′的变化.如图9 所示,休止角α随着μr和γ′的增大而增大,即钨颗粒的流动性随滚动摩擦系数μr和缩减表面能密度γ′的增大而减小.图9 (b) 的模拟结果与已报道的实验结果一致[14].此外,休止角α对γ′变化的敏感性随γ′的增大而迅速减小.

图9 休止角α 随滚动摩擦系数μr 和缩减表面能密度γ′的变化Fig.9 Repose angle α versus coefficient μr of rolling friction and reduced surface energy density γ′

图10 绘制了致密度ρr随滚动摩擦系数μr和缩减表面能密度γ′的变化.如图10 所示,致密度ρr随着μr和γ′的增大而减小.

图10 粉堆致密度ρr 随μr 和γ′的变化Fig.10 Relative density ρr of the powder pile versus coefficient μr of rolling friction and reduced surface energy density γ′

图11 给出了致密度ρr随滚动摩擦系数μr和缩减表面能密度γ′变化的云图.云图上叠加了曲线ρr=0.572 (虚线) 和α=32.81° (实线).虚线ρr=0.572 和实线α=32.81°的交点所对应的滚动摩擦系数μr=0.031 5 和缩减表面能密度γ′=0.15 mJ/m2即为通过实验校准得到模拟用μr和γ′.

图11 粉堆致密度ρr 随滚动摩擦系数μr 和缩减表面能密度γ′变化的云图Fig.11 Contour plot of relative density ρr of the powder pile varying with coefficient μr of rolling friction and reduced surface energy density γ′

至此,离散元模拟需要的所有参数已全部确定并汇总于表3.

表3 铺粉离散元模拟所需参数Table 3 Parameters for discrete element modeling of powder spreading

4 单层铺粉离散元模拟

4.1 单层铺粉离散元计算模型

图12 给出了单层铺粉的离散元计算模型,包括6 mm×0.4 mm 的基板,0.1 mm×0.4 mm×1 mm 的刮刀和 0.6 mm×0.4 mm×1 mm 的储粉器.刮刀与基板间的距离为粉层厚度δ.在离散元模拟的初始时刻,在储粉器中生成约2.7 万个钨颗粒,并在重力作用下沉积静置.然后,移除储粉器右侧壁面,与此同时位于储粉器左侧的刮刀开始向右以速率V移动,进行铺粉.

图12 单层粉层离散元计算模型 (单位: mm)Fig.12 Discrete element modeling of single-layer powder bed spreading (unit: mm)

本节研究铺粉速率V和粉层厚度δ对粉层质量的影响.设置V的取值为

设置粉层厚度δ的取值为

为了减少粉层起始端x=0 mm 和终止端x=6 mm 对粉层质量的影响,待铺粉终止且颗粒静止后,选择x∈[2.75,5.25] mm 内的粉层作为观测区进行粉层质量表征.

4.2 粉层质量表征

图13(a) 展示了铺粉速率V和粉层厚度δ对粉层形貌的影响.从图13 可以发现,增大铺粉速率V和减小粉层厚度δ都会减少铺粉过程中沉积的颗粒数量,增加粉层的空斑缺陷.减小粉层厚度会增加粉层空斑缺陷的模拟结果与图13(b) 所示的已报道实验结果吻合良好[39].

图14 展示了铺粉速率V和粉层厚度δ对粉层密度的影响.这里粉层密度表示粉层内金属颗粒总体积与粉层总体积的比值,即

图14 粉层密度 随粉层厚度δ 与铺粉速率V 的变化Fig.14 Packing density of the powder layer versus layer thickness δ and spreading speed V

这里,N表示观测区内金属颗粒的数量,Vi表示观测区粉层内第i个金属颗粒的体积,L,W和δ分别表示观测区的长、宽和高.

配位数分布是用来反映粉层内颗粒排布紧凑程度的常用指标.对任意给定的一个金属颗粒,其配位数为与其接触的其余颗粒数量.一般来说,配位数越大,粉层内颗粒排布越紧密,打印时出现缺陷的概率越小[32].图15 给出了铺粉速率V和粉层厚度δ对颗粒配位数分布的影响.如图15(a) 所示,对于任意给定的铺粉速率V,出现概率最大的配位数随着粉层厚度δ的增大而有所增大,且增大的幅度随着铺粉速率V的减小而增大.如图15(b) 所示,对于任意给定的粉层厚度δ,出现概率最大的配位数随着铺粉速率V的减小而有所增大,且增大的幅度随着粉层厚度δ的增大而增大.当粉层厚度δ=25 μm,铺粉速率V=0.2 m/s 时,金属颗粒配位数分布的峰值为0,表示大部分颗粒都不与其他颗粒接触.当粉层厚度δ=75 μm,铺粉速率V=0.02 m/s 时,金属颗粒配位数分布的峰值为3,颗粒排布变得更加紧密.

图15 粉层厚度δ 与铺粉速率V 对配位数分布的影响Fig.15 Effects of the layer thickness δ and spreading speed V on the distribution of coordination number

在激光选区熔化中,粉层的均匀性通常指不同尺寸金属颗粒分布的随机性[25].粉层的均匀性对打印时的能量传输和熔池形貌都有重要影响.这里使用三维Delaunay 四面体来表征粉层的均匀性.Delaunay 四面体是指连接粉层内临近4 个金属颗粒得到的四面体,每个四面体内除顶点处的金属颗粒,不包含其他金属颗粒.在一个均匀的粉层中,较小的金属颗粒应均匀地分布在较大的金属颗粒间,使Delaunay 四面体具有较窄的等效半径分布[25].这里的等效半径是指与Delaunay 四面体体积相同的球的半径.图16 展示了Delaunay 四面体等效半径分布随粉层厚度δ与铺粉速率V的变化.如图16(a) 所示,当铺粉速率V≤ 0.08 m/s 时,增大粉层厚度δ可以使四面体等效半径的分布曲线更陡峭,显著提高粉层均匀性;当V＞ 0.08 m/s 时,δ对四面体等效半径的分布基本无影响.如图16(b) 所示,当δ≥ 45 μm时,减小V可以使四面体等效半径的分布曲线更陡峭,显著提高粉层均匀性;当δ＜ 45 μm 时,V对四面体等效半径的分布基本无影响.

图16 Delaunay 四面体等效半径分布随 (a) 粉层厚度δ 与 (b) 铺粉速率V 的变化Fig.16 Equivalent radius of Delaunay tetrahedrons versus (a) layer thickness δ and (b) spreading speed V

图17 展示了铺粉速率V和粉层厚度δ对Delaunay四面体球形度分布的影响.Delaunay 四面体球形度是指与 Delaunay 四面体体积相同的球的面积与该四面体的面积的比值.当Delaunay 四面体球形度趋近于0.671 时[41],Delaunay 四面体趋近于正四面体,粉层内金属颗粒相互堆叠.当Delaunay 四面体球形度趋近于0 时,Delaunay 四面体扁平,粉层内金属颗粒基本没有堆叠.如图17(a) 所示,当铺粉速率V≤0.08 m/s 时,增大粉层厚度δ可以使金属颗粒大量堆叠;当V＞ 0.08 m/s 时,δ对金属颗粒的堆叠无显著影响.如图17(b) 所示,当δ≥ 45 μm 时,减小V可以使金属颗粒大量堆叠;当δ＜ 45 μm 时,V对金属颗粒的堆叠无显著影响.

图17 Delaunay 四面体球形度分布随 (a) 粉层厚度δ 与 (b) 铺粉速率V 的变化Fig.17 Sphericity of Delaunay tetrahedrons versus (a) layer thickness δ and (b) spreading speed V

综合考虑铺粉速率V和粉层厚度δ对配位数分布、Delaunay 四面体等效半径分布和球形度分布的影响,可以发现铺粉速率V≤ 0.08 m/s,粉层厚度δ≥ 45 μm,是粉层颗粒有效堆叠与紧密均匀排布的必要条件.

已有的研究表明,粗糙的粉层表面会增加成形零件的孔隙率与表面粗糙度[56-58].这里使用射线追踪法确定粉层的表面轮廓,进而计算粉层表面粗糙度.如图18 所示,基板表面位于z=0,沿-z方向传播的射线与粉层表面相交时,记录交点高度z.观测区内激光射线沿x和y方向均匀分布,间距为2 μm,共采集表面轮廓高度数据1250×200 个.取粉层表面轮廓高度的标准差作为粉层表面粗糙度Sa,即

图18 射线追踪法确定粉层高度示意图Fig.18 Schematic of determining the height of the powder layer using ray tracing approach

图19 展示了粉层厚度δ和铺粉速率V对粉层表面粗糙度Sa的影响.如图19 所示,粉层表面粗糙度Sa随粉层厚度δ的增大而增大,随扫描速率V的增大而减小.通常情况下,粉层表面粗糙度Sa随粉层厚度δ和铺粉速率V变化的趋势受刮刀类型、刮刀形状和基板粗糙度等众多因素的影响.因此,不同的铺粉工况可能造成不同的变化趋势.例如,Zhao等[59]基于离散元模拟研究了刮刀在水平基板上移动的铺粉过程,发现Sa随δ的增大而增大,与本文的结论一致.Xiang 等[36]模拟了刮刀在粗糙基板上的铺粉过程,发现Sa随着V的增大先增大后减小,随粉层厚度δ的增大而增大.Haeri[60]基于离散元模拟研究了形状优化后的刮刀在水平基板上进行铺粉的过程,发现Sa随V的增大而增大.

图19 粉层表面粗糙度Sa 随粉层厚度δ 与铺粉速率V 的变化Fig.19 Surface roughness Sa of the powder layer versus layer thickness δ and spreading speed V

粉层颗粒偏析是铺粉过程中的一种常见现象,特别是当所使用的金属粉颗尺寸分布范围较大时[15,31,34].颗粒偏析改变了粉末层中颗粒尺寸的空间分布,直接影响了选区熔化过程中激光-颗粒的相互作用以及随后的熔化凝固过程[61].为了分析颗粒偏析,将观测区内沉积的粉层进行网格划分.在x方向划分为 10 格,y和z方向不进行划分,x方向网格尺寸为 250 μm,y方向网格尺寸为400 μm,z方向网格尺寸为粉层厚度δ.定义网格偏析指数SI为该网格中金属颗粒的平均半径与初始时刻离散元计算模型中所有金属颗粒的平均半径之比.偏析指数SI表达式如下[34]

其中,n表示网格中金属颗粒数量,Ri表示网格中第i个金属颗粒的半径,表示初始时刻离散元计算模型中所有金属颗粒的平均半径.当SI＜ 1 时,小颗粒更容易在基板上沉积;当SI＞ 1 时,大颗粒更容易在基板上沉积.

图20 展示了粉层厚度δ和铺粉速率V对偏析指数SI分布的影响.如图20(a) 所示,增大铺粉速率V可以显著降低粉层的颗粒偏析现象,与已有的研究结果一致[30-31].在刮刀铺粉过程中,小尺寸金属颗粒相较于大尺寸金属颗粒具有更大的沉积速率,更容易在粉堆中下降并沉积到基板上,形成颗粒偏析.然而,铺粉速率V的增加减小了小尺寸颗粒在粉堆中下降并先于大尺寸颗粒沉积的时间窗口,在一定程度上减弱了颗粒偏析[30-31].如图20(b) 所示,增大粉层厚度δ同样降低了粉层的颗粒偏析现象.在铺粉过程中,较小的粉层厚度δ容易导致颗粒堵塞.在颗粒堵塞期间,刮刀与基板间狭缝入口前的粉堆中存在强接触力拱,导致粉堆中金属颗粒无法进入狭缝.随着刮刀的前进,颗粒堵塞状态被打破,强接触力拱中存储的弹性能释放,狭缝前的金属颗粒穿过狭缝并在基板上沉积.小尺寸金属颗粒质量小,在强接触力拱坍塌后,更容易穿过狭缝,并最终造成颗粒偏析 (SI＜ 1)[62].此外,粉层厚度δ越小,颗粒堵塞发生的频率与持续时间越长.因此,当粉层厚度δ=25 μm 时,SI＜ 1 (除V=0.2 m/s 之外),颗粒偏析现象明显.随着粉层厚度δ增大,偏析现象迅速衰减消失.

图20 粉层偏析指数SI 分布随粉层厚度δ 与铺粉速率V 的变化Fig.20 Index of particle segregation in powder layer versus layer thickness δ and spreading speed V

图20 粉层偏析指数SI 分布随粉层厚度δ 与铺粉速率V 的变化 (续)Fig.20 Index of particle segregation in powder layer versus layer thickness δ and spreading speed V (continued)

5 高效多层铺粉离散元新模型的构建

如图21(a)和图21(b) 所示,通过3D 模型stl 文件,可以将激光选区熔化传热传质模拟中已部分熔化的粉层的全部几何形貌信息(图21(a))作为固定边界条件,输出至多层铺粉离散元计算模型中(图21(b)).此时,在如图21(b)所示的多层铺粉离散元模型中,未熔化金属颗粒、部分熔化金属颗粒和熔迹的几何形貌都由大量三角面片拼接而成.例如图21(c) 中的球形金属颗粒,即是由约1000 个三角面片拼成.这些构成未熔化球形金属颗粒的大量三角面片极大地降低了多层铺粉离散元模型的计算效率.

图21 高效多层铺粉离散元新模型的构建示意图Fig.21 Efficient discrete element modeling of multi-layer spreading

为了提高计算效率,我们构建了高效多层铺粉离散元计算新模型 (详细构建步骤见下述讨论).该新模型可以识别粉层中的未熔化球形金属颗粒,并直接使用可移动的球颗粒 (图21(f)),而不是采用大量的三角面片来表示未熔化球形金属颗粒,更加符合真实铺粉过程,增加了模拟的拟真程度.

高效多层铺粉离散元新模型的构建步骤如下.

(1) 在激光束选区熔化前,记录粉层中每一个金属颗粒在激光选区熔化细观尺度模型中的球心坐标和半径.

(2) 在激光束选区熔化前,记录激光选区熔化细观尺度模型中每一个单元的液体体积分数fb.这里的液体包括视为凝固液体的未熔化金属.液体体积分数表示单元中液体体积与单元体积的比值.

(3) 在激光束选区熔化后,记录激光选区熔化细观尺度模型中每一个单元的液体体积分数fa.

(4) 对于任意一个金属颗粒,根据金属颗粒的球心坐标、半径以及激光选区熔化细观尺度模型网格划分情况,计算激光选区熔化细观尺度模型中分布有该金属颗粒的单元,并记单元数量为N0.逐一对比这些单元的fb和fa.引入人为设定的参数α1和α2,其取值为较小的正数 (如0.005).当|1-fb/fa| ＞α1时,认为该单元内的金属颗粒熔化并流出该单元,或有熔融金属流入并黏附于该单元内的金属颗粒;当|1-fb/fa| ≤α1时,认为该单元内没有熔融金属流入流出.记录满足|1-fb/fa| ＞α1的单元数量为N.当N/N0＞α2时,认为该金属颗粒熔化 (包含部分熔化) 或被熔融金属黏附,在新高效多层铺粉离散元模型中需要使用三角面片复原其表面形貌.当N/N0≤α2时,认为该金属颗粒完全未熔化且未被熔融金属黏附,可以在新高效多层铺粉离散元模型中使用球颗粒表示该金属颗粒.

(5) 在激光选区熔化细观尺度模型中,令仅分布有未熔化且未被熔融金属黏附金属颗粒单元的液体体积分数fa=0,则可以得到如图21(d) 所示的仅含有熔迹、部分熔化金属颗粒以及被熔融金属黏附金属颗粒的激光选区熔化细观尺度模型.

(6) 通过stl 文件,将如图21(d) 所示激光选区熔化细观尺度模型的几何形貌信息输出至高效多层铺粉离散元新模型中,并通过金属颗粒的球心坐标和半径在新模型中生成代表未熔化且未被熔融金属黏附的金属球形颗粒 (图21(e),红色球形颗粒),则可以得到如图21(e) 所示的高效多层铺粉离散元计算新模型.

图21(e) 所示的高效多层铺粉离散元新模型包含球颗粒346 个,三角面片116 480 个,模型每秒可完成24～25 个时间步的计算.图21(b) 所示的多层铺粉离散元模型包含球颗粒0 个,三角面片275 904个,模型每秒可完成3～4 个时间步的计算.通过减少三角面片数量,多层铺粉离散元新模型将计算效率提升了6 倍以上.此外,从高效多层铺粉离散元新模型的构建步骤可以看出,选区熔化结束后,粉层中未熔化金属颗粒比例越大,减少的三角面片数量越多,新模型的计算效率提升越显著.

激光选区熔化工艺包括铺粉和选区熔化两个过程.除第一层铺粉外,其余铺粉过程都是多层铺粉.在成形区面积不变的情况下,随着铺粉层数的增加,非成形区未熔化金属颗粒的数量和表面积线性增加.因此,铺粉层数越大,未熔化金属颗粒比例越大,新模型的计算效率提升越显著.

高效多层铺粉离散元计算新模型使用球颗粒表示未熔化金属颗粒,消除了多层铺粉模型中,未熔化金属颗粒不能移动的非物理限制,提高了模型的计算效率.成形区表面形貌、非成形区颗粒与刮刀和粉堆颗粒的相互作用等是多层铺粉过程特有的因素,且对粉层质量有重要影响.高效多层铺粉离散元计算新模型在具有较高计算效率的同时,能够研究上述影响因素在多层铺粉过程中的作用机理,对深入研究激光选区熔化的铺粉机制具有重要意义.

6 结论

本文首先通过对比休止角和致密度的实验与模拟结果,获得了钨金属颗粒的滚动摩擦系数与表面能密度.然后,采用单层铺粉离散元模拟定量研究了粉层厚度和铺粉速率对粉层质量 (包括粉层形貌、粉层密度、配位数分布、Delaunay 四面体体积等效半径分布、球形度分布和粉层表面粗糙度) 的影响,确定了使粉层颗粒有效堆叠与紧密均匀排布的铺粉工艺窗口.最后,构建了高效多层铺粉离散元计算新模型,新模型使用可移动球颗粒表示完全未熔化金属颗粒,显著提高了多层铺粉离散元模拟的计算效率与拟真程度.