十几年前我刚接触COMSOL时觉得相变模拟只是给传热方程加一个潜热项。后来真正用COMSOL模拟石蜡相变储能单元才发现“加一个潜热项”的背后全是坑。第一次跑熔化模拟眼看温度云图已经非常漂亮液相界面却怎么都不肯移动模型的探测器曲线像心电图一样上下跳。后来一步步排查才明白是液相分数与网格尺度不匹配等效热容的尖峰把时间步长拉爆了。那段时间我几乎把COMSOL的文档、案例库和社区问答翻了个遍慢慢形成了自己的相变模拟方法论。这篇博文不打算讲太多理论推导重点分享这几年来用COMSOL做石蜡、熔盐、金属三类相变材料的实操经验包括材料参数怎么设、移动网格怎么调、单元活化怎么用、收敛失败怎么排查以及热管和电磁耦合等进阶场景的建模思路。无论你是刚下载COMSOL想跑通第一个相变案例还是已经在做激光焊接熔池仿真下面的内容应该能帮你少走不少弯路。1. 相变模拟到底在求什么从“表观热容”到“移动界面”的两条路线1.1 潜热到底怎么进控制方程相变模拟的第一步不是急着打开COMSOL选物理场而是先弄清楚“相变”这两个字在方程里到底意味着什么。常规传热方程只考虑显热也就是材料温度升高所吸收的热量。但发生固液相变时材料在相变点还会额外吸收或释放一份潜热这部分能量不会直接体现为温度变化而是体现在焓的变化上。正确做法是把控制方程写成焓形式对单位质量焓做时间导数再与热传导项平衡。这个念叨起来像教科书的话实际上决定了你在COMSOL里应该使用哪种材料定义方式。COMSOL用户最常接触的“表观热容法”本质上就是把焓方程重新改写。具体来说把潜热折合成一项附加热容叠加在真实热容上于是就有了 C_eff C_p L·df/dT其中 f 是液相分数从0变化到1。这个做法最大的好处是可以直接沿用固体传热接口的所有数值算法不需要额外写质量守恒或动量方程稳定性也更好。可它也有一个固有缺陷潜热并不是在一个温度点集中释放而是被“抹”在一个有限宽度的温度区间里。区间设置得太窄df/dT 会非常大时间步长会被逼到纳秒量级区间设置得太宽熔化过程被严重“稀释”温度平台失真。我自己在COMSOL里定义 f 时最常用的是平滑阶跃函数flc2hs(T - Tm, dT)比如石蜡熔化区间 45~48℃可以取 Tm46.5、dT2。这个函数的好处是导数连续不会像线性分段函数那样在固相线和液相线处造成热容跳变。如果你在材料属性里用的是手动表达式一定要用d(f(T), T)来求导数而不是自己写近似的差分公式否则求解器内部的一致解算会出现偏差。1.2 等效热容法之外的另一条路线移动界面等效热容法听上去很美但它把相变界面模糊化了。对于石蜡、熔盐这类商业材料熔化区间本来就比较宽界面其实是一个“糊状区”用等效热容法处理非常合适。可是金属熔化不一样尤其是激光焊接里的熔池固液界面锐利到只有微米级等效热容法会把界面抹在一个网格单元里熔池边界根本看不准。这时候就需要第二条路线空间显式追踪界面。在COMSOL中显式追踪界面主要有两种实现方式。一种是用“变形几何”接口让网格节点跟随固液界面移动在每个时间步都保持界面为一个清晰的几何边界另一种是引入水平集或相场变量用颜色函数隐式捕捉界面位置。前者适合界面位移不是特别剧烈、拓扑不变的场景后者适合气泡、液滴分裂合并这种大变形问题。选择时不要贪图“看起来高级”要问自己真正关注的是什么如果只要温度场和熔化时间等效热容法够了如果要做熔池形貌和匙孔演化移动网格或水平集逃不掉。方法适用材料COMSOL入口优势风险等效热容法石蜡、熔盐、合金固体传热 相变材料特征实现简单、计算稳定界面被抹平无法追踪变形几何/移动网格纯金属、激光焊接熔池固体传热 变形几何界面清晰可输出边界位移网格大变形易纠缠水平集/相场法蒸发冷凝、枝晶生长流体流动接口 水平集可处理拓扑变化计算量大参数多1.3 建模前先回答的四个问题我在帮别人排查相变模型时发现大量返工模型都是因为一开始没有想清楚四个问题。第一相变温度是单值还是区间纯物质是单值合金和商业石蜡是区间这决定了 f(T) 的表达形式。第二熔化前后密度变化大不大会不会引起自然对流如果液相密度差产生的浮力流不可忽略就必须把层流接口加进来。第三固相和液相的导热系数差多少这个差距影响等效热容法中的导热插值也影响自然对流的强度。第四你是否关心材料内部的应力变形如果要焊接残余应力就不能只做传热必须顺序耦合固体力学。这四步思考看起来简单但能帮你省下大量盲目尝试的时间。比如石蜡储能盒如果只关心总储热量完全可以不开层流用纯导热模型快速估算但如果关心局部热分层、熔化时间那必须开层流非等温流动。再比如激光焊接如果只关心热影响区宽度纯传热模型可以做初步估算但要做熔池形状就得移动网格。建模没有银弹先想清楚输出要求才能决定输入复杂度。2. 石蜡与熔盐模拟实测材料参数、边界条件和那些看不见的坑2.1 石蜡储能单元为什么纯导热会严重低估熔化速度我的第一个相变项目是圆柱形石蜡储能单元外壁温度恒定加热。当时为了省事只在固体传热接口里定义了石蜡的等效热容和热导率没有开流体计算。结果模拟出来的熔化完成时间比实验整整慢了约30%而且最后阶段的温度曲线明显过冷。原因不复杂石蜡熔化后液相区的温度不均匀会产生密度差热流体上升、冷流体下降形成自然对流。这个对流对换热的影响极大液态石蜡的有效传热能力可以比纯分子导热高出数倍不把它算进去熔化速度当然会严重偏低。在COMSOL中把自然对流加进去并不难。建立一个“层流”接口与“固体传热”形成非等温流动耦合然后在体积力节点中添加浮力项-rho0 * g * beta * (T - T_ref)。石蜡的热膨胀系数大约在 0.001/K 量级比水还大浮力项会显著影响压力场所以直接全量加载时压力残差很容易爆掉。一个有效改法是先跑一段纯导热让温度场有个基础分布再把浮力项乘以一个随计算时间从0到1线性增加的斜坡函数让自然对流逐渐参与进来。这样给非线性迭代器一个缓冲模型往往能顺利跨过初始阶段。另外在定义固液区域的流动物理场时不要贪图简单把整个域都设成层流。固相区还是固体不该有流动自由度。工程上常用的做法是给液相定义一个“糊状区孔隙率”当温度低于固相线时Darcy阻尼系数设成很大迫使流速趋于零当温度高于液相线时阻尼系数设成0。这个阻尼系数我曾经从 10^5 调节到 10^7 kg/(m³·s)数值差异非常大需要结合你模拟对象的几何尺寸去试。2.2 熔盐蓄热系统温度相关物性和空心圆柱几何熔盐和石蜡虽然都属于潜热储能材料但熔盐的工作温度高液相温度范围很宽物性随温度的变化也远比石蜡强烈。我们在做熔盐蓄热模拟时最忌讳把黏度、导热系数和比热设成常数。熔盐的动力黏度随温度升高呈指数式下降密度也是如此。如果这些参数用常数值自然对流流场会失真冷热分层的位置、边界换热系数全都会偏。我在COMSOL里推荐的做法是在“材料”节点中把热导率、密度、比热、动力黏度都定义成温度的函数可以直接用实验数据的插值表也可以用拟合公式。比如有些熔盐黏度符合 Andrade 公式eta A * exp(B / T)只需要两个拟合常数。这样简单但有物理依据比自己拍脑袋取值放心得多。熔盐储罐的几何常做成空心圆柱。做这类模型时很多人会纠结在COMSOL里怎么画空心圆柱。最稳妥、也最利于网格收敛的方法不是用三维圆柱差集而是用二维轴对称画一个矩形截面然后利用“旋转”操作生成空心圆柱体。这样网格可以沿径向和轴向自由控制边界层能精确加在内壁和外壁。如果你真的需要在三维里做空心圆柱建议先用草图画一个圆环截面再执行“拉伸”而不是先建实心圆柱再布尔减掉内芯后者经常会产生退化面让边界条件选择出现多余的面编号。2.3 一个可以直接套用的等效热容参数模板下面这套参数定义我在多个项目里用过适合熔化区间较宽的石蜡和熔盐。第一步在“全局定义”里创建参数Tm、dT、L、rho_s、rho_l、k_s、k_l、Cp_s、Cp_l。第二步创建液相分数解析函数f flc2hs(T - Tm, dT)注意这里的 T 是COMSOL中的温度变量你可以用T直接在解析函数表达式里引用。第三步定义等效热容、等效密度和等效导热系数。COMSOL表达式可以写成rho_eff rho_s (rho_l - rho_s) * f cp_eff cp_s (cp_l - cp_s) * f L * d(f, T) k_eff k_s (k_l - k_s) * f第四步把材料属性里的密度、恒压热容和导热系数分别换成上述表达式。注意d(f,T)在COMSOL符号中表示解析函数 f 对温度 T 的导数不要自己在旁边写一个积分或差分版本。如果用内置的“相变材料”特征可以直接填潜热和相变区间内部也是这一套逻辑但使用自定义表达式的好处是你能完全控制 f 的形状甚至在 f 里加入滞后效应。还有一个容易被忽视的问题固相和液相的密度差到底怎么处理。你如果把密度做成 f 的线性插值那在熔化区间内体积会变化热膨胀应力计算就会自动包含体积变化。可如果你同时用不可压缩层流计算自然对流密度差还会产生浮力项两个效果叠加时容易矛盾。我的建议是传热计算中的密度插值只用于焓和热容计算自然对流的浮力项单独用 Boussinesq 近似的参考密度和膨胀系数不要混用同一组密度表达式。3. 金属熔化与激光焊接移动网格、单元活化和凝固收缩3.1 为什么金属相变必须换一套思路金属和石蜡、熔盐最大的不同在于金属熔化温度高、热导率高、固液界面非常锐利而且液态金属的表面张力和 Marangoni 效应会直接影响熔池形状。当你做激光焊接或电子束焊接仿真时光靠等效热容法是不够的。它虽然能给出一个大致温度场但熔池的轮廓、表面凹陷、匙孔演化都需要一个真正能跟踪界面的数值方法。我早期做过一个激光点焊模型只用等效热容法观察熔化深度。结果熔池宽深比明显失真实验焊缝截面呈现“钉头”形状而模拟结果却是半圆形。后来把模型升级为“固体传热 变形几何 层流”再配上热流密度源项熔池形状才和实验焊缝截面基本吻合。这个经历让我意识到金属相变模拟的核心不是“潜热”而“界面”。潜热只决定能量守恒界面形态才决定最终焊接质量。3.2 激光焊接仿真里怎么实现“单元活化”很多人问COMSOL有没有像ANSYS生死单元那样的“单元活化”功能。严格来说COMSOL没有一个独立的一键“单元生死”开关但在激光焊接、增材制造仿真里完全可以实现等效效果而且实现思路比生死单元更细腻。我常用的套路是把所有可能熔化的金属都提前建在一个域里不要做几何布尔减除然后在层流接口中把液相体积分数 alpha_L 定义成温度的函数当alpha_L 0.5时该区域的动量方程被一个巨大的阻尼力冻结流体速度强制为0看起来就像单元被“杀死”了当温度超过液相线后阻尼消失流体自由度重新激活等效于单元被“活化”。具体到COMSOL操作先建立一个“层流”接口然后在流体属性节点中选择“孔隙率”或自定义阻力源项。阻力项可以写成S -K_darcy(T) * u其中 K_darcy 在低温固相区取 10^6 量级高温液相区取 0中间用平滑阶跃函数过渡。这样就不需要真正的网格拓扑变化数值稳定性也好很多。有人觉得这样做不是真正的“单元活化”但工程上关心的不是算法名字而是最终得到的熔池边界。只要液相分数场的梯度足够陡你就能从结果中提取出清晰的熔池轮廓。如果你是冲着COMSOL案例库里“激光焊接单元活化”这个关键词来的我想额外提醒一句链接里很多工业级案例都使用了“网格自适应与重新划分”功能也就是说当移动网格变形达到一定阈值后求解器会暂停自动重画网格再把旧网格上的物理场插值到新网格上。这个功能对捕捉匙孔非常有帮助但也要付出额外代价那就是计算时间成倍增加。从我的经验看如果没有强热源导致自由表面剧烈变形先用固定网格加阻尼法往往就够用了。3.3 凝固收缩、热应力和顺序耦合金属熔化凝固过程中的体积变化会在焊接件里留下残余应力甚至产生凝固裂纹。若要评估这个问题COMSOL里需要把固体力学耦合成一个结构分析。这里有一个必须避开的坑熔化状态下的金属几乎没有剪切刚度若材料弹性模量直接设为0刚度矩阵会奇异求解器直接报错。工程做法是给高温区一个很小的弹性模量比如 1 Pa而不是0这样方程仍然满秩计算能跑下去而应力结果在熔池内部会自动趋近于0。焊接残余应力计算我不建议做全耦合瞬态分析。温度场、流场、结构场三步同时求解会让非线性程度爆炸时间步长会被压得非常小。实践中我更推荐顺序耦合第一步用“固体传热层流”计算整个焊接过程的瞬态温度场并存储足够密的时程结果第二步把温度场作为载荷映射到结构网格上关闭流动只做固体力学求解。这样既可以得到应力演化也避免在一个方程里同时处理熔池流动和固体刚度突变。代价是错过了液态金属流动对最终应力分布的直接影响但在焊缝应力评估的精度范围内这个误差通常可以接受。4. 仿真跑不动/不收敛网格、时间步长和求解器配置的排查清单4.1 网格加密策略先用粗网格找趋势再加密界面相变模拟遇到不收敛时很多人第一反应是减小时间步长但真正问题往往出在网格。等效热容法下潜热只在一个窄温度区间内释放如果这个区间内的网格太大液相分数根本无法形成连续移动的界面。典型表现是温度场已经越过相变点好几次但液相分数仍为0或者熔化界面像拼图一样一块一块地跳跃。你需要在预计的界面移动区域加密网格至少保证界面宽度内有 5~10 层网格。COMSOL的网格细划可以直接针对域或边界设置“最大单元尺寸”。我的原则是先用粗网格跑少量时间确认温度分布趋势和加热功率匹配再用网格细化研究界面区域的局部加密。不建议一开始就用自适应网格和移动网格叠加否则网格变形和自适应误差会混在一起出了问题很难定位。先把固定网格版本调到稳定收敛再引入移动网格。4.2 时间步长和求解器配置最大时间步不是越小越好COMSOL的瞬态求解器默认使用自适应时间步长但“自适应”不代表你可以完全不管。相变模拟时时间步长如果太大求解器会跳过潜热释放峰如果太小计算时间又会爆炸。我习惯在“研究”设置里手动指定最大时间步取特征热扩散时间除以50。特征热扩散时间大约是特征长度的平方除以热扩散系数。对石蜡储能单元特征长度10 cm热扩散系数约 0.1 mm²/s特征时间约 10000 秒最大时间步控制在 200 秒左右就够但激光焊接熔池特征长度只有几毫米特征时间非常短最大时间步常常要设到 10^-4 秒量级。另一个容易被忽略的是“非线性方法”设置。COMSOL默认用Newton迭代和阻尼因子遇到收敛困难时会自动减小步长。但有些情况下你会发现残差卡在某个量级怎么都降不下去尤其是浮力项和Darcy阻尼项同时作用时。这时可以在“全耦合”求解器配置中把“最大迭代次数”从默认的 4 增加到 8并打开“自动补偿”选项往往能帮求解器多挣扎一下。4.3 一次典型不收敛问题的完整排查链路有一段时间我在跑熔盐储热模型压力残差上升到 10^3 量级再也不降。我按这个顺序排查最后终于找到问题。首先检查网格把自然对流边界层加细压力振荡稍好但没有根治。第二步检查物性发现熔盐黏度在低温区指数式增长导致局部区域雷诺数极低压力方程几乎奇异加密网格后反而更差。第三步修改浮力我把 Boussinesq 近似的 beta 从线性差值改成平滑阶跃让浮力在固液界面处不会突变残差立刻下降。第四步调整时间步把最大时间步缩小为原来的1/5让熔化界面在每个步长内只移动不到一个网格单元。最终模型稳定收敛计算结果和实验对上了。这个排查过程想说明一件事相变模拟里的不收敛绝大多数不是单一原因不要指望只改一个参数就能过关。先网格后材料再物理场耦合最后时间步这是最有效的排查顺序。现象可能原因首选排查温度曲线锯齿振荡等效热容尖峰过高增大 dT熔化界面不移动网格太粗加密界面区域压力残差不降Boussinesq 浮力突变平滑 beta液相区速度回旋Darcy 阻尼不合理调整 K_darcy移动网格纠缠网格随动范围过大限制动网格域5. 相变模拟的边界拓展热管仿真、电磁热源耦合和弱形式方程的启示5.1 热管仿真里的蒸发冷凝相变不能照搬等效热容热管仿真看起来很接近相变储能但它的难点在蒸发段和冷凝段的相变并不是简单的固液相变而是液气两相流动和相间传质。等效热容法可以把蒸发吸热作为体积热源加入但无法描述液膜收缩和蒸汽流动整体热阻会严重失真。我建议第一步先做一维或轴对称的简化模型把蒸发端的吸热通量和冷凝端的放热通量用传热边界条件代替校验整体温度分布。第二步再考虑用水平集或相场方法模拟气液界面。千万不要一上来就建一个完整的三维热管模型模型复杂度会迅速失控最后什么也没算出来。COMSOL中做蒸发冷凝还有一个常用技巧在气液界面加载一个质量源项并用一个平滑函数将质量源铺在一个很小的过渡区内而不是集中在边界上。这样做可以避免相变界面处的速度突变引发压力振荡。如果你看过一些COMSOL热管案例会发现很多案例用的其实是“有效导热系数法”把蒸发冷凝巨大换热量折算成等效热导率。这个方法适合工程快速评估适合做系统级仿真但不适合研究热管内部机理。5.2 电磁热源与相变耦合注意趋肤深度感应加热、微波加热这类场景里相变模拟和电磁场仿真需要耦合。COMSOL的“电磁场固体传热”多物理场耦合可以直接把电磁损耗作为热源但很多人忘了检查趋肤深度和网格的关系。交流电磁场在导体中的趋肤深度 δ sqrt(2/(ω μ σ))在趋肤深度层内至少要保证 3 层网格否则热源分布完全失真相变界面位置也会跟着错。比如感应加热一个空心圆柱金属坯料电磁场集中在表面如果网格只在内部加密、表面只有一层粗网格算出来的表面温升会严重偏低。我通常会在表面边界层加 5~8 层网格每层厚度按趋肤深度的 1/5 递减。另外如果金属材料在相变点附近电导率和磁导率发生突变电磁损耗的分布也会突变这时更要缩小时间步长否则热源阶梯变化会导致温度场振荡。5.3 弱形式方程能给相变模拟带来什么看完“基于comsol弱形式方程求解色散光子晶体能带”这个搜索词可能会觉得弱形式离相变模拟很远。但我的体会是理解了弱形式就等于理解了COMSOL所有物理接口背后的“编译器”。相变模拟中经常出现自定义本构模型比如记忆合金的焓-温度滞后曲线或者非傅里叶热传导方程预置接口不一定能直接支持。这时候会写弱形式你就能把方程以最自由的方式嵌入COMSOL。初学者不要拿相变方程来练弱形式应该先用最简单的一维热传导方程验证自己的弱形式推导和COMSOL内置结果是否一致。验证通过了再加相变潜热源项。弱形式的好处是你可以随意把液相分数、移动边界条件写进去它的缺点也很明显调试困难一旦不收敛你无法判断是方程写错还是数值算法问题。所以我的态度是能用预置物理接口解决的优先用预置接口只有当预置接口真的满足不了才上弱形式。5.4 与实验对标不要只盯着云图好不好看相变模拟的最终评判标准永远是实验数据而不是云图的配色。每次跑完模型我至少会做三组对标熔化过程监测点温度曲线、熔化完成时间、以及最终液相分数空间分布。如果模拟结果和实验误差超过15%我第一怀疑材料参数——石蜡的相变潜热、熔盐的黏度、金属的高温导热系数这些参数不同批次可能差很多第二怀疑边界条件——对流换热系数的经验值经常不准最后才怀疑模型维度和数值设置。按这个顺序查基本能定位问题。还有一个容易被忽略的细节实验样品里的热膨胀和容器间隙会造成接触热阻但仿真模型默认是理想接触。这个差异在小尺寸样品中会非常明显。所以我会在实验中额外记录样品表面温度和外壁温度差用这个温差反推接触热阻再代入仿真。这个方法比盲目把换热系数调大要靠谱得多。我有一个习惯每次跑相变模拟前都会先用一个纯热扩散模型跑一遍基础加热过程确认网格和时间步长满足要求再往里面加入潜热、自然对流或移动网格。这一步看起来浪费时间实际上能避免至少一半的收敛问题。相变模拟涉及大量非线性和多物理场耦合任何额外的不稳定因素都可能被放大。如果你正卡在COMSOL相变模拟不收敛或者结果不合理不妨回到起点把物理过程一层一层地加回去真相往往就藏在最后那一步。