简介面向增材制造与固体力学方向研究者及工程技术人员的 COMSOL 6.1 激光选区融化SLM仿真资料围绕单道与多道成形中的温度场、应力场展开覆盖热固耦合、活化功能及圆柱形增材等建模要点适合希望复现 SLM 热力演化过程、评估工艺参数影响并拓展到其他增材制造场景的读者。压缩包共 3 个文件含 1 份 PDF 说明、1 份 Markdown 实践文档和 1 份 HTML 页面整体约 480KB分别承担模型讲解、操作要点整理与在线浏览用途体积轻便却信息集中。已有 276 人学习说明其在相关主题中具备一定参考价值。资料重点呈现模型注释、可修改性与收敛性处理读者可据此理解温度场与应力场的计算流程借鉴参数设置、边界条件与求解策略并在此基础上开展单道多道对比分析为增材制造热固耦合仿真提供可复用的建模思路。1. 熔池只有几十微米应力却能把零件掰弯SLM 仿真里最反直觉的一件事激光选区熔化SLM的熔池直径通常只有几十微米激光扫过一点的时间以微秒计。按直觉这么局部的热扰动应该瞬间被周围金属吸走可实际打印出来的零件常在取件那一刻肉眼可见地翘曲甚至直接开裂。温度场究竟怎么把应力场塑造到这种程度是增材制造数值仿真里最值得深入分析的一环。热固耦合的难点不在方程本身——传热方程和弹性力学方程课本里都有——而在尺度跨度熔池几十微米零件几十毫米中间还夹着固液相变、材料参数随温度剧烈变化、成千上万次热循环。COMSOL 6.1 这类多物理场平台的价值是把这些耦合关系显式摆到同一张网格上让每一次热循环留下的塑性应变都能被累加。下面按“温度场怎么算准 → 应力场怎么接上 → 铺粉怎么在模型里长出来 → 哪里最容易翻车”的顺序拆。2. 移动热源与相变潜热温度场为什么难算2.1 高斯热源在网格上“走”一步物理量发生了什么热源模型的选择直接决定了后面所有结果的可信度。SLM 里最常见的是面内高斯分布、深度方向按指数衰减的体热源它的物理含义是把激光能量在光斑半径内按高斯权重分配同时在铺粉层厚度内被金属粉末和已凝固层吸收。写成表达式就是功率除以光斑面积再乘一个归一化因子 2/π这样在无穷大平面上积分正好等于总功率乘以吸收率。吸收率这个参数在工程上通常取 0.3 到 0.5粉末床因为多次反射会比实体略高但这个值的不确定性本身就能让熔池温度摆动上百开尔文。真正的问题出在离散化上。热源在解析解里是连续移动的在有限元网格上却是“跳”着走的每个时间步高斯中心平移一个 v·dt 的距离。如果这个距离大于热源附近单元的特征尺寸高斯峰就会从一个单元直接跳到下一个单元中间那个单元永远等不到峰值温度。温度曲线上表现为明显的锯齿熔池边界忽大忽小后处理出来的熔池长度完全是随机的。常见的实现思路是给热源区域设一个硬约束v·dt ≤ h/2h 是激光路径上最小的单元尺寸。举一组实际会遇到的数光斑特征半径 r0 40 μm热源附近单元尺寸 h 15 μm扫描速度 v 800 mm/s那么时间步上限约 1.9×10⁻⁵ s。一条 10 mm 长的扫描线要走 12.5 ms也就是六百多个时间步一个 3 mm×3 mm×2 mm 的小试样几百层几百道步数轻易上百万。这就是为什么多数 SLM 仿真最后都在“减少层数、缩短扫描线”上妥协。提示网格在热源附近至少要能分辨 r0 的三分之一深度方向铺粉层厚度内至少放三层单元否则潜热释放和温度梯度都会失真。下面这段是 LiveLink for MATLAB 里定义移动体热源的大致写法方法名随版本略有差异% 移动高斯体热源面内高斯深度方向指数衰减 model.variable(var1).set(P_laser, 200[W]); % 激光功率 model.variable(var1).set(eta, 0.35); % 粉末床有效吸收率 model.variable(var1).set(r0, 40[um]); % 光斑特征半径 model.variable(var1).set(v_scan, 800[mm/s]); % 扫描速度 model.variable(var1).set(d_opt, 35[um]); % 光学穿透深度 model.variable(var1).set(t0, 0[s]); % 起始时刻 % 热源中心随时间的 x 坐标y 为道间距方向 model.variable(var1).set(x_src, 0[mm] v_scan*(t-t0)); model.variable(var1).set(q_laser, ... 2*P_laser*eta/(pi*r0^2)*exp(-2*((x-x_src)^2y^2)/r0^2)*exp(-z/d_opt));q_laser的量纲是 W/m³直接作为广义源项塞进传热方程的右端。d_opt不是随便给的数它决定了能量在深度方向分给几层单元取值和铺粉层厚同量级时热源才真正是“体”的而不是“面”的。把d_opt设得过小模型退化成面热源熔池深度会被系统性低估设得过大激光会“穿透”到基板上熔池深宽比完全失控。2.2 潜热为什么必须用表观热容而不是直接给焓熔化潜热的量级经常被低估。以常见的不锈钢粉末为例比热容约 700 J/(kg·K)熔点附近潜热约 270 kJ/kg。如果把这 270 kJ/kg 摊到一个只有 20 K 宽的相变区间里等效比热容要额外增加 13500 J/(kg·K)——是固相比热的近二十倍。忽略潜热的效果不是“稍微热一点”而是熔池温度虚高、凝固阶段冷却速率被高估、温度梯度被放大最终算出来的残余应力可能比实际高出一大截。处理相变有两条路。焓法把内能写成温度的函数 H(T)导热方程改写成 ∂H/∂t 的形式物理上最干净但 COMSOL 里要自己构造 H(T) 并处理它的反函数实操上麻烦。更常见的是表观热容法在相变区间内给比热容叠加上一个与液相分数导数成正比的项数学上和焓法等价但形式上仍然是一个 Cp(T) 的表达式直接填进材料节点就行。% 液相分数用平滑阶跃避免导数在 Ts / Tl 处不连续 model.variable(var1).set(dT_ph, 20[K]); % 相变区间展宽 model.variable(var1).set(phi_l, ... flc2hs((T-(Ts-dT_ph/2))/dT_ph, dT_ph)); % 0 - 1 平滑过渡 model.variable(var1).set(Cp_eff, ... Cp_solid(T) Lf*d(phi_l,T)); % 表观热容flc2hs的第二个参数是过渡宽度它的作用是把这个本来是阶跃的函数磨圆让d(phi_l,T)有界。这里有个纯粹的数值取舍dT_ph取太小比如 1 K表观热容的峰值会高到让雅可比矩阵条件数爆炸牛顿迭代开始反复被截断取太大比如 100 K相变被抹平到接近熔点以下的温度区间熔池边界位置提前熔池尺寸偏大。实践里 20 到 50 K 是个比较稳的窗口具体值还要和材料熔点、时间步一起调。Ts和Tl也不该取手册上的平衡值。SLM 的冷却速率在 10⁵ 到 10⁶ K/s 量级凝固过程中的固液相线间隔会被显著拉宽有些合金的实际凝固区间比平衡相图宽出几十开尔文。用平衡值算出来的熔池凝固时间偏短热影响区的温度历史会失真。2.3 材料参数随温度跳变带来的非线性热导率、比热容、密度、热膨胀系数、屈服强度这些参数在熔点附近全部是剧烈变化的。热导率从固相的十几 W/(m·K) 到液相略有上升但液相里真正的传热机制是对流——马兰戈尼效应把熔池里的热量从中心往边缘搬。不显式求解流场的话工程上惯常的做法是把液相热导率人为放大 2 到 3 倍用一个等效导热系数把对流的效果折进去。这个处理当然不严格但它把熔池形状拉回到和实验接近的量级比完全不管要靠谱得多。真正的麻烦是数据缺失。商业材料库给出的热物性通常只到 1000 K 到 1500 K而熔池峰值温度可能到 2500 K 以上。外推的方式有好几种——线性延长、取常数值、或者用液相的值硬顶——不同选择对熔池过热度的影响是几百开尔文的差异。我的习惯是固相区严格用实测数据做分段线性插值熔点以上取常数并明确记录这个假设因为那个区间本来就没什么可信数据与其编一条曲线不如承认它的不确定性。参数跳变还会直接打击求解器。热导率在固液相之间可能有两倍以上的差别如果材料节点里直接写成阶跃函数每次牛顿迭代跨过相变温度都会让残差剧烈变化。用和潜热同样的思路所有在熔点附近有跳变的参数都该用flc2hs平滑一下过渡宽度和dT_ph对齐。2.4 时间离散隐式步长为什么不能贪大传热方程是抛物型的隐式格式本身是无条件稳定的所以直觉上时间步应该可以随便取。问题出在“稳定”和“准确”是两回事。移动热源经过某个位置的时间尺度是 r0/v对上面那组参数就是 5×10⁻⁵ s 量级如果时间步取到 10⁻³ s热源一步跨过二十个光斑那个位置根本感受不到峰值温度。温度历史被严重平滑峰值温度偏低冷却速率偏离真实值。COMSOL 里的瞬态求解器默认用自由步长加 BDF 时间积分内部会根据局部误差估计自动调整。这个机制在纯传热问题里工作得不错但一旦耦合了相变和随温度变化的材料参数误差估计会变得非常保守步长被压到极小计算卡死。更可控的做法是给求解器设一个输出时间列表同时限制内部最大步长让求解器在“必须输出”和“可以大跨步”之间有个明确的边界。还有一个容易被忽略的事件激光开启和关闭的瞬间。这两个时刻热源功率是阶跃变化的求解器如果恰好跨过这个时刻会把它当成一个光滑变化处理结果就是开头几步的温度完全不对。手动在这些时刻插入输出点或者用事件接口显式标记都能避免这个问题。3. 热-固耦合为什么顺序耦合成了工程默认解3.1 全耦合与顺序耦合的自由度规模差从方程上看热固耦合是双向的温度场通过热应变影响位移场位移场通过变形做功反过来影响温度场。COMSOL 里确实可以把传热和固体力学放进同一个全耦合求解器一次性解出所有自由度。问题在于规模。三维传热每个节点一个自由度三维固体力学每个节点三个。十万节点的模型全耦合系统的规模是四十万而且雅可比矩阵是非对称的不能再用对称求解器那套省内存的技巧。更致命的是条件数传热矩阵的量级和弹性矩阵的量级差了好几个数量级放进同一个矩阵里预条件子很难做好。顺序耦合的假设是力学变形对传热的影响可以忽略。这个假设在 SLM 里几乎总是成立——熔池的流动和传热由激光能量主导固态区的变形量在微米级对热传导路径的影响小到测不出来唯一的例外是零件和基板之间出现间隙或接触状态变化。所以在工程实践里先解一个纯瞬态传热把每个时间步的温度场存下来再把温度作为体载荷插进力学场是最常见的路线。代价是插值误差。两套网格如果不一致温度要从传热网格映射到力学网格上映射过程必须用保守型插值否则总能量不守恒局部温度会凭空多出来或少掉。COMSOL 的耦合算子默认会处理这件事但如果两套网格疏密差别很大比如传热网格在热源附近细到 10 μm、力学网格为了省算力粗到 100 μm映射出来的温度梯度会被严重平滑热应力随之被低估。3.2 热应变怎么进入力学方程小变形框架下总应变分解成弹性、塑性、热应变三部分热应变张量是各向同性的等于热膨胀系数乘以温度相对参考温度的差值。这里面有两个坑。第一个是热膨胀系数的定义。材料手册给的 α 通常是相对于室温的割线值也就是 (L(T)-L(293K))/(L(293K)·(T-293K))而本构方程里需要的是切线值 dL/dT。这两者在温度跨度大的时候差别不小。COMSOL 的热膨胀节点默认按割线处理如果你的材料数据是切线值直接填进去会系统性偏差。第二个是参考温度Tref。它不是一个随便填的数物理上对应“零件处于无应力状态”的那个温度。如果仿真从室温粉末床开始Tref就该是室温如果工艺里有基板预热Tref应该取预热温度否则整个零件从第一步就带着一个均匀的热应变残余应力里混进一个本不该有的偏置。% 顺序耦合传热先算完力学场读入温度历史 model.physics(solid).feature(lemm1).feature(the1).set(Tref, 293.15[K]); model.physics(solid).feature(lemm1).feature(the1).set(alpha, alpha_solid(T)); % 温度来自传热场的插值而不是重新解一次 model.physics(solid).feature(lemm1).feature(the1).set(T, T_ht);T_ht是传热场映射过来的温度变量力学场只把它当成一个已知的载荷不参与求解。塑性部分用温度相关的屈服面和强化曲线SLM 的循环加热冷却特征意味着等向强化不够用——材料在反复的拉压交替下会出现包辛格效应随动强化模型比如 Chaboche 的多背应力形式更贴近实际但参数标定成本高多数工程仿真还是用等向强化加一个折减的屈服强度来近似。3.3 约束条件决定了残余应力的空间分布残余应力不是材料属性是“约束”的产物。同样一条扫描路径基板夹得死和夹得松算出来的应力云图能差出一倍。最常见的错误是把基板底面全固定。这个边界在数学上没有问题物理上却是错的真实工艺里基板靠几个螺栓压在平台上允许一定的转动和面内位移。全固定相当于给零件加了一个无限刚性的约束热膨胀被完全挡住应力被人为放大。合理的做法是只约束刚体运动的六个模态——比如底面三个点分别约束法向、两个切向、再加两个点约束转动——剩下的让它自由变形。基板本身的厚度也影响巨大。薄基板在热循环下会跟着翘吸收掉一部分热应变厚基板相当于刚性平台热量往下传得快冷却速率高应力反而更大。仿真里如果把基板简化成一层薄壳或者干脆省略结果基本不可用。取件之后应力释放、零件翘曲是另一个必须单独算的工况。做法是在第一个 study 算完之后把基板区域的单元钝化掉再跑一个静力平衡看零件在新的约束状态下的变形。这一步经常被跳过但它恰恰是现场最关心的那个量。4. 生死单元与多层多道铺粉过程怎么在模型里“长”出来4.1 单元激活的数值本质SLM 是一层一层长出来的仿真里要让材料“凭空出现”。COMSOL 的增材制造模块提供了单元激活功能底层机制其实很朴素给每个单元定义一个激活指示变量未激活时为 0激活后为 1然后把材料属性——热导率、比热容、密度、弹性模量——全部乘上这个指示变量。这里有个绝对不能踩的坑比热容和密度不能真的乘 0。传热方程里的瞬态项是 ρ·Cp·∂T/∂t如果这两个都变成 0那一块区域的矩阵对角线就是零整体矩阵奇异求解器直接报错。惯常做法是乘一个极小值比如 10⁻⁶让矩阵保持非奇异同时它对结果的影响小到可以忽略。第二个细节是激活瞬间的温度初始化。新激活的单元如果初始温度取 0 K 或者绝对零度附近的默认值激活那一刻会出现一个巨大的温度跳变热流密度瞬间冲到非物理的量级时间步被砍到 10⁻¹² 秒计算直接卡死。正确做法是让新单元的初始温度继承下方已凝固层的温度或者干脆用一个平滑的斜坡在几百微秒内把材料属性从“几乎不存在”过渡到“完全存在”。第三个细节是激活顺序和网格的匹配。如果一层单元在同一时刻全部激活相当于整层粉末瞬间变成金属并瞬间被加热熔池会变成一个扁平的方块。真实过程是激光逐点扫描所以激活必须跟着热源位置走单元激活的判断条件通常写成热源中心到单元中心的距离小于某个阈值或者按扫描路径参数化地给出激活时刻表。4.2 扫描策略对残余应力的影响扫描策略是 SLM 里少数几个“不改材料、不改设备、只改路径就能显著改善残余应力”的杠杆。常见的有长边扫描、短边扫描、棋盘格分区、层间旋转 67°。长边扫描的扫描线长温度梯度沿着线的方向拉得很开冷却时收缩量大端部的应力集中明显。短边扫描把长线切成短段单段的热输入小温度梯度缓和但道数变多道间的重熔区增加。棋盘格把一层切成若干小岛每个岛单独扫描岛与岛之间留出搭接区好处是每一小块的热积累都很有限坏处是岛边界处容易出现未熔合和孔洞。层间旋转 67° 是工业上最常见的做法它的作用不是让应力均匀而是打断热流路径的周期性。如果每层都沿同一方向扫热量会沿着扫描线方向一直往前传形成一条贯穿多层的“热通道”冷却后这条通道上下的收缩量和其他区域不一致容易诱发起裂。旋转一个非特殊角度67°、45° 都是常见值前者因为不与层内道间距形成整数倍关系而更受欢迎能让热流方向在层间被打散。仿真里复现这些策略的难点在于路径的表达。把扫描路径写成参数化曲线用弧长参数驱动热源位置是最干净的做法如果路径是离散的点序列就要用插值函数或者分段线性函数把时间映射到坐标上段与段之间的衔接处容易出现速度不连续导致热源在拐角处“停顿”局部过热。这是很多仿真结果里熔池在扫描线端点异常膨胀的原因。4.3 时间步、层间冷却与算力的三方拉锯一层扫描的时间和层间冷却的时间通常是同一个量级甚至后者更长。真实工艺里铺粉加刮平要几秒到十几秒而扫描一层只要几十毫秒。如果老老实实把冷却时间也算进去计算量的绝大部分都花在“温度几乎不变”的阶段性价比极低。常见的加速手段有三种。第一种是把层间冷却时间压缩比如真实 10 s 压到 0.1 s代价是热积累被低估层间温度偏低残余应力偏小。第二种是用“跳层”策略算几层跳几层跳过的层用插值补上这对温度场的宏观分布影响不大但对熔池尺度的细节完全失真。第三种是只算一个代表性体积单元把多道多层的边界条件用对称或者周期性条件近似这种做法能保住熔池细节但没法回答零件级的变形问题。这三种手段本质上都是在“熔池尺度”和“零件尺度”之间做选择没有哪个是普适的。选哪个取决于你要回答什么问题关心熔池形貌和孔隙缺陷就得保熔池、牺牲层数关心整体翘曲和残余应力分布就得保零件、牺牲熔池细节。这个取舍在第 6 章还会再展开。5. 排查清单SLM 热固耦合最容易踩的五个坑5.1 坑一熔池凭空消失温度曲线爬不到熔点现象很典型求解顺利完成后处理一看最高温度只有 1200 K连熔点都没到熔池根本不存在。原因通常有两个。一是热源表达式里的单位错位比如光斑半径填了40[um]但公式里忘了平方或者功率填成200而没带[W]COMSOL 会按 SI 基本单位解释成 200 W看起来没错但如果几何单位是 mm 而表达式里混用了就会出现量级偏差。二是吸收率和穿透深度共同作用把能量摊得太薄——d_opt设成几百微米时同样功率分布在几十倍的体积里峰值温度自然上不去。解决办法是先做一个量级校验用解析的 Rosenthal 解或者简单的一维估算算一下给定功率和扫描速度下熔池大致能到多少度和仿真结果对一遍。差一个数量级一定是单位或者公式的问题差百分之二三十才是模型精度的问题。5.2 坑二温度曲线锯齿振荡时间步被砍到极小现象是某个节点的温度历史不是光滑的峰而是一串尖刺求解日志里时间步长一路被压缩到 10⁻¹⁰ 秒以下进度条几乎不动。原因多半是时间步和网格不匹配热源一步跨过了多个单元高斯峰在网格上“跳跃”每个单元交替经历“被峰扫到”和“没被扫到”。另一种可能是材料参数里存在阶跃比如热导率在熔点处直接从一个值跳到另一个值牛顿迭代每次跨过这个温度都会剧烈震荡。解决办法分两步先把最大时间步限制在h/(2v)以内确认锯齿消失如果还有残留把所有在熔点附近跳变的材料参数用flc2hs平滑过渡宽度不小于相变区间的展宽。这两步做完九成以上的锯齿问题都能解决。5.3 坑三迭代不收敛报“非线性求解器不收敛”现象是求解器在某个时刻反复回退步长最后还是失败日志里写的是牛顿迭代达到最大次数。原因往往出在激活瞬间或者相变区间。新单元激活时材料属性从 10⁻⁶ 跳到 1如果激活的单元数量在同一时间步里很多刚度矩阵和传导矩阵会同时发生大变化牛顿迭代的初值离解太远。另一个常见原因是塑性迭代和传热迭代在同一个耦合步里互相打架尤其是热膨胀系数取切线值的时候。解决办法是给激活加时间斜坡把属性变化摊到几十个时间步上给非线性求解器加阻尼把更新量的步长限制在合理范围内实在不行就把传热和力学彻底分开成两个 study中间只传温度场不让它们在同一个牛顿循环里互相牵制。5.4 坑四位移全为零或者零件直接“飞走”位移全零通常意味着力学场的约束给多了——所有边界都被固定热应变无处释放位移自然是零但应力大得离谱。零件“飞走”则相反约束不足刚体模态没被消掉求解器解出一个巨大的刚体位移。这两种现象的根因是同一个约束条件的物理意义没想清楚。合理的做法是先问“真实工艺里这块材料被什么挡住了”只保留那些真实存在的约束剩下的用弱弹簧或者刚体运动抑制来消掉数值上的奇异性。刚体运动抑制不会引入额外的应力这是它比全固定优越的地方。5.5 坑五结果看起来很合理数值却差一个量级最危险的一种情况云图颜色漂亮熔池形状也对但残余应力的数值和实测差了一倍。这种“看起来对”的错误最难查。常见根源有三个吸收率取值偏离实际粉末床的真实吸收率随粉末粒径、表面氧化状态变化很大0.3 和 0.5 能让应力差 40%参考温度设错导致整体带一个偏置应力塑性模型里用的屈服强度是室温值而熔池附近温度高达上千度实际屈服强度只有室温的十分之一。核查办法是固定看三个量熔池的最高温度和尺寸、冷却速率对凝固组织有直接影响、以及取件后零件的最大翘曲量。这三个量和实验或文献一对比量级对不对立刻现形。只看应力云图的颜色分布是最容易自欺欺人的做法。6. 设计取舍二维切片、单道外推还是硬算整个零件仿真做久了会发现SLM 建模的核心矛盾不是“会不会算”而是“算哪个尺度”。熔池尺度需要微米级网格和微秒级时间步零件尺度需要毫米级网格和秒级时间步两者相差三到四个数量级任何试图同时覆盖的模型都会被计算量压垮。一个实用的判断方式是先明确模型要回答的问题属于哪一类。关心熔池形貌、球化、匙孔、未熔合就做单道或者几道的二维/三维小模型网格细到 5 μm时间步细到微秒层数不重要。关心残余应力分布和零件翘曲就用“等效热源”或者“逐层激活 粗网格”的方式把熔池细节抹平用体积热源或者分层的温度边界条件代替重点放在约束条件和材料本构上。关心扫描策略的优劣对比用一个中等规模的板状试样十几层、几十道网格 20 μm 量级能分辨道间搭接但分辨不了熔池内部。模型层级典型网格尺度时间步能回答的问题算不出来的东西单道熔池3–10 μm1–10 μs熔池形貌、深宽比、冷却速率零件变形、层间热积累多道多层小试样15–30 μm10–50 μs道间搭接、层间热积累、局部应力整体翘曲、基板效应零件级等效100–500 μm0.1–1 ms整体温度分布、残余应力趋势、翘曲熔池细节、孔隙缺陷这张表最有用的地方不是选哪一行而是提醒你任何一行都在放弃某些信息而这些被放弃的信息恰恰可能在真实工艺里是关键。二维切片能省掉一个数量级的算力但它隐含了“面外方向均匀”的假设对于扫描线端部、零件边缘、薄壁结构这个假设直接失效。单道外推到多道忽略了道间的重熔和热积累应力会被系统性低估。我自己的习惯是任何模型开算之前先用解析估算把熔池温度和尺寸对一遍再用一个极简的二维模型把热源参数标定好最后才上三维。这个流程看起来多花时间但它把“参数错了”和“模型选错了”这两类问题分开了——前者在二维阶段就能发现后者要到三维才会暴露早发现早止损。那次教训来自一个薄壁件为了省时间直接用了二维切片扫出来应力很小结论是工艺可行。实际打印出来壁厚方向翘得没法用。后来才想明白二维切片把面外方向的约束当成了无限刚体热应变全被压住了应力当然小。从那以后我每次建模前都强制走一遍“这个假设在哪个方向失效”的自检不再让网格疏密决定结论。希望这些经验帮到你少走一点弯路。本文还有配套的精品资源点击获取 SLM热固耦合仿真：移动热源、相变潜热与残余应力控制
简介面向增材制造与固体力学方向研究者及工程技术人员的 COMSOL 6.1 激光选区融化SLM仿真资料围绕单道与多道成形中的温度场、应力场展开覆盖热固耦合、活化功能及圆柱形增材等建模要点适合希望复现 SLM 热力演化过程、评估工艺参数影响并拓展到其他增材制造场景的读者。压缩包共 3 个文件含 1 份 PDF 说明、1 份 Markdown 实践文档和 1 份 HTML 页面整体约 480KB分别承担模型讲解、操作要点整理与在线浏览用途体积轻便却信息集中。已有 276 人学习说明其在相关主题中具备一定参考价值。资料重点呈现模型注释、可修改性与收敛性处理读者可据此理解温度场与应力场的计算流程借鉴参数设置、边界条件与求解策略并在此基础上开展单道多道对比分析为增材制造热固耦合仿真提供可复用的建模思路。1. 熔池只有几十微米应力却能把零件掰弯SLM 仿真里最反直觉的一件事激光选区熔化SLM的熔池直径通常只有几十微米激光扫过一点的时间以微秒计。按直觉这么局部的热扰动应该瞬间被周围金属吸走可实际打印出来的零件常在取件那一刻肉眼可见地翘曲甚至直接开裂。温度场究竟怎么把应力场塑造到这种程度是增材制造数值仿真里最值得深入分析的一环。热固耦合的难点不在方程本身——传热方程和弹性力学方程课本里都有——而在尺度跨度熔池几十微米零件几十毫米中间还夹着固液相变、材料参数随温度剧烈变化、成千上万次热循环。COMSOL 6.1 这类多物理场平台的价值是把这些耦合关系显式摆到同一张网格上让每一次热循环留下的塑性应变都能被累加。下面按“温度场怎么算准 → 应力场怎么接上 → 铺粉怎么在模型里长出来 → 哪里最容易翻车”的顺序拆。2. 移动热源与相变潜热温度场为什么难算2.1 高斯热源在网格上“走”一步物理量发生了什么热源模型的选择直接决定了后面所有结果的可信度。SLM 里最常见的是面内高斯分布、深度方向按指数衰减的体热源它的物理含义是把激光能量在光斑半径内按高斯权重分配同时在铺粉层厚度内被金属粉末和已凝固层吸收。写成表达式就是功率除以光斑面积再乘一个归一化因子 2/π这样在无穷大平面上积分正好等于总功率乘以吸收率。吸收率这个参数在工程上通常取 0.3 到 0.5粉末床因为多次反射会比实体略高但这个值的不确定性本身就能让熔池温度摆动上百开尔文。真正的问题出在离散化上。热源在解析解里是连续移动的在有限元网格上却是“跳”着走的每个时间步高斯中心平移一个 v·dt 的距离。如果这个距离大于热源附近单元的特征尺寸高斯峰就会从一个单元直接跳到下一个单元中间那个单元永远等不到峰值温度。温度曲线上表现为明显的锯齿熔池边界忽大忽小后处理出来的熔池长度完全是随机的。常见的实现思路是给热源区域设一个硬约束v·dt ≤ h/2h 是激光路径上最小的单元尺寸。举一组实际会遇到的数光斑特征半径 r0 40 μm热源附近单元尺寸 h 15 μm扫描速度 v 800 mm/s那么时间步上限约 1.9×10⁻⁵ s。一条 10 mm 长的扫描线要走 12.5 ms也就是六百多个时间步一个 3 mm×3 mm×2 mm 的小试样几百层几百道步数轻易上百万。这就是为什么多数 SLM 仿真最后都在“减少层数、缩短扫描线”上妥协。提示网格在热源附近至少要能分辨 r0 的三分之一深度方向铺粉层厚度内至少放三层单元否则潜热释放和温度梯度都会失真。下面这段是 LiveLink for MATLAB 里定义移动体热源的大致写法方法名随版本略有差异% 移动高斯体热源面内高斯深度方向指数衰减 model.variable(var1).set(P_laser, 200[W]); % 激光功率 model.variable(var1).set(eta, 0.35); % 粉末床有效吸收率 model.variable(var1).set(r0, 40[um]); % 光斑特征半径 model.variable(var1).set(v_scan, 800[mm/s]); % 扫描速度 model.variable(var1).set(d_opt, 35[um]); % 光学穿透深度 model.variable(var1).set(t0, 0[s]); % 起始时刻 % 热源中心随时间的 x 坐标y 为道间距方向 model.variable(var1).set(x_src, 0[mm] v_scan*(t-t0)); model.variable(var1).set(q_laser, ... 2*P_laser*eta/(pi*r0^2)*exp(-2*((x-x_src)^2y^2)/r0^2)*exp(-z/d_opt));q_laser的量纲是 W/m³直接作为广义源项塞进传热方程的右端。d_opt不是随便给的数它决定了能量在深度方向分给几层单元取值和铺粉层厚同量级时热源才真正是“体”的而不是“面”的。把d_opt设得过小模型退化成面热源熔池深度会被系统性低估设得过大激光会“穿透”到基板上熔池深宽比完全失控。2.2 潜热为什么必须用表观热容而不是直接给焓熔化潜热的量级经常被低估。以常见的不锈钢粉末为例比热容约 700 J/(kg·K)熔点附近潜热约 270 kJ/kg。如果把这 270 kJ/kg 摊到一个只有 20 K 宽的相变区间里等效比热容要额外增加 13500 J/(kg·K)——是固相比热的近二十倍。忽略潜热的效果不是“稍微热一点”而是熔池温度虚高、凝固阶段冷却速率被高估、温度梯度被放大最终算出来的残余应力可能比实际高出一大截。处理相变有两条路。焓法把内能写成温度的函数 H(T)导热方程改写成 ∂H/∂t 的形式物理上最干净但 COMSOL 里要自己构造 H(T) 并处理它的反函数实操上麻烦。更常见的是表观热容法在相变区间内给比热容叠加上一个与液相分数导数成正比的项数学上和焓法等价但形式上仍然是一个 Cp(T) 的表达式直接填进材料节点就行。% 液相分数用平滑阶跃避免导数在 Ts / Tl 处不连续 model.variable(var1).set(dT_ph, 20[K]); % 相变区间展宽 model.variable(var1).set(phi_l, ... flc2hs((T-(Ts-dT_ph/2))/dT_ph, dT_ph)); % 0 - 1 平滑过渡 model.variable(var1).set(Cp_eff, ... Cp_solid(T) Lf*d(phi_l,T)); % 表观热容flc2hs的第二个参数是过渡宽度它的作用是把这个本来是阶跃的函数磨圆让d(phi_l,T)有界。这里有个纯粹的数值取舍dT_ph取太小比如 1 K表观热容的峰值会高到让雅可比矩阵条件数爆炸牛顿迭代开始反复被截断取太大比如 100 K相变被抹平到接近熔点以下的温度区间熔池边界位置提前熔池尺寸偏大。实践里 20 到 50 K 是个比较稳的窗口具体值还要和材料熔点、时间步一起调。Ts和Tl也不该取手册上的平衡值。SLM 的冷却速率在 10⁵ 到 10⁶ K/s 量级凝固过程中的固液相线间隔会被显著拉宽有些合金的实际凝固区间比平衡相图宽出几十开尔文。用平衡值算出来的熔池凝固时间偏短热影响区的温度历史会失真。2.3 材料参数随温度跳变带来的非线性热导率、比热容、密度、热膨胀系数、屈服强度这些参数在熔点附近全部是剧烈变化的。热导率从固相的十几 W/(m·K) 到液相略有上升但液相里真正的传热机制是对流——马兰戈尼效应把熔池里的热量从中心往边缘搬。不显式求解流场的话工程上惯常的做法是把液相热导率人为放大 2 到 3 倍用一个等效导热系数把对流的效果折进去。这个处理当然不严格但它把熔池形状拉回到和实验接近的量级比完全不管要靠谱得多。真正的麻烦是数据缺失。商业材料库给出的热物性通常只到 1000 K 到 1500 K而熔池峰值温度可能到 2500 K 以上。外推的方式有好几种——线性延长、取常数值、或者用液相的值硬顶——不同选择对熔池过热度的影响是几百开尔文的差异。我的习惯是固相区严格用实测数据做分段线性插值熔点以上取常数并明确记录这个假设因为那个区间本来就没什么可信数据与其编一条曲线不如承认它的不确定性。参数跳变还会直接打击求解器。热导率在固液相之间可能有两倍以上的差别如果材料节点里直接写成阶跃函数每次牛顿迭代跨过相变温度都会让残差剧烈变化。用和潜热同样的思路所有在熔点附近有跳变的参数都该用flc2hs平滑一下过渡宽度和dT_ph对齐。2.4 时间离散隐式步长为什么不能贪大传热方程是抛物型的隐式格式本身是无条件稳定的所以直觉上时间步应该可以随便取。问题出在“稳定”和“准确”是两回事。移动热源经过某个位置的时间尺度是 r0/v对上面那组参数就是 5×10⁻⁵ s 量级如果时间步取到 10⁻³ s热源一步跨过二十个光斑那个位置根本感受不到峰值温度。温度历史被严重平滑峰值温度偏低冷却速率偏离真实值。COMSOL 里的瞬态求解器默认用自由步长加 BDF 时间积分内部会根据局部误差估计自动调整。这个机制在纯传热问题里工作得不错但一旦耦合了相变和随温度变化的材料参数误差估计会变得非常保守步长被压到极小计算卡死。更可控的做法是给求解器设一个输出时间列表同时限制内部最大步长让求解器在“必须输出”和“可以大跨步”之间有个明确的边界。还有一个容易被忽略的事件激光开启和关闭的瞬间。这两个时刻热源功率是阶跃变化的求解器如果恰好跨过这个时刻会把它当成一个光滑变化处理结果就是开头几步的温度完全不对。手动在这些时刻插入输出点或者用事件接口显式标记都能避免这个问题。3. 热-固耦合为什么顺序耦合成了工程默认解3.1 全耦合与顺序耦合的自由度规模差从方程上看热固耦合是双向的温度场通过热应变影响位移场位移场通过变形做功反过来影响温度场。COMSOL 里确实可以把传热和固体力学放进同一个全耦合求解器一次性解出所有自由度。问题在于规模。三维传热每个节点一个自由度三维固体力学每个节点三个。十万节点的模型全耦合系统的规模是四十万而且雅可比矩阵是非对称的不能再用对称求解器那套省内存的技巧。更致命的是条件数传热矩阵的量级和弹性矩阵的量级差了好几个数量级放进同一个矩阵里预条件子很难做好。顺序耦合的假设是力学变形对传热的影响可以忽略。这个假设在 SLM 里几乎总是成立——熔池的流动和传热由激光能量主导固态区的变形量在微米级对热传导路径的影响小到测不出来唯一的例外是零件和基板之间出现间隙或接触状态变化。所以在工程实践里先解一个纯瞬态传热把每个时间步的温度场存下来再把温度作为体载荷插进力学场是最常见的路线。代价是插值误差。两套网格如果不一致温度要从传热网格映射到力学网格上映射过程必须用保守型插值否则总能量不守恒局部温度会凭空多出来或少掉。COMSOL 的耦合算子默认会处理这件事但如果两套网格疏密差别很大比如传热网格在热源附近细到 10 μm、力学网格为了省算力粗到 100 μm映射出来的温度梯度会被严重平滑热应力随之被低估。3.2 热应变怎么进入力学方程小变形框架下总应变分解成弹性、塑性、热应变三部分热应变张量是各向同性的等于热膨胀系数乘以温度相对参考温度的差值。这里面有两个坑。第一个是热膨胀系数的定义。材料手册给的 α 通常是相对于室温的割线值也就是 (L(T)-L(293K))/(L(293K)·(T-293K))而本构方程里需要的是切线值 dL/dT。这两者在温度跨度大的时候差别不小。COMSOL 的热膨胀节点默认按割线处理如果你的材料数据是切线值直接填进去会系统性偏差。第二个是参考温度Tref。它不是一个随便填的数物理上对应“零件处于无应力状态”的那个温度。如果仿真从室温粉末床开始Tref就该是室温如果工艺里有基板预热Tref应该取预热温度否则整个零件从第一步就带着一个均匀的热应变残余应力里混进一个本不该有的偏置。% 顺序耦合传热先算完力学场读入温度历史 model.physics(solid).feature(lemm1).feature(the1).set(Tref, 293.15[K]); model.physics(solid).feature(lemm1).feature(the1).set(alpha, alpha_solid(T)); % 温度来自传热场的插值而不是重新解一次 model.physics(solid).feature(lemm1).feature(the1).set(T, T_ht);T_ht是传热场映射过来的温度变量力学场只把它当成一个已知的载荷不参与求解。塑性部分用温度相关的屈服面和强化曲线SLM 的循环加热冷却特征意味着等向强化不够用——材料在反复的拉压交替下会出现包辛格效应随动强化模型比如 Chaboche 的多背应力形式更贴近实际但参数标定成本高多数工程仿真还是用等向强化加一个折减的屈服强度来近似。3.3 约束条件决定了残余应力的空间分布残余应力不是材料属性是“约束”的产物。同样一条扫描路径基板夹得死和夹得松算出来的应力云图能差出一倍。最常见的错误是把基板底面全固定。这个边界在数学上没有问题物理上却是错的真实工艺里基板靠几个螺栓压在平台上允许一定的转动和面内位移。全固定相当于给零件加了一个无限刚性的约束热膨胀被完全挡住应力被人为放大。合理的做法是只约束刚体运动的六个模态——比如底面三个点分别约束法向、两个切向、再加两个点约束转动——剩下的让它自由变形。基板本身的厚度也影响巨大。薄基板在热循环下会跟着翘吸收掉一部分热应变厚基板相当于刚性平台热量往下传得快冷却速率高应力反而更大。仿真里如果把基板简化成一层薄壳或者干脆省略结果基本不可用。取件之后应力释放、零件翘曲是另一个必须单独算的工况。做法是在第一个 study 算完之后把基板区域的单元钝化掉再跑一个静力平衡看零件在新的约束状态下的变形。这一步经常被跳过但它恰恰是现场最关心的那个量。4. 生死单元与多层多道铺粉过程怎么在模型里“长”出来4.1 单元激活的数值本质SLM 是一层一层长出来的仿真里要让材料“凭空出现”。COMSOL 的增材制造模块提供了单元激活功能底层机制其实很朴素给每个单元定义一个激活指示变量未激活时为 0激活后为 1然后把材料属性——热导率、比热容、密度、弹性模量——全部乘上这个指示变量。这里有个绝对不能踩的坑比热容和密度不能真的乘 0。传热方程里的瞬态项是 ρ·Cp·∂T/∂t如果这两个都变成 0那一块区域的矩阵对角线就是零整体矩阵奇异求解器直接报错。惯常做法是乘一个极小值比如 10⁻⁶让矩阵保持非奇异同时它对结果的影响小到可以忽略。第二个细节是激活瞬间的温度初始化。新激活的单元如果初始温度取 0 K 或者绝对零度附近的默认值激活那一刻会出现一个巨大的温度跳变热流密度瞬间冲到非物理的量级时间步被砍到 10⁻¹² 秒计算直接卡死。正确做法是让新单元的初始温度继承下方已凝固层的温度或者干脆用一个平滑的斜坡在几百微秒内把材料属性从“几乎不存在”过渡到“完全存在”。第三个细节是激活顺序和网格的匹配。如果一层单元在同一时刻全部激活相当于整层粉末瞬间变成金属并瞬间被加热熔池会变成一个扁平的方块。真实过程是激光逐点扫描所以激活必须跟着热源位置走单元激活的判断条件通常写成热源中心到单元中心的距离小于某个阈值或者按扫描路径参数化地给出激活时刻表。4.2 扫描策略对残余应力的影响扫描策略是 SLM 里少数几个“不改材料、不改设备、只改路径就能显著改善残余应力”的杠杆。常见的有长边扫描、短边扫描、棋盘格分区、层间旋转 67°。长边扫描的扫描线长温度梯度沿着线的方向拉得很开冷却时收缩量大端部的应力集中明显。短边扫描把长线切成短段单段的热输入小温度梯度缓和但道数变多道间的重熔区增加。棋盘格把一层切成若干小岛每个岛单独扫描岛与岛之间留出搭接区好处是每一小块的热积累都很有限坏处是岛边界处容易出现未熔合和孔洞。层间旋转 67° 是工业上最常见的做法它的作用不是让应力均匀而是打断热流路径的周期性。如果每层都沿同一方向扫热量会沿着扫描线方向一直往前传形成一条贯穿多层的“热通道”冷却后这条通道上下的收缩量和其他区域不一致容易诱发起裂。旋转一个非特殊角度67°、45° 都是常见值前者因为不与层内道间距形成整数倍关系而更受欢迎能让热流方向在层间被打散。仿真里复现这些策略的难点在于路径的表达。把扫描路径写成参数化曲线用弧长参数驱动热源位置是最干净的做法如果路径是离散的点序列就要用插值函数或者分段线性函数把时间映射到坐标上段与段之间的衔接处容易出现速度不连续导致热源在拐角处“停顿”局部过热。这是很多仿真结果里熔池在扫描线端点异常膨胀的原因。4.3 时间步、层间冷却与算力的三方拉锯一层扫描的时间和层间冷却的时间通常是同一个量级甚至后者更长。真实工艺里铺粉加刮平要几秒到十几秒而扫描一层只要几十毫秒。如果老老实实把冷却时间也算进去计算量的绝大部分都花在“温度几乎不变”的阶段性价比极低。常见的加速手段有三种。第一种是把层间冷却时间压缩比如真实 10 s 压到 0.1 s代价是热积累被低估层间温度偏低残余应力偏小。第二种是用“跳层”策略算几层跳几层跳过的层用插值补上这对温度场的宏观分布影响不大但对熔池尺度的细节完全失真。第三种是只算一个代表性体积单元把多道多层的边界条件用对称或者周期性条件近似这种做法能保住熔池细节但没法回答零件级的变形问题。这三种手段本质上都是在“熔池尺度”和“零件尺度”之间做选择没有哪个是普适的。选哪个取决于你要回答什么问题关心熔池形貌和孔隙缺陷就得保熔池、牺牲层数关心整体翘曲和残余应力分布就得保零件、牺牲熔池细节。这个取舍在第 6 章还会再展开。5. 排查清单SLM 热固耦合最容易踩的五个坑5.1 坑一熔池凭空消失温度曲线爬不到熔点现象很典型求解顺利完成后处理一看最高温度只有 1200 K连熔点都没到熔池根本不存在。原因通常有两个。一是热源表达式里的单位错位比如光斑半径填了40[um]但公式里忘了平方或者功率填成200而没带[W]COMSOL 会按 SI 基本单位解释成 200 W看起来没错但如果几何单位是 mm 而表达式里混用了就会出现量级偏差。二是吸收率和穿透深度共同作用把能量摊得太薄——d_opt设成几百微米时同样功率分布在几十倍的体积里峰值温度自然上不去。解决办法是先做一个量级校验用解析的 Rosenthal 解或者简单的一维估算算一下给定功率和扫描速度下熔池大致能到多少度和仿真结果对一遍。差一个数量级一定是单位或者公式的问题差百分之二三十才是模型精度的问题。5.2 坑二温度曲线锯齿振荡时间步被砍到极小现象是某个节点的温度历史不是光滑的峰而是一串尖刺求解日志里时间步长一路被压缩到 10⁻¹⁰ 秒以下进度条几乎不动。原因多半是时间步和网格不匹配热源一步跨过了多个单元高斯峰在网格上“跳跃”每个单元交替经历“被峰扫到”和“没被扫到”。另一种可能是材料参数里存在阶跃比如热导率在熔点处直接从一个值跳到另一个值牛顿迭代每次跨过这个温度都会剧烈震荡。解决办法分两步先把最大时间步限制在h/(2v)以内确认锯齿消失如果还有残留把所有在熔点附近跳变的材料参数用flc2hs平滑过渡宽度不小于相变区间的展宽。这两步做完九成以上的锯齿问题都能解决。5.3 坑三迭代不收敛报“非线性求解器不收敛”现象是求解器在某个时刻反复回退步长最后还是失败日志里写的是牛顿迭代达到最大次数。原因往往出在激活瞬间或者相变区间。新单元激活时材料属性从 10⁻⁶ 跳到 1如果激活的单元数量在同一时间步里很多刚度矩阵和传导矩阵会同时发生大变化牛顿迭代的初值离解太远。另一个常见原因是塑性迭代和传热迭代在同一个耦合步里互相打架尤其是热膨胀系数取切线值的时候。解决办法是给激活加时间斜坡把属性变化摊到几十个时间步上给非线性求解器加阻尼把更新量的步长限制在合理范围内实在不行就把传热和力学彻底分开成两个 study中间只传温度场不让它们在同一个牛顿循环里互相牵制。5.4 坑四位移全为零或者零件直接“飞走”位移全零通常意味着力学场的约束给多了——所有边界都被固定热应变无处释放位移自然是零但应力大得离谱。零件“飞走”则相反约束不足刚体模态没被消掉求解器解出一个巨大的刚体位移。这两种现象的根因是同一个约束条件的物理意义没想清楚。合理的做法是先问“真实工艺里这块材料被什么挡住了”只保留那些真实存在的约束剩下的用弱弹簧或者刚体运动抑制来消掉数值上的奇异性。刚体运动抑制不会引入额外的应力这是它比全固定优越的地方。5.5 坑五结果看起来很合理数值却差一个量级最危险的一种情况云图颜色漂亮熔池形状也对但残余应力的数值和实测差了一倍。这种“看起来对”的错误最难查。常见根源有三个吸收率取值偏离实际粉末床的真实吸收率随粉末粒径、表面氧化状态变化很大0.3 和 0.5 能让应力差 40%参考温度设错导致整体带一个偏置应力塑性模型里用的屈服强度是室温值而熔池附近温度高达上千度实际屈服强度只有室温的十分之一。核查办法是固定看三个量熔池的最高温度和尺寸、冷却速率对凝固组织有直接影响、以及取件后零件的最大翘曲量。这三个量和实验或文献一对比量级对不对立刻现形。只看应力云图的颜色分布是最容易自欺欺人的做法。6. 设计取舍二维切片、单道外推还是硬算整个零件仿真做久了会发现SLM 建模的核心矛盾不是“会不会算”而是“算哪个尺度”。熔池尺度需要微米级网格和微秒级时间步零件尺度需要毫米级网格和秒级时间步两者相差三到四个数量级任何试图同时覆盖的模型都会被计算量压垮。一个实用的判断方式是先明确模型要回答的问题属于哪一类。关心熔池形貌、球化、匙孔、未熔合就做单道或者几道的二维/三维小模型网格细到 5 μm时间步细到微秒层数不重要。关心残余应力分布和零件翘曲就用“等效热源”或者“逐层激活 粗网格”的方式把熔池细节抹平用体积热源或者分层的温度边界条件代替重点放在约束条件和材料本构上。关心扫描策略的优劣对比用一个中等规模的板状试样十几层、几十道网格 20 μm 量级能分辨道间搭接但分辨不了熔池内部。模型层级典型网格尺度时间步能回答的问题算不出来的东西单道熔池3–10 μm1–10 μs熔池形貌、深宽比、冷却速率零件变形、层间热积累多道多层小试样15–30 μm10–50 μs道间搭接、层间热积累、局部应力整体翘曲、基板效应零件级等效100–500 μm0.1–1 ms整体温度分布、残余应力趋势、翘曲熔池细节、孔隙缺陷这张表最有用的地方不是选哪一行而是提醒你任何一行都在放弃某些信息而这些被放弃的信息恰恰可能在真实工艺里是关键。二维切片能省掉一个数量级的算力但它隐含了“面外方向均匀”的假设对于扫描线端部、零件边缘、薄壁结构这个假设直接失效。单道外推到多道忽略了道间的重熔和热积累应力会被系统性低估。我自己的习惯是任何模型开算之前先用解析估算把熔池温度和尺寸对一遍再用一个极简的二维模型把热源参数标定好最后才上三维。这个流程看起来多花时间但它把“参数错了”和“模型选错了”这两类问题分开了——前者在二维阶段就能发现后者要到三维才会暴露早发现早止损。那次教训来自一个薄壁件为了省时间直接用了二维切片扫出来应力很小结论是工艺可行。实际打印出来壁厚方向翘得没法用。后来才想明白二维切片把面外方向的约束当成了无限刚体热应变全被压住了应力当然小。从那以后我每次建模前都强制走一遍“这个假设在哪个方向失效”的自检不再让网格疏密决定结论。希望这些经验帮到你少走一点弯路。本文还有配套的精品资源点击获取