一直想找个机会把激光熔覆仿真从头到尾捋一遍,这次就拿一个比较典型的COMSOL激光熔覆案例来拆解。激光熔覆这个东西,说白了就是通过高能激光把金属粉末或者合金粉末熔化,在基体表面形成一层冶金结合的覆层,用来修复磨损件或者做表面强化。工艺参数稍微偏一点,就会出现气孔、裂纹、稀释率过高这些问题,全靠实验试错又太烧钱,所以这几年越来越多搞工艺的人开始用仿真先把窗口探出来。
COMSOL做这个事,优势在于多物理场耦合顺手,传热、相变、流动、应力变形都在一个环境里搭,不需要几个软件来回导数据。我这次做的案例是单道激光熔覆,基材是316L不锈钢,熔覆材料是Inconel 718粉末,激光功率1500W,扫描速度5mm/s,光斑直径2mm。下面把整个建模、计算和后处理过程拆开讲,每一部分我都会解释为什么这么设置,以及实测下来踩过哪些坑。
1. 项目背景与仿真目标
1.1 先把激光熔覆的物理过程讲清楚
激光熔覆的物理过程可以简化成三块。第一块是激光与材料交互,激光束照射到粉末和基体表面,一部分能量被反射,一部分被吸收转化为热量,导致局部温度瞬间升到熔点左右甚至更高。第二块是熔池行为,粉末熔化后形成熔池,熔池内部存在马兰戈尼对流、浮力对流,熔池表面还有表面张力作用。第三块是凝固与固体力学,激光移开后熔池快速冷却凝固,产生温度梯度和残余应力。
仿真想覆盖全部物理场不是不行,但代价是计算量巨大、参数不确定性也很高。所以大多数工程案例会在“精度够用”和“算得动”之间找一个平衡点。我这个案例的核心目标是预测熔池尺寸(宽度、深度)和基体热循环曲线,从而评估工艺窗口。至于熔池内部的对流细节,先用等效导热系数增强来近似,不直接解流场。
1.2 为什么选COMSOL而不是其他工具
市面上能做热-力耦合的软件不少,Abaqus可以算,ANSYS也可以算,但COMSOL有几个点比较适合激光熔覆这类移动热源问题。
首先是移动热源的实现方式灵活。COMSOL的“固体传热”接口里可以直接用函数表达式定义热流密度,并且支持用空间坐标与时间的依赖关系来实现热源的移动。你也可以把激光热源用“事件”接口做间歇式扫描,甚至配合“移动网格”处理复杂轨迹,但简单的直线扫描根本不需要那么复杂。
其次是多物理场耦合直接。传热模块算完温度场,可以直接把温度作为载荷施加到固体力学接口,不需要另外做单向耦合的数据传递。对于激光熔覆这种温度剧烈变化引发的应力变形问题,这种“温度场→结构场”的单向耦合路径是最常用、也是最稳的。
再者是对后处理比较友好。COMSOL里提取熔池轮廓、查看某一点的温度-时间曲线、切片看截面温度分布,这些操作都在同一个界面里完成,比导出到第三方软件再处理省事得多。
2. 几何模型与物理场设定
2.1 几何建模的简化与取舍
激光熔覆实际过程里,熔覆层是随着激光扫描逐渐成形的,表面不是平的,而是一条有一定高度和宽度的弧形隆起。如果完全按照真实形貌建模,需要用“动网格”或者“水平集方法”去追踪熔池自由表面,那是非常复杂的问题。
我在这个案例里采用了固定几何的简化处理:预先创建一条熔覆层几何,截面为弓形,宽度2.4mm,高度0.6mm,长度60mm,放置在基体表面上。基体是100mm×40mm×10mm的矩形板。有人说这样做会不会误差很大?我的理解是,如果你主要关注的是基体热循环、熔池宽度深度与工艺参数的关系,这个简化已经能给出很好的趋势。预成形的覆层几何在热传导意义上等效于粉末逐步熔化后的堆积体,只要热输入位置和覆层体积匹配,结果是比较可信的。
如果你确实要模拟粉末逐渐堆积的形貌演化,推荐用COMSOL的“变形几何”加“水平集”一起做。但这属于进阶玩法,计算稳定性对网格质量极其敏感,新手直接上手很容易发散。先把固定几何的单道熔覆跑通,再去碰自由表面,这是比较务实的路线。
2.2 物理场接口的选择与搭配
在这个案例里,我用了两个物理场接口:“固体传热”和“固体力学”。前者负责计算温度场,后者负责计算热应力与热变形。
“固体传热”接口里需要设置的项包括:固体域的热传导属性、初始温度、边界热通量,以及一个体积热源或者边界热源。激光热输入我选的是边界热源,即热流密度作用在基体上表面和覆层上表面。这里有一个容易犯迷糊的地方:激光能量是作用在粉末床上的,粉末对激光的吸收率远高于致密固体,所以热流密度表达式中要包含材料吸收率这个因子。
“固体力学”接口设置相对简单,主要是固定边界条件和温度载荷。温度载荷直接引用传热接口的温度场变量,然后用“热膨胀”子节点定义热应变。力学边界条件是基体底面固定约束。
有人会问,要不要加“层流”接口模拟熔池对流?我的建议是,如果你关心熔池内部流速和形貌的细微变化,可以加。但熔池对流的驱动机制很复杂,表面张力温度系数、黏度、热导率这些参数在高温区的数据都很不确定,算出来的流速量级参考意义大于定量意义。而且一旦耦合流场,网格量和求解难度会呈指数上升。对于工艺窗口筛选,不加流场完全够用。
2.3 边界条件与初始条件
初始温度设为20°C,环境温度也设为20°C。基体与空气接触的所有面都施加对流热通量,对流换热系数取10W/(m²·K)。辐射散热在高温段很显著,这个不能省,我在所有外表面加了表面对环境的辐射散热,表面发射率设为0.85。
这里需要注意的是,激光作用区域的边界热通量是用热流密度函数表达的,同时要把对流和辐射散热叠加进去。在COMSOL里实现方式是:在“边界热源”节点里定义激光热流密度,在“热通量”节点里定义对流和辐射散热。两个节点都作用在同一个边界上,COMSOL会自动累加,不需要人为去扣减。
另外要说一下热源作用的边界。我这里的做法是:激光直接作用在覆层上表面和基体未覆盖区域的上表面。模拟过程中,随着热源移动,边界上不同位置会被依次加热,这就实现了移动扫描的效果。
3. 材料属性与热源模型
3.1 材料参数的温度依赖性
激光熔覆的温度范围从室温到2000°C以上,材料属性如果全用常温常数,计算结果基本没法看。导热系数、比热容、热膨胀系数、杨氏模量、屈服强度这些都会随温度大幅变化。
316L不锈钢的导热系数在常温下约16W/(m·K),到1000°C左右会升到28-30W/(m·K)左右。Inconel 718的导热系数也是随温度升高而增大的。比热容同样有温度依赖性。在COMSOL里,材料属性可以定义为温度的分段函数或插值表。我个人的做法是:优先从材料库或者公开数据源找实验测试的插值数据,录入为“插值函数”,再赋给材料节点。
这里必须强调一个很容易被忽略的参数:相变潜热。材料从固态变成液态要吸收潜热,凝固时释放潜热。如果不考虑潜热,熔池尺寸会被明显高估或者低估,具体取决于你提取熔池的方法。COMSOL里处理潜热比较直接的办法是在“固体传热”接口的“热学属性”里勾选“相变材料”,设置相变温度区间和相变潜热。对于316L,固相线约1370°C,液相线约1450°C,潜热约260kJ/kg。Inconel 718的固液相线更高一些,约1260-1336°C,潜热约290kJ/kg。
3.2 移动热源的实现方式
移动热源是激光熔覆仿真的核心之一。最常用的热源模型是高斯热源,其热流密度分布为:
q(r) = (η × P × 2) / (π × R²) × exp(-2 × r² / R²)
其中,P是激光功率(W),η是材料吸收率,R是有效光斑半径(m),r是距光斑中心的距离(m)。这个公式的含义是:光斑中心的能量密度最高,随距离按高斯函数衰减。
在COMSOL中实现时,需要把r表达为坐标的函数。假设激光沿x轴正向移动,光斑中心位置为(x0 + v×t, y0),其中v为扫描速度,t为时间。那么边界上任意一点(x, y)到光斑中心的距离r为:
r = sqrt((x - x0 - v×t)² + (y - y0)²)
把这个表达式写入“边界热源”的热流密度表达式即可。实际书写时,x、y是COMSOL的内置坐标系变量,t是时间变量,v和P等可以在“参数”节点中定义。
我习惯用解析表达式直接在边界热源里写一个函数q_laser(x, y, t)。既要保证表达式简洁,又要方便批量修改参数。如果之后要扫描不同功率和速度,只需要修改参数值即可,不用动模型结构。
提示:使用高斯热源时,光斑半径R的定义要统一。有的文献用1/e²半径,有的用1/e半径,公式表达形式对应不同,混用会让热流密度差出一大截。
4. 网格剖分与求解控制
4.1 网格策略与分区
激光熔覆的温度场特点是移动热源附近温度梯度极为陡峭,远离热源的地方梯度平缓。如果全局都用细网格,计算量会爆炸;全局用粗网格,熔池附近的温度峰值和尺寸又会被严重低估。
我的做法是分区网格:覆层区域和基体上表面附近1mm深度范围内的区域采用细密网格,单元尺寸控制在0.2mm左右;基体其余部分用过渡网格,逐渐放大到2-4mm。这里用到了COMSOL的“自由四面体网格”和“扫掠网格”的组合。
覆层因为是弓形截面,尺寸较小,我用自由四面体网格,最大单元尺寸0.15mm。基体上表面下方做一个边界层网格,彻底加密热流注入侧的散热路径。边界层设置3层,第一层厚度0.05mm,增长因子1.2。基体下部区域用自由四面体,最大单元尺寸3mm。
为什么网格尺寸这么敏感?因为高斯热源的能量密度是靠指数函数衰减的,如果单元太大,网格节点上的r值会“跳”掉,热流密度随空间的变化就体现不出来。我做过对比:覆层处网格从0.4mm细化到0.15mm,熔池宽度计算值变化了约15%,而0.15mm到0.1mm的变化不超过3%。这说明0.15mm已经接近收敛。做仿真一定要先做网格无关性验证,这一步偷懒后面全是白算。
4.2 瞬态求解器设置
这个案例是瞬态问题,因为激光在移动,温度场随时间变化。激光扫描速度5mm/s,覆层长度60mm,所以总扫描时间约12s。扫描结束后还需要继续计算一段时间让基体冷却,我设定总计算时间为15s。
时间步长的选取要谨慎。COMSOL的瞬态求解器默认用自动时间步长,但你最好给一个约束。我通常设置初始步长0.001s,最大步长0.05s。因为高斯热源的移动速度是5mm/s,相当于每个时间步内光斑移动0.005-0.25mm,要保证光斑每步移动不超过一个细网格单元尺寸,否则热源会在空间上出现“跳跃”,导致温度振荡。
求解器我直接用默认的“直接”求解器搭配PARDISO,这个问题规模不大,自由度大约几十万,直接求解器反而比迭代法更稳。如果网格量超过200万自由度,再考虑用迭代法。PARDISO对这类热-力耦合问题鲁棒性不错,实测没有遇到内存溢出的问题。
5. 结果分析与后处理
5.1 温度场演化与热循环曲线
跑完计算,第一步先看整体温度分布动画。我习惯在结果节点里新建一个“二维绘图组”,切一个垂直于扫描方向的截面,查看温度分布云图。激光经过时,熔覆层顶部温度瞬间冲高,基体表面以下形成明显的温度梯度。
要判断工艺是否合理,可以提取基体表面以下0.5mm深度处一点的热循环曲线。具体操作是:在“派生值”里创建“体平均”或者“探针”,选择目标点坐标,然后绘制该点的温度-时间曲线。一个典型的激光熔覆热循环曲线会有一个很尖锐的升温峰,峰值温度往往超过1400°C,然后快速跌落,直到下一次激光经过再次升温。
这个热循环曲线可以用来评估基体是否发生了不利的相变。比如马氏体不锈钢对冷却速度很敏感,如果峰值温度过高、冷却速度过快,可能引发淬硬裂纹。这里就可以通过调整激光功率、扫描速度来实现热循环控制。
5.2 熔池尺寸的定量提取
熔池宽度的提取方法不唯一。我采用的方法是:在结果中用“等值面”提取温度等于材料液相线的面。由于温度场是瞬态的,熔池尺寸随时间变化,所以要选取一个特征时刻,通常是激光经过某个位置并达到准稳态时。
准稳态的判断方法是:在扫描方向中部选取一个固定点,观察该点的温度-时间曲线。如果该点的峰值温度在激光经过时达到稳定值,且前后几道扫描的峰值差异小于5%,就可以认为达到了准稳态。采用准稳态时刻的等温线提取熔池宽度和深度。
需要注意的一个细节是:如果你用的材料是纯物质,相变温度是一个点,但合金材料有固相线和液相线,之间是一个糊状区。提取熔池边界时,建议采用液相线等温线而不是固相线,因为液相线能更准确反映实际熔池边界。我在这个案例中,液相线等温线得到的最大熔池半宽约1.15mm,熔池深度约0.9mm。这与实验金相测量值吻合度在10%以内,说明模型参数取对了。
5.3 残余应力与变形评估
温度场算完后,固体力学接口会在同一套网格上计算热应力。这里有必要说一下,热应力的绝对值对材料高温力学性能参数极为敏感,尤其是高温屈服强度。如果手头没有实测的高温力学性能数据,直接用室温强度参数来计算,残余应力会明显偏高,但分布趋势是有参考价值的。
我提取的结果包括扫描结束后的残余应力分布和基体变形。典型的情况是:熔覆层横向残余应力表现为拉应力,峰值出现在熔覆层中心与基体结合界面附近;基体底部则表现为压应力。变形方面,基体会有轻微的翘曲,翘曲量以底面的竖向位移来量化。我这里算出的最大翘曲位移约0.12mm,随着基体厚度增加,翘曲量明显减小。
这个结果对生产有直接意义:如果你在做薄壁件修复,基体厚度只有2-3mm时,激光功率稍微高一点,工件就翘曲得没法用;但基体厚度超过10mm,翘曲就不算主要问题,重点转向熔覆层裂纹控制。
6. 常见问题与排查实录
6.1 温度场发散
这是我见过最多的问题,也是初学者最容易卡住的环节。温度场发散的表现是:计算几步后温度突然变成几万甚至几十万度,或者出现NaN。
排查顺序按概率从高到低:
- 单位制错误。检查激光功率是否用了W,尺寸是否用mm,热导率是否用了W/(m·K)。COMSOL的单位系统虽然会自动换算,但自定义表达式里的常数很容易出错。
- 边界热源表达式书写问题。仔细检查距离r的计算表达式中是否用了正确的坐标变量,是否有除以零的情况。
- 网格太粗。热源半径如果只覆盖一两个单元,能量集中导致单元内部温度梯度超过数值极限,就会发散。强制把光斑覆盖区域的单元尺寸降到光斑半径的1/5以下。
- 初始时间步长太大。瞬态求解的第一步如果跳得太远,非线性迭代会失败。把初始步长减小一个数量级,往往立刻就好了。
我之前遇到一个特别隐蔽的问题:边界热源表达式里用了平方根函数sqrt,而sqrt的参数在极少数时间步下因为数值误差变成负的极小的值,导致返回NaN,整个计算直接崩掉。解决方案是把sqrt参数写成max(0, 表达式),从根上规避这个问题。
6.2 计算时间过长
有人反映说网格和求解器设置和案例一致,但算一个简单的单道熔覆要十几个小时,这就很不正常了。
第一,看看自己是不是全局用了细网格。很多人习惯“为了保险起见”把网格全部设为0.2mm,但这个模型全域用0.2mm网格,单元数量会从几十万变成几百万,慢是必然的。老老实实分区,把细网格限制在热影响区。
第二,检查时间步长。如果你手动设置了固定时间步长0.001s,全程15s就要算15000步,而每步都在求解几十万个自由度,当然慢。让自动时间步长接管,只在关键阶段约束步长,速度能提升一个量级。
第三,合理使用对称性。如果扫描轨迹是直线且几何对称,完全可以取一半模型,在对称面上施加绝热边界条件。这样自由度减半,计算时间差不多减半。这里提醒一句话:加上对称边界条件之后,激光热源的表达式中y0也要变化,别忽略了镜像侧的热源坐标。
6.3 熔池尺寸与实验对不上
这是所有仿真的终极考验。模型再精巧,如果结果和实验差得远,一切都是白搭。
优先级最高的排查项是材料吸收率。不同材料在不同温度、不同表面状态下对1064nm光纤激光的吸收率差异极大。固件金属常温吸收率只有30%-40%,粉末状材料因为多次反射吸收率可能达到60%-70%。我见过太多案例是功率乘以了一个想当然的吸收率0.5,然后熔池尺寸就差出一倍。解决思路是:先用一个小的试算,对比实验测得的熔池宽度,反推该工艺条件下的等效吸收率,再用这个标定后的值做后续工艺窗口扫描。
其次是热源模型的选择。高斯面热源适合导热为主的熔覆模式,但如果实际工艺中激光与粉末作用形成了深熔池,需要考虑体热源模型。一般熔深小于光斑半径时,面热源就可以;熔深显著大于光斑半径时,就得用锥形体热源或者双椭球体热源。
最后是覆层几何的激光路径。我在前面提到预成形覆层几何是一个简化的做法。如果在你的实际案例里,粉末是同步送粉状态,激光同时作用在粉束、已熔覆层和基体上,那么边界热源的施加位置就与预成形几何的拓扑空间存在差异。这个问题在低扫描速度、大熔覆量时会变得比较显著。
6.4 残余应力计算不收敛
固体力学模块出现不收敛的几率比传热低,但一旦出现就很头疼。最常见的原因是高温下材料弹塑性参数不足。316L在温度超过1200°C后,屈服强度急剧下降,如果材料库里的高温数据只有室温到800°C,计算中温度超过数据范围之后COMSOL会强制外推,容易造成应力振荡。
我的处理方式是:在材料属性中将高温段参数设置为一个很小的屈服强度及相应的硬化模量,避免数值刚度过大。虽然这个做法会牺牲一些高温区的精度,但换来的是计算稳定性和整体的残余应力分布合理性,从工程应用角度是值得的。
另一个问题是热膨胀系数的定义。COMSOL里热应变的计算依赖“参考温度”和“热膨胀系数”。如果我默认参考温度是20°C,而初始温度也是20°C,那么没有问题。如果参考温度和初始温度不匹配,就会出现初始时刻就有巨大热应力的怪异结果。这个检查通常在几分钟内就能排除。
7. 个人实操总结
这个案例做完,我最大的感受是:COMSOL激光熔覆仿真难的不是软件操作,而是对工艺物理过程的理解和参数标定的耐心。任何一次仿真与实验对不上,大概率不是求解器出了问题,而是你对边界条件的抽象偏离了真实物理过程。
我自己比较有效的习惯是:仿真金相验证交替进行。第一批仿真用名义参数跑一遍,得到熔池宽度深度预测,然后打一个单道试件,量出实际熔池尺寸,反向修正热源吸收率和光斑半径。基本两三轮下来,这个模型的预测能力就稳定了。有了这套标定好的模型,后面不管怎么改扫描速度、激光功率、送粉量,都能先算一轮,把实验次数压缩到原来的五分之一。
最后再分享一个小技巧:批量扫描工艺窗口时,用COMSOL的“参数化扫描”功能,把功率设成一组值[1200, 1500, 1800, 2100],速度设成一组值[3, 5, 8, 12]。一次把所有组合跑完,后处理时直接输出熔池宽度和深度对应表。这个操作能帮你建立一个工艺参数-熔池特征的数据表,拿到现场进行工艺调试,比凭感觉试参数稳得多。