第一次用 COMSOL 做结晶凝固仿真时我天真地以为这不过是一个带相变的传热问题。结果第一轮瞬态计算跑了两百步温度云图直接出现 -230°C残差曲线比心电图还刺激。那一刻我才意识到结晶凝固看着像纯传热骨子里却是一个多物理场问题传热、流动、溶质再分配甚至固体力学都在同一个时空尺度上互相影响。这篇文章就把我做凝固仿真的建模思路、方程处理、实操细节和仿真发散的排查经验整理出来给正在做相变仿真、或者被不收敛问题折磨的同行一个可上手的参考也希望能帮你少走几条弯路。1. 凝固相变仿真的难点不在熔点而在三个“暗礁”1.1 一个直观例子冰水相变的能量账本很多人第一次接触 COMSOL 凝固仿真反应是“不就是温度到了 0°C 就变成冰吗”。如果真这么处理结果会错得离谱。我习惯先算一笔能量账。取 1 kg 水从 1°C 冷却到 0°C需要放出的显热大约是 4.186 kJ。但水在 0°C 完全结冰时潜热是 333.5 kJ/kg——显热的 80 倍。如果把这个热量漏掉那“完全凝固需要多长时间”这个最核心的指标计算误差不是百分之几而是几倍甚至几十倍。所以模拟凝固问题第一个暗礁就是潜热怎么处理。它不是一个小修正项而是控制整个凝固时间尺度的主导项。用恒定的比热去硬算一个含相变的过程本质上是把能量账本做错了。1.2 固液前沿的对流为何不能忽略第二个暗礁是流动。凝固过程中只要有温度梯度就有密度差重力一作用自然对流就起来了。凝固前沿释放潜热会让界面附近的液体温度比远处更高于是形成局部的浮力驱动流这个流动反过来又推送热量让凝固前沿变得不均匀甚至影响晶粒生长结构。我早期犯过一个典型错误为了省事只做了纯导热模型算出来的凝固壳层非常均匀、非常“漂亮”。后来用冷水做对照实验发现实际凝固边界是不规则的角落的地方冰层明显更厚。原因就是自然对流把低温流体带到了某个区域加速了局部结晶。在 COMSOL 里做金属铸造仿真时这个现象更突出。液态铝的动力粘度大约 1.2×10⁻³ Pa·s比水还低自然对流一旦成立流速可以达到厘米每秒量级这种流动对凝固时间的影响绝不是什么二阶修正。1.3 多物理场版图传热、流动、溶质、应力谁主谁次把几个物理场摆在一起看你会发现它们之间是一个互相咬合的闭环温度场通过液相线、固相线决定固相分数相变潜热作为源项反过来再影响温度场温度梯度和浓度梯度通过浮力项驱动流动流场通过对流项输运热量和溶质凝固收缩、热应变产生应力应力又可能改变构件与模具的接触状态接触热阻变了传热边界就跟着变。这就是“多物理场的奇妙交织”的本质。做工程仿真不是所有环节都要开全而是先判断哪个物理场对你的目标问题起主导作用。只关心铸件冷却时间那传热和潜热是主演关心缩孔缩松就得把流场和糊状区流动阻力算进去关心宏观偏析溶质场必须出场关心热裂纹应力场就躲不掉了。所以我建议上手时不要一上来就把五个物理场全叠进去那大概率会得到一个发散得毫无头绪的结果。先把主演拉住再把配角一层层加进来反而是最快到终点的路径。2. 把物理写进 COMSOL控制方程与材料参数的处理技巧2.1 潜热处理等效热容法的思想与阈值COMSOL 里处理潜热最常用也最稳妥的是等效热容法也叫表观热容法。思路很简单把潜热 ΔH 摊到一个很窄的相变温度区间 [Ts, Tl] 内在这个区间内的等效比热写成cp_eff cp ΔH / (Tl - Ts)从公式能直接看出为什么“相变区间”这个参数这么敏感。用铝来算一笔账固相线 655°C液相线 660°C熔融潜热约 397 kJ/kg区间宽度 5 K那等效比热就是900 397000 / 5 ≈ 80,300 J/(kg·K)铝的基础比热才 900 J/(kg·K)相变区间的等效比热直接变成 89 倍。你会看到材料属性在相变点附近出现一个极其尖利的峰——这个峰不是物理上不存在而是数值上极度难收敛。如果区间缩到 2 K峰值会到 199,400 J/(kg·K)基础值的 220 多倍。所以我在实际项目里的经验是金属凝固的相变区间通常至少给到 3-5 K水/冰这类纯物质可以放在 0.5-1 K但在 COMSOL 里建议也在 2 K 左右用足够平滑的过渡函数去削弱尖峰。精度损失不大但收敛性会显著好转。在 COMSOL 材料定义里这个表达式可以写成cp_eff C_p0 L_ph / (T_l - T_s) * flc2hs((T - T_mid) / (T_l - T_s), 0.5)这里的 flc2hs 是 COMSOL 内置的连续平滑 Heaviside 函数用来把阶跃式的相变“抹”成连续过渡。这也是实现等效热容法最省事的手段。2.2 流动方程怎么选层流、Boussinesq 还是全可压缩流场部分大多数凝固场景都可以用层流假设。凝固过程的液态流速极低特别是在糊状区里流速基本在毫米每秒以下雷诺数很小一上来就上湍流模型属于给自己找罪受。浮力项的处理要留意。COMSOL 的层流接口默认不可压缩流动直接把密度当常数浮力当然也就被忽略了。要计入自然对流常见做法有两种第一种是 Boussinesq 近似在动量方程的体积力上添加F ρ_ref * β * (T - T_ref) * g这个近似假设密度只在重力项里随温度变化其余地方都当常数。对温差不大、密度变化幅度小的过程很有效水从室温到接近冰点差十几度用 Boussinesq 没问题。第二种是使用“非等温流”耦合节点把层流和流体传热真正联立起来再勾选考虑密度变化。凝固过程其实有固体和液体两相密度差异往往不能忽略水变冰密度从 1000 掉到 917 kg/m³铝液凝固后也会收缩。如果流动和界面形态是研究重点建议直接用 COMSOL 的“非等温流”多物理场节点它会自动把传热对流体属性的影响带进去。2.3 溶质场与偏析什么时候非加不可溶质再分配是凝固仿真里最容易被忽略、但工程影响极大的一块。如果只算温度场和流场你只能回答“哪里先凝固”回答不了“哪里成分偏析了”。宏观偏析的经典描述是 Scheil-Gulliver 方程Cs k0 * C0 * (1 - fs)^(k0 - 1)其中 k0 是平衡分配系数。这个式子表达了一个物理事实先凝固的固体它的溶质含量和母液不一样剩下的液相会被不断富集凝固到最后的地方成分往往最“偏”。在 COMSOL 里实现宏观溶质输运要添加“稀物质传递”接口并额外加一个与固相分数变化率成正比的源项模拟凝固过程中液相浓度被排斥富集的过程。很多教材不会讲这个源项怎么加我的经验是先在传热接口里用变量方式导出固相分数 fs然后在稀物质传递的“反应”节点里写R (k0 - 1) * c * d(fs, t) / (1 - fs eps)配合合适的初始浓度和扩散系数就能在宏观尺度上看到通道偏析这类现象的雏形。但老实说这种连续介质方法对糊状区内的枝晶间偏析并不够精确。真要做微观枝晶生长得换相场法。工程上评估“哪里容易偏析”用稀物质传递加 Scheil 级别的源项已经能给出很好的倾向性判断。2.4 材料参数的常见翻车操作做 COMSOL 仿真的朋友应该都知道参数翻车导致的发散和假结果比求解器问题还常见。我自己踩过和见过别人踩的典型有这么几种用室温固态的物性替代液相物性。最典型的就是导热系数和粘度铝固态导热系数 237 W/(m·K)液态只有 100 左右差了一倍多硬用固态值算出来的温度场必然偏。潜热单位搞错。COMSOL 里材料属性是严格按单位处理的如果从论文里看到 J/mol没有换算成 J/kg 就填进去潜热直接放大几十倍结果自然是温度平台长时间不动。液相线和固相线填反。看起来很低级但我在多组分合金里确实犯过一不留神把 Ts 填成了 Tl等效热容被分配到错误的温度区域整个凝固顺序完全乱掉。相变区间设太窄。这个问题在 2.1 已经说过我这里再强调一次纯理论区间 1 K 看着严谨数值上基本是给自己挖坑。工程仿真不是一个纯粹的“精确”问题模型的数值可求解性也是方案的一部分。忽略接触热阻。如果是砂型铸造或者有模具铸件和模具之间不是理想接触传热系数可能只有几百 W/(m²·K)而不是理想接触下的“连续温度”。这个边界条件怎么设对凝固时间的影响可能比材料导热还大。3. 从几何到求解器搭建凝固模型的实际操作流程3.1 几何与网格先根据热扩散尺度估算搭模型之前先别急着画图而是估算一下这个问题的热扩散特征长度。扩散尺度公式很简单L ≈ sqrt(α * t)其中 α 是热扩散率t 是你关心的物理时间。以冰蓄冷为例水的热扩散率约 1.34×10⁻⁷ m²/s取 1 小时的凝固时间扩散长度大概是 sqrt(1.34e-7 × 3600) ≈ 0.022 m也就是 22 毫米。那你在冷壁面附近、凝固前沿会经过的区域网格就不能比这个尺度相差太大否则界面推进一次就跨过好几个单元温度场和固相分数的分辨率完全不够。我的经验是在扩散长度内至少布置 10 到 20 个单元也就是说这个尺度下界面附近网格要 1-2 mm 量级。而一个 200 mm 的铸铝件铝液热扩散率约 3.56×10⁻⁵ m²/s冷却 60 秒扩散长度也接近 46 mm。所以粗算下来界面附近网格 3-6 mm 才说得过去其余区域可以放宽到 10 mm 以上再用边界层网格照顾壁面。几何维度的选择也很关键。能对称就对称能用二维或轴对称就不要一上来开三维全尺寸模型。很多圆坯、板材问题用二维轴对称模型就能拿到足够工程精度的结果计算量却小一两个数量级。先二维跑通物理过程再决定要不要升级到三维这个顺序永远不会浪费。3.2 多物理场节点的装配流体传热、层流与溶质的联动COMSOL 里做凝固仿真我的默认起点是“流体传热”加“层流”。在物理场接口里分别添加之后直接在“多物理场”节点下自动生成“非等温流”耦合这一步能省很多事它会正确地处理流速对温度方程的对流输运以及温度对流体属性的影响。需要提醒的是传热接口里的对流项不是自动全部打开的。“流体传热”在纯固体区域会自动退化掉对流在流体区域才有 u·∇T 项。如果你的模型里有一个区域是固体金属它就不应该有流场所以你需要用一个“指定位移/变形”或干脆把流场速度在固相区强制归零。最省事的做法是我稍后要讲的“糊状区阻力法”用一个随固相分数变化的阻力项让速度在凝固后自动变成零。溶质场的接入在物理场列表里选“稀物质传递”并指定扩散系数和对流速度来自层流。这样三套方程就联动起来了传热决定固相分数固相分数通过阻力项和源项分别作用回动量和溶质方程。3.3 凝固前沿处理动网格的边界感和“糊状区阻力法”很多人想到追踪凝固界面第一反应是 COMSOL 的“移动网格”或“变形几何”。我的建议是宏观大尺度凝固问题不要用移动网格除非你对界面几何有极其明确的定义和把握。原因很简单真实凝固不是一个尖锐的几何边界而是一个宽度可观的模糊区。糊状区里液固共存流动阻力连续变化。硬要用动网格去追踪一条“边界”不仅几何不好定义网格变形到一定程度还会翻折导致雅可比行列式变成负值整个求解直接崩掉。我在铸造仿真里最常用的替代方案是湖状区阻力法F_drag -A_mushy * (1 - f_l)² / (f_l³ eps) * v其中 f_l 是液相分数A_mushy 是一个人工阻力系数典型的量级在 10⁵ 到 10⁸ kg/(m³·s)。这个式子表达的是当液相分数接近 1 时阻力接近 0速度自由发展当凝固程度升高、液相分数趋近 0 时阻力趋于无穷大速度被“冻住”。我在层流动量方程的源项里加这个表达式让流动区在凝固前沿自然消失完全不用动网格数值稳定性也远好于几何边界追踪。如果想要更连续还可以用 flc2hs 平滑这个过渡让阻力变化不会陡成一面墙。3.4 瞬态研究设置BDF、全耦合与步长上限求解器设置是另一个大坑我见过太多人拿着默认设置直接跑相变问题然后对着发散报错发呆。下面这张表是我经过多次实际项目后一套比较稳的默认起点按网格规模可以微调设置项推荐值原因时间步进BDF最大阶 2-5对非线性瞬态问题稳定适合相变这种强非线性初始步长自动或 1e-3 s防止一开局就跨过相变尖峰最大时间步长相变区间宽度 / 最大温变速率保证每个相变平台里有足够采样点全耦合阻尼0.7-0.8压低迭代振荡收敛慢一点但稳非线性求解器最大迭代数25默认可能偏小复杂耦合下不够用线性求解器PARDISO 或 MUMPS直接求解器在强非线性下更稳内存允许就优先最大时间步长这条值得单独解释如果你设置了 5 K 的相变区间而局部冷却速率最猛的地方有 2 K/s那么一个时间步最多走 2.5 秒就会跨完整个相变峰。为了保证在等效比热峰值附近有足够多迭代采样点最大步长最好控制在 0.5-1 s。这一步是在“精度”和“不收敛”之间做的最有效平衡之一。全耦合阻尼也是一个关键旋钮。COMSOL 默认的全耦合阻尼比较激进相变这种强非线性场景很容易振荡。我把阻尼降到 0.7 之后很多原来看起来不可能收敛的模型都安静下来了。4. 仿真发散排查从异常温度到残差暴走的完整定位链4.1 发散的症状残差、负温度和 NaN 的解读先说说发散之前你会看到什么。我总结过几个典型症状残差曲线反复振荡下不去。这说明非线性迭代在有效热容尖峰附近来回弹跳像过阻尼弹簧停不下来。常见原因就是相变区间太窄、时间步太大或者全耦合阻尼太小。温度出现负几百摄氏度。抱错发生在传热方程里负温度就是迭代过程某一步跨过了极端值然后被“放大”到全解域。很多时候是因为潜热等效比热突然变成基础比热的几十上百倍导致线性化失败。NaN 突然出现。这是最彻底的崩溃矩阵已经出现了无穷大或 0/0。常见原因湖状区阻力法的分母 (f_l³ eps) 里 eps 设得太小当 f_l 0 时分母趋近 0阻力项直接爆掉或者稀物质传递的浓度在迭代中变成负值扩散方程取了对数之类运算。报错信息本身往往只有一句“Failed to find consistent solution”或“Divergence detected in nonlinear solver”价值不大真正有用的是看“哪一个物理场先起药效”。4.2 我的排查顺序按“物理-网格-数值”三层推进遇到发散不要慌也不要瞎调。我固定的排查顺序是三层第一步砍物理场。先把流场关掉只做纯导热加等效热容看看这个最简模型能不能跑通。如果最简模型都不稳那问题大概率在材料参数或网格不在耦合。如果最简模型稳了再把层流加回去问题大概率出在流动与传热的耦合步进上。第二步观察相变区和时间步。检查等效热容曲线是否平滑把相变区间从 2 K 扩到 5 K 再试。同时手动限制最大时间步长到 0.1 s先跑过最开始几秒物理时间把初始瞬态尖峰消化掉再放开步长。第三步调求解器。全耦合阻尼降到 0.5最大迭代数调 25不行就换成分离式求解器把“流体传热”和“稀物质传递”分开求解。面这种多物理场强耦合问题分离式有时反而比全耦合更容易收敛只是每个物理场之间的信息交换会滞后一步。我做过一个小统计身边的同事遇上的凝固仿真发散大约 70% 是相变区间过窄或时间步过大20% 是材料参数单位搞错只有 10% 才轮得到真正的求解器配置问题。所以排查顺序一定是先物理再数值别一上来就去动求解器。4.3 真正收敛的标准能量守恒、网格无关性和物理指标“软件没报错”不等于“收敛了”。这是很多初学者最大的误区。我自己判断一个凝固仿真是否真正收敛至少看三样东西。第一是能量守恒。COMSOL 后处理里可以对整个域做热量积分的验证边界流入的总热量加上初始内能应该等于最终内能加上累计相变热。如果这账对不上五个百分点以上结果基本不能信。第二是网格无关性。把界面附近网格加密 1.4 倍把最大网格尺寸从 2 mm 缩到 1.4 mm再看关键指标——比如某个节点温度达到固相线的时间——变化是否小于 1%。如果变化很大说明网格还没足够细继续加密。第三是物理合理性。凝固时间如果比理论估计或实验值差了数量级那不管残差曲线多漂亮都是白搭。我一直强调仿真器的友好界面会给人“算出来就是对的”的错觉但真实工程世界里判断结果合理性的能力才是核心技能。5. 后处理里的工程信号固相分数、温度回升与缺陷倾向5.1 固相分数判断凝固进度比温度更可靠很多人看云图习惯只看温度但在凝固仿真里真正有意义的是固相分数 fs。温度只能告诉你“有没有到相变温度”而 fs 告诉你“这个位置已经完全凝固还是仍处于糊状区”。在 COMSOL 里我通常定义一个变量fs flc2hs((Tl - T) / (Tl - Ts), 0.5)然后在后处理里画 fs 的等值线或者做全域平均。当这个平均值等于 1 时代表区域全部凝固完成。这个“全部凝固时间”比“某个点降到多少度”更能代表一个铸件的整体凝固节奏。判断一个铸件哪里是最后凝固位置就画 fs 最小的区域往往就是热节的位置。热节处最容易形成缩孔缩松如果不做补缩措施问题大概率会出现在那儿。5.2 潜热释放与温度回升别被大时间步骗过去一个很有意思的物理现象是温度回升。纯物质凝固时当局部形核开始释放潜热如果散热速度跟不上潜热释放速度局部温度不仅不下降反而会往上回升一小截。在温度曲线上这个区域会呈现一个平台甚至一个微小的驼峰。我见过很多仿真报告里的降温曲线是一条特别光滑的直线完全没有这个平台那基本可以断定时间步设得太大把温度回升抹掉了。这个平台的意义在于它代表了潜热释放与外部散热的拉锯过程抹掉它凝固时间的计算就有偏差。我的习惯是额外输出 dT/dt 的曲线。温度回升在 dT/dt 上表现为一个尖锐的回峰一眼就能看出模型是否捕获了这个关键物理细节。5.3 缩孔、冷隔与偏析在仿真结果里怎么预告仿真做出来之后怎么把这些场变量翻译成“这里可能会出缺陷”缩孔倾向看压力场和最后凝固区域。如果某块区域液相分数最后才到 1同时该处压力下降、流动补缩路径被切断那么缩孔缩松的风险就非常高。我常用“最后凝固时间等值线图”加“压力最小值位置”叠加来判断补缩通道是否畅通。冷隔倾向看两股凝固前沿的汇合状态。铝合金压铸里经常出现两股低温熔体相遇温度已经降到液相线以下界面没能融合形成氧化皮般的冷隔。在仿真里我会看两个熔体前沿到达同一位置时各自的温度如果都低于液相线冷隔风险就大了。这个不能只看最终温度场要看不同时刻的前沿位置。宏观偏析则看溶质场。稀物质传递的浓度云图会出现明显的高浓度区就是被富集液相冲刷过的通道。这类结果在宏观尺度上已经能指导浇注系统和冒口位置的调整。6. 结晶凝固仿真的迁移能力从金属铸件到电池热管理6.1 电池热管理中的 PCM等效热容法的同款迁移做完金属凝固模型后你会发现这套打法几乎可以直接迁移到电池热管理里的相变材料PCM仿真。石蜡这类 PCM 的潜热大约 180-220 kJ/kg比金属低一些但相变区间往往比金属宽所以等效热容峰不会那么尖数值稳定性反而好。锂电池散热中常见的设计是在电池模组间隙填充石蜡/石墨复合 PCM靠相变吸收脉冲热量。我在 COMSOL 里做的第一个电池热管理仿真用的就是和铸件凝固几乎一模一样的代码结构有效热容法加自然对流只是相变点换成了 28-42°C 的温区。区别在于PCM 凝固过程中会有体积膨胀密度变化更大这时候 Boussinesq 近似可能需要换成弱可压缩流动另外 PCM 的导热系数往往很低只有 0.2-0.3 W/(m·K)所以要加高导热填料来补材料的等效导热系数必须重测不能指望纯石蜡的数据库值有多靠谱。6.2 焊接与增材制造里的极端凝固条件如果你想碰更极端的场景焊接熔池和增材制造熔池是同一个体系里的“偏科生”。熔池尺寸小冷却速度高达 10³-10⁶ K/s温度梯度极大温度回升现象可能非常短暂如果时间步不够小完全看不见。这里有两个和前面完全不同的重点。一是热源模型高斯热源或双椭球热源需要写成随时间和位置变化的体热源函数二是熔池内的对流驱动机制不只是浮力还有表面张力梯度驱动的马兰戈尼对流后者在激光焊接里经常占据主导。对于这类问题我的建议是千万别一上来就开全物理场。先做纯导热等效热容模型跑出熔池形状再逐步加上流场和表面张力每加一个物理场都重新验证一遍温度场。这个循序渐进的思路在凝固这类强非线性问题上比任何高级求解器设置都管用。6.3 案例库拆解与个人学习路径COMSOL 自带的案例库和官网案例下载区其实是一个很被低估的学习资源。与凝固、相变相关的最直接案例有相变材料储热、材料熔化、铸造凝固等你不需要照着抄重点拆解它的“多物理场”节点是怎么装配的、材料属性里相变是怎么定义的、求解器的阻尼和步长是怎么设的。我每次接触一个新物理现象第一件事都是去案例库找最接近的那个模型把它从头到尾拆一遍尤其是看“变量”和“表达式”页签里的写法往往能发现很多官方文档里没写透的细节。一个比较顺畅的学习路径是先用流体传热接口把等效热容法的水结冰模型跑通然后加层流做自然对流再引入稀物质传递做溶质输运最后加固体力学看热应力。每一步只改动一个物理场其余保持不动。这样一旦出现发散你可以立即定位到是新引入的那个物理场带来的问题。现在我做任何凝固或相变课题第一步永远不会变先跑纯导热的等效热容模型拿到固相分数随时间变化这条曲线再往里面加流动、溶质和应力。这套流程看着不够炫但它是所有复杂多物理场模型的地基。地基打稳上面的大楼才不至于在第一次迭代时就塌成一片残差曲线。
COMSOL凝固仿真全攻略：从等效热容法到多物理场收敛排查
第一次用 COMSOL 做结晶凝固仿真时我天真地以为这不过是一个带相变的传热问题。结果第一轮瞬态计算跑了两百步温度云图直接出现 -230°C残差曲线比心电图还刺激。那一刻我才意识到结晶凝固看着像纯传热骨子里却是一个多物理场问题传热、流动、溶质再分配甚至固体力学都在同一个时空尺度上互相影响。这篇文章就把我做凝固仿真的建模思路、方程处理、实操细节和仿真发散的排查经验整理出来给正在做相变仿真、或者被不收敛问题折磨的同行一个可上手的参考也希望能帮你少走几条弯路。1. 凝固相变仿真的难点不在熔点而在三个“暗礁”1.1 一个直观例子冰水相变的能量账本很多人第一次接触 COMSOL 凝固仿真反应是“不就是温度到了 0°C 就变成冰吗”。如果真这么处理结果会错得离谱。我习惯先算一笔能量账。取 1 kg 水从 1°C 冷却到 0°C需要放出的显热大约是 4.186 kJ。但水在 0°C 完全结冰时潜热是 333.5 kJ/kg——显热的 80 倍。如果把这个热量漏掉那“完全凝固需要多长时间”这个最核心的指标计算误差不是百分之几而是几倍甚至几十倍。所以模拟凝固问题第一个暗礁就是潜热怎么处理。它不是一个小修正项而是控制整个凝固时间尺度的主导项。用恒定的比热去硬算一个含相变的过程本质上是把能量账本做错了。1.2 固液前沿的对流为何不能忽略第二个暗礁是流动。凝固过程中只要有温度梯度就有密度差重力一作用自然对流就起来了。凝固前沿释放潜热会让界面附近的液体温度比远处更高于是形成局部的浮力驱动流这个流动反过来又推送热量让凝固前沿变得不均匀甚至影响晶粒生长结构。我早期犯过一个典型错误为了省事只做了纯导热模型算出来的凝固壳层非常均匀、非常“漂亮”。后来用冷水做对照实验发现实际凝固边界是不规则的角落的地方冰层明显更厚。原因就是自然对流把低温流体带到了某个区域加速了局部结晶。在 COMSOL 里做金属铸造仿真时这个现象更突出。液态铝的动力粘度大约 1.2×10⁻³ Pa·s比水还低自然对流一旦成立流速可以达到厘米每秒量级这种流动对凝固时间的影响绝不是什么二阶修正。1.3 多物理场版图传热、流动、溶质、应力谁主谁次把几个物理场摆在一起看你会发现它们之间是一个互相咬合的闭环温度场通过液相线、固相线决定固相分数相变潜热作为源项反过来再影响温度场温度梯度和浓度梯度通过浮力项驱动流动流场通过对流项输运热量和溶质凝固收缩、热应变产生应力应力又可能改变构件与模具的接触状态接触热阻变了传热边界就跟着变。这就是“多物理场的奇妙交织”的本质。做工程仿真不是所有环节都要开全而是先判断哪个物理场对你的目标问题起主导作用。只关心铸件冷却时间那传热和潜热是主演关心缩孔缩松就得把流场和糊状区流动阻力算进去关心宏观偏析溶质场必须出场关心热裂纹应力场就躲不掉了。所以我建议上手时不要一上来就把五个物理场全叠进去那大概率会得到一个发散得毫无头绪的结果。先把主演拉住再把配角一层层加进来反而是最快到终点的路径。2. 把物理写进 COMSOL控制方程与材料参数的处理技巧2.1 潜热处理等效热容法的思想与阈值COMSOL 里处理潜热最常用也最稳妥的是等效热容法也叫表观热容法。思路很简单把潜热 ΔH 摊到一个很窄的相变温度区间 [Ts, Tl] 内在这个区间内的等效比热写成cp_eff cp ΔH / (Tl - Ts)从公式能直接看出为什么“相变区间”这个参数这么敏感。用铝来算一笔账固相线 655°C液相线 660°C熔融潜热约 397 kJ/kg区间宽度 5 K那等效比热就是900 397000 / 5 ≈ 80,300 J/(kg·K)铝的基础比热才 900 J/(kg·K)相变区间的等效比热直接变成 89 倍。你会看到材料属性在相变点附近出现一个极其尖利的峰——这个峰不是物理上不存在而是数值上极度难收敛。如果区间缩到 2 K峰值会到 199,400 J/(kg·K)基础值的 220 多倍。所以我在实际项目里的经验是金属凝固的相变区间通常至少给到 3-5 K水/冰这类纯物质可以放在 0.5-1 K但在 COMSOL 里建议也在 2 K 左右用足够平滑的过渡函数去削弱尖峰。精度损失不大但收敛性会显著好转。在 COMSOL 材料定义里这个表达式可以写成cp_eff C_p0 L_ph / (T_l - T_s) * flc2hs((T - T_mid) / (T_l - T_s), 0.5)这里的 flc2hs 是 COMSOL 内置的连续平滑 Heaviside 函数用来把阶跃式的相变“抹”成连续过渡。这也是实现等效热容法最省事的手段。2.2 流动方程怎么选层流、Boussinesq 还是全可压缩流场部分大多数凝固场景都可以用层流假设。凝固过程的液态流速极低特别是在糊状区里流速基本在毫米每秒以下雷诺数很小一上来就上湍流模型属于给自己找罪受。浮力项的处理要留意。COMSOL 的层流接口默认不可压缩流动直接把密度当常数浮力当然也就被忽略了。要计入自然对流常见做法有两种第一种是 Boussinesq 近似在动量方程的体积力上添加F ρ_ref * β * (T - T_ref) * g这个近似假设密度只在重力项里随温度变化其余地方都当常数。对温差不大、密度变化幅度小的过程很有效水从室温到接近冰点差十几度用 Boussinesq 没问题。第二种是使用“非等温流”耦合节点把层流和流体传热真正联立起来再勾选考虑密度变化。凝固过程其实有固体和液体两相密度差异往往不能忽略水变冰密度从 1000 掉到 917 kg/m³铝液凝固后也会收缩。如果流动和界面形态是研究重点建议直接用 COMSOL 的“非等温流”多物理场节点它会自动把传热对流体属性的影响带进去。2.3 溶质场与偏析什么时候非加不可溶质再分配是凝固仿真里最容易被忽略、但工程影响极大的一块。如果只算温度场和流场你只能回答“哪里先凝固”回答不了“哪里成分偏析了”。宏观偏析的经典描述是 Scheil-Gulliver 方程Cs k0 * C0 * (1 - fs)^(k0 - 1)其中 k0 是平衡分配系数。这个式子表达了一个物理事实先凝固的固体它的溶质含量和母液不一样剩下的液相会被不断富集凝固到最后的地方成分往往最“偏”。在 COMSOL 里实现宏观溶质输运要添加“稀物质传递”接口并额外加一个与固相分数变化率成正比的源项模拟凝固过程中液相浓度被排斥富集的过程。很多教材不会讲这个源项怎么加我的经验是先在传热接口里用变量方式导出固相分数 fs然后在稀物质传递的“反应”节点里写R (k0 - 1) * c * d(fs, t) / (1 - fs eps)配合合适的初始浓度和扩散系数就能在宏观尺度上看到通道偏析这类现象的雏形。但老实说这种连续介质方法对糊状区内的枝晶间偏析并不够精确。真要做微观枝晶生长得换相场法。工程上评估“哪里容易偏析”用稀物质传递加 Scheil 级别的源项已经能给出很好的倾向性判断。2.4 材料参数的常见翻车操作做 COMSOL 仿真的朋友应该都知道参数翻车导致的发散和假结果比求解器问题还常见。我自己踩过和见过别人踩的典型有这么几种用室温固态的物性替代液相物性。最典型的就是导热系数和粘度铝固态导热系数 237 W/(m·K)液态只有 100 左右差了一倍多硬用固态值算出来的温度场必然偏。潜热单位搞错。COMSOL 里材料属性是严格按单位处理的如果从论文里看到 J/mol没有换算成 J/kg 就填进去潜热直接放大几十倍结果自然是温度平台长时间不动。液相线和固相线填反。看起来很低级但我在多组分合金里确实犯过一不留神把 Ts 填成了 Tl等效热容被分配到错误的温度区域整个凝固顺序完全乱掉。相变区间设太窄。这个问题在 2.1 已经说过我这里再强调一次纯理论区间 1 K 看着严谨数值上基本是给自己挖坑。工程仿真不是一个纯粹的“精确”问题模型的数值可求解性也是方案的一部分。忽略接触热阻。如果是砂型铸造或者有模具铸件和模具之间不是理想接触传热系数可能只有几百 W/(m²·K)而不是理想接触下的“连续温度”。这个边界条件怎么设对凝固时间的影响可能比材料导热还大。3. 从几何到求解器搭建凝固模型的实际操作流程3.1 几何与网格先根据热扩散尺度估算搭模型之前先别急着画图而是估算一下这个问题的热扩散特征长度。扩散尺度公式很简单L ≈ sqrt(α * t)其中 α 是热扩散率t 是你关心的物理时间。以冰蓄冷为例水的热扩散率约 1.34×10⁻⁷ m²/s取 1 小时的凝固时间扩散长度大概是 sqrt(1.34e-7 × 3600) ≈ 0.022 m也就是 22 毫米。那你在冷壁面附近、凝固前沿会经过的区域网格就不能比这个尺度相差太大否则界面推进一次就跨过好几个单元温度场和固相分数的分辨率完全不够。我的经验是在扩散长度内至少布置 10 到 20 个单元也就是说这个尺度下界面附近网格要 1-2 mm 量级。而一个 200 mm 的铸铝件铝液热扩散率约 3.56×10⁻⁵ m²/s冷却 60 秒扩散长度也接近 46 mm。所以粗算下来界面附近网格 3-6 mm 才说得过去其余区域可以放宽到 10 mm 以上再用边界层网格照顾壁面。几何维度的选择也很关键。能对称就对称能用二维或轴对称就不要一上来开三维全尺寸模型。很多圆坯、板材问题用二维轴对称模型就能拿到足够工程精度的结果计算量却小一两个数量级。先二维跑通物理过程再决定要不要升级到三维这个顺序永远不会浪费。3.2 多物理场节点的装配流体传热、层流与溶质的联动COMSOL 里做凝固仿真我的默认起点是“流体传热”加“层流”。在物理场接口里分别添加之后直接在“多物理场”节点下自动生成“非等温流”耦合这一步能省很多事它会正确地处理流速对温度方程的对流输运以及温度对流体属性的影响。需要提醒的是传热接口里的对流项不是自动全部打开的。“流体传热”在纯固体区域会自动退化掉对流在流体区域才有 u·∇T 项。如果你的模型里有一个区域是固体金属它就不应该有流场所以你需要用一个“指定位移/变形”或干脆把流场速度在固相区强制归零。最省事的做法是我稍后要讲的“糊状区阻力法”用一个随固相分数变化的阻力项让速度在凝固后自动变成零。溶质场的接入在物理场列表里选“稀物质传递”并指定扩散系数和对流速度来自层流。这样三套方程就联动起来了传热决定固相分数固相分数通过阻力项和源项分别作用回动量和溶质方程。3.3 凝固前沿处理动网格的边界感和“糊状区阻力法”很多人想到追踪凝固界面第一反应是 COMSOL 的“移动网格”或“变形几何”。我的建议是宏观大尺度凝固问题不要用移动网格除非你对界面几何有极其明确的定义和把握。原因很简单真实凝固不是一个尖锐的几何边界而是一个宽度可观的模糊区。糊状区里液固共存流动阻力连续变化。硬要用动网格去追踪一条“边界”不仅几何不好定义网格变形到一定程度还会翻折导致雅可比行列式变成负值整个求解直接崩掉。我在铸造仿真里最常用的替代方案是湖状区阻力法F_drag -A_mushy * (1 - f_l)² / (f_l³ eps) * v其中 f_l 是液相分数A_mushy 是一个人工阻力系数典型的量级在 10⁵ 到 10⁸ kg/(m³·s)。这个式子表达的是当液相分数接近 1 时阻力接近 0速度自由发展当凝固程度升高、液相分数趋近 0 时阻力趋于无穷大速度被“冻住”。我在层流动量方程的源项里加这个表达式让流动区在凝固前沿自然消失完全不用动网格数值稳定性也远好于几何边界追踪。如果想要更连续还可以用 flc2hs 平滑这个过渡让阻力变化不会陡成一面墙。3.4 瞬态研究设置BDF、全耦合与步长上限求解器设置是另一个大坑我见过太多人拿着默认设置直接跑相变问题然后对着发散报错发呆。下面这张表是我经过多次实际项目后一套比较稳的默认起点按网格规模可以微调设置项推荐值原因时间步进BDF最大阶 2-5对非线性瞬态问题稳定适合相变这种强非线性初始步长自动或 1e-3 s防止一开局就跨过相变尖峰最大时间步长相变区间宽度 / 最大温变速率保证每个相变平台里有足够采样点全耦合阻尼0.7-0.8压低迭代振荡收敛慢一点但稳非线性求解器最大迭代数25默认可能偏小复杂耦合下不够用线性求解器PARDISO 或 MUMPS直接求解器在强非线性下更稳内存允许就优先最大时间步长这条值得单独解释如果你设置了 5 K 的相变区间而局部冷却速率最猛的地方有 2 K/s那么一个时间步最多走 2.5 秒就会跨完整个相变峰。为了保证在等效比热峰值附近有足够多迭代采样点最大步长最好控制在 0.5-1 s。这一步是在“精度”和“不收敛”之间做的最有效平衡之一。全耦合阻尼也是一个关键旋钮。COMSOL 默认的全耦合阻尼比较激进相变这种强非线性场景很容易振荡。我把阻尼降到 0.7 之后很多原来看起来不可能收敛的模型都安静下来了。4. 仿真发散排查从异常温度到残差暴走的完整定位链4.1 发散的症状残差、负温度和 NaN 的解读先说说发散之前你会看到什么。我总结过几个典型症状残差曲线反复振荡下不去。这说明非线性迭代在有效热容尖峰附近来回弹跳像过阻尼弹簧停不下来。常见原因就是相变区间太窄、时间步太大或者全耦合阻尼太小。温度出现负几百摄氏度。抱错发生在传热方程里负温度就是迭代过程某一步跨过了极端值然后被“放大”到全解域。很多时候是因为潜热等效比热突然变成基础比热的几十上百倍导致线性化失败。NaN 突然出现。这是最彻底的崩溃矩阵已经出现了无穷大或 0/0。常见原因湖状区阻力法的分母 (f_l³ eps) 里 eps 设得太小当 f_l 0 时分母趋近 0阻力项直接爆掉或者稀物质传递的浓度在迭代中变成负值扩散方程取了对数之类运算。报错信息本身往往只有一句“Failed to find consistent solution”或“Divergence detected in nonlinear solver”价值不大真正有用的是看“哪一个物理场先起药效”。4.2 我的排查顺序按“物理-网格-数值”三层推进遇到发散不要慌也不要瞎调。我固定的排查顺序是三层第一步砍物理场。先把流场关掉只做纯导热加等效热容看看这个最简模型能不能跑通。如果最简模型都不稳那问题大概率在材料参数或网格不在耦合。如果最简模型稳了再把层流加回去问题大概率出在流动与传热的耦合步进上。第二步观察相变区和时间步。检查等效热容曲线是否平滑把相变区间从 2 K 扩到 5 K 再试。同时手动限制最大时间步长到 0.1 s先跑过最开始几秒物理时间把初始瞬态尖峰消化掉再放开步长。第三步调求解器。全耦合阻尼降到 0.5最大迭代数调 25不行就换成分离式求解器把“流体传热”和“稀物质传递”分开求解。面这种多物理场强耦合问题分离式有时反而比全耦合更容易收敛只是每个物理场之间的信息交换会滞后一步。我做过一个小统计身边的同事遇上的凝固仿真发散大约 70% 是相变区间过窄或时间步过大20% 是材料参数单位搞错只有 10% 才轮得到真正的求解器配置问题。所以排查顺序一定是先物理再数值别一上来就去动求解器。4.3 真正收敛的标准能量守恒、网格无关性和物理指标“软件没报错”不等于“收敛了”。这是很多初学者最大的误区。我自己判断一个凝固仿真是否真正收敛至少看三样东西。第一是能量守恒。COMSOL 后处理里可以对整个域做热量积分的验证边界流入的总热量加上初始内能应该等于最终内能加上累计相变热。如果这账对不上五个百分点以上结果基本不能信。第二是网格无关性。把界面附近网格加密 1.4 倍把最大网格尺寸从 2 mm 缩到 1.4 mm再看关键指标——比如某个节点温度达到固相线的时间——变化是否小于 1%。如果变化很大说明网格还没足够细继续加密。第三是物理合理性。凝固时间如果比理论估计或实验值差了数量级那不管残差曲线多漂亮都是白搭。我一直强调仿真器的友好界面会给人“算出来就是对的”的错觉但真实工程世界里判断结果合理性的能力才是核心技能。5. 后处理里的工程信号固相分数、温度回升与缺陷倾向5.1 固相分数判断凝固进度比温度更可靠很多人看云图习惯只看温度但在凝固仿真里真正有意义的是固相分数 fs。温度只能告诉你“有没有到相变温度”而 fs 告诉你“这个位置已经完全凝固还是仍处于糊状区”。在 COMSOL 里我通常定义一个变量fs flc2hs((Tl - T) / (Tl - Ts), 0.5)然后在后处理里画 fs 的等值线或者做全域平均。当这个平均值等于 1 时代表区域全部凝固完成。这个“全部凝固时间”比“某个点降到多少度”更能代表一个铸件的整体凝固节奏。判断一个铸件哪里是最后凝固位置就画 fs 最小的区域往往就是热节的位置。热节处最容易形成缩孔缩松如果不做补缩措施问题大概率会出现在那儿。5.2 潜热释放与温度回升别被大时间步骗过去一个很有意思的物理现象是温度回升。纯物质凝固时当局部形核开始释放潜热如果散热速度跟不上潜热释放速度局部温度不仅不下降反而会往上回升一小截。在温度曲线上这个区域会呈现一个平台甚至一个微小的驼峰。我见过很多仿真报告里的降温曲线是一条特别光滑的直线完全没有这个平台那基本可以断定时间步设得太大把温度回升抹掉了。这个平台的意义在于它代表了潜热释放与外部散热的拉锯过程抹掉它凝固时间的计算就有偏差。我的习惯是额外输出 dT/dt 的曲线。温度回升在 dT/dt 上表现为一个尖锐的回峰一眼就能看出模型是否捕获了这个关键物理细节。5.3 缩孔、冷隔与偏析在仿真结果里怎么预告仿真做出来之后怎么把这些场变量翻译成“这里可能会出缺陷”缩孔倾向看压力场和最后凝固区域。如果某块区域液相分数最后才到 1同时该处压力下降、流动补缩路径被切断那么缩孔缩松的风险就非常高。我常用“最后凝固时间等值线图”加“压力最小值位置”叠加来判断补缩通道是否畅通。冷隔倾向看两股凝固前沿的汇合状态。铝合金压铸里经常出现两股低温熔体相遇温度已经降到液相线以下界面没能融合形成氧化皮般的冷隔。在仿真里我会看两个熔体前沿到达同一位置时各自的温度如果都低于液相线冷隔风险就大了。这个不能只看最终温度场要看不同时刻的前沿位置。宏观偏析则看溶质场。稀物质传递的浓度云图会出现明显的高浓度区就是被富集液相冲刷过的通道。这类结果在宏观尺度上已经能指导浇注系统和冒口位置的调整。6. 结晶凝固仿真的迁移能力从金属铸件到电池热管理6.1 电池热管理中的 PCM等效热容法的同款迁移做完金属凝固模型后你会发现这套打法几乎可以直接迁移到电池热管理里的相变材料PCM仿真。石蜡这类 PCM 的潜热大约 180-220 kJ/kg比金属低一些但相变区间往往比金属宽所以等效热容峰不会那么尖数值稳定性反而好。锂电池散热中常见的设计是在电池模组间隙填充石蜡/石墨复合 PCM靠相变吸收脉冲热量。我在 COMSOL 里做的第一个电池热管理仿真用的就是和铸件凝固几乎一模一样的代码结构有效热容法加自然对流只是相变点换成了 28-42°C 的温区。区别在于PCM 凝固过程中会有体积膨胀密度变化更大这时候 Boussinesq 近似可能需要换成弱可压缩流动另外 PCM 的导热系数往往很低只有 0.2-0.3 W/(m·K)所以要加高导热填料来补材料的等效导热系数必须重测不能指望纯石蜡的数据库值有多靠谱。6.2 焊接与增材制造里的极端凝固条件如果你想碰更极端的场景焊接熔池和增材制造熔池是同一个体系里的“偏科生”。熔池尺寸小冷却速度高达 10³-10⁶ K/s温度梯度极大温度回升现象可能非常短暂如果时间步不够小完全看不见。这里有两个和前面完全不同的重点。一是热源模型高斯热源或双椭球热源需要写成随时间和位置变化的体热源函数二是熔池内的对流驱动机制不只是浮力还有表面张力梯度驱动的马兰戈尼对流后者在激光焊接里经常占据主导。对于这类问题我的建议是千万别一上来就开全物理场。先做纯导热等效热容模型跑出熔池形状再逐步加上流场和表面张力每加一个物理场都重新验证一遍温度场。这个循序渐进的思路在凝固这类强非线性问题上比任何高级求解器设置都管用。6.3 案例库拆解与个人学习路径COMSOL 自带的案例库和官网案例下载区其实是一个很被低估的学习资源。与凝固、相变相关的最直接案例有相变材料储热、材料熔化、铸造凝固等你不需要照着抄重点拆解它的“多物理场”节点是怎么装配的、材料属性里相变是怎么定义的、求解器的阻尼和步长是怎么设的。我每次接触一个新物理现象第一件事都是去案例库找最接近的那个模型把它从头到尾拆一遍尤其是看“变量”和“表达式”页签里的写法往往能发现很多官方文档里没写透的细节。一个比较顺畅的学习路径是先用流体传热接口把等效热容法的水结冰模型跑通然后加层流做自然对流再引入稀物质传递做溶质输运最后加固体力学看热应力。每一步只改动一个物理场其余保持不动。这样一旦出现发散你可以立即定位到是新引入的那个物理场带来的问题。现在我做任何凝固或相变课题第一步永远不会变先跑纯导热的等效热容模型拿到固相分数随时间变化这条曲线再往里面加流动、溶质和应力。这套流程看着不够炫但它是所有复杂多物理场模型的地基。地基打稳上面的大楼才不至于在第一次迭代时就塌成一片残差曲线。
