做相变仿真这行久了你会发现最绕不开的数学模型就是焓法。不管你是做金属凝固、石蜡储能、还是电池热管理里的熔融相变表面上看问题是“材料随着温度变了物理状态”但真正落到方程里难点全部集中在相变界面的处理上——界面在移动、潜热在释放、流场跟着温度走如果直接用传统的比热容法去硬解数值上有一万种方式让你怀疑人生。这篇文章我就把自己在多个项目里用焓法耦合温度场和流场做相变模型的核心经验拆开讲从模型思路、数学表达、实操设置到问题排查尽量把每一步为什么这么做、有哪些坑、参数怎么选说透。适合正在被相变收敛问题折磨的仿真工程师也适合刚入门想做材料相变过程模拟的同学参考。1. 相变模型为什么绕不开焓法先搞清楚问题本质1.1 相变问题的难点在于移动边界先说个最简单的物理图景。一块材料放在热源旁边温度升高到熔点固体开始变成液体。这个过程里固液界面是不断移动的——外圈先熔化熔化的液体如果受热不均还会流动流动又反过来参与传热界面就变成不规则的、随时间变化的曲面。你让仿真软件去追踪这个界面本质上是在解一个移动边界问题。移动边界问题最让人头疼的地方在于界面位置不是事先知道的结果而是求解过程中由能量平衡自动决定的。如果你用界面追踪法比如水平集法在相变里的一些应用每一步都要显式地重构界面、设置边界条件、更新网格复杂几何下很容易把网格扭成一团而且两相密度差明显时界面附近的质量和能量守恒很难同时保证。焓法不太一样。它不直接追踪界面而是把一个称为“总焓”的物理量作为基本未知量。固相、液相、两相区统一用一套控制方程描述相变潜热被吸收进焓值的变化里界面位置在求解之后从焓值反推出来不需要单独处理。1.2 焓法和等效热容法的本质区别很多初学者一开始接触的是等效热容法也叫表观热容法。思路是把相变潜热折算到相变温度区间内的比热容上比如材料在固相时比热是 2000 J/(kg·K)在液相时是 2100 J/(kg·K)相变潜热是 200 kJ/kg相变区间是 40℃ 到 60℃那就把这段 20℃ 区间内的比热加成(200000 / 20) 10000 J/(kg·K)等效热容一下子翻了五倍。这个方法好用是好用但它有个致命弱点等效热容在相变区间形成了很高的尖峰如果网格不够密、时间步长不够小温度场在穿越相变区间时会发生振荡甚至出现温度来回抖动但不收敛的情况。而且尖峰区间内的积分精度对网格极其敏感你在 A 方案里算得挺好把网格一加密等效热容的积分反而出偏差了。焓法不这样做。它的控制方程里直接把焓 h 作为待求变量潜热以焓的连续变化形式存在不存在高耸的尖峰。从数值稳定性上讲焓法的非线性主要体现在焓到温度的反解上处理起来要温和得多。我个人的建议很直接只要你的模型里有明显的相变现象且相变区间相对较小直接用焓法。等效热容法只适合粗略估算稳态温度分布的情况动态过程和时间相关的相变比例焓法的结果可靠得多。2. 焓法数学模型搭建几个关键表达式必须吃透2.1 液相分数的定义与插值方式焓法的核心变量是总焓 h它的物理定义是显热和潜热之和。对于纯物质或者窄相变区间材料通常把焓写成h ∫ c_p dT f_L · L其中 c_p 是比热容L 是单位质量潜热f_L 是液相分数。f_L 的取值从 0全固相到 1全液相连续变化。这个 f_L 是整个模型里的灵魂变量温度场、流场、材料参数都要挂在它身上。对于纯物质理想情况相变发生在单一温度点 T_mf_L 是一个阶跃函数。但数值计算里没法处理无穷陡的阶跃必须给它一个小的过渡区间。风扇常见做法是定义固相线温度 T_s 和液相线温度 T_l然后在线性插值假设下f_L (T - T_s) / (T_l - T_s)当 T_s ≤ T ≤ T_lf_L 0当 T T_sf_L 1当 T T_l这是工程中最常见的线性插值。如果相变不是线性的也可以用 smoother 的 S 形函数比如基于误差函数或者正弦的函数补插但对于大多数工程模拟线性插值已经能给出足够准确的结果而且数值稳定性更好。这里有个细节值得注意相变区间的宽度不是随便设的。区间设置太窄比如 T_l - T_s 只有 1℃焓法模型里热容变化剧烈求解器每一步都要处理强非线性极易不收敛。区间设置太宽相变过程被人为拉长潜热释放被抹平均摊结果温度曲线失真。通常建议把区间宽度设置为材料实际相变范围的 1.5 到 2 倍或者至少大于二维/三维网格尺寸对应的最小特征时间内的温度增量。2.2 焓值反解温度的逻辑焓法里第二个关键点是如何从焓值反推出温度。因为控制方程解出来的是焓 h但材料参数、换热系数、浮力项都需要温度 T所以必须建立 h 到 T 的逆关系。在分段线性假设下反解逻辑是这样的如果 h h_s ∫_0^T_s c_p,s dT说明处于全固相T h / c_p,s如果 h h_l h_s L ∫_T_s^T_l c_p,l dT说明处于全液相T T_l (h - h_l) / c_p,l如果 h_s ≤ h ≤ h_l说明处于两相区T T_s (h - h_s) / (h_l - h_s) · (T_l - T_s)这个分段逻辑在软件里用条件表达式写直接用 h 作为自变量。实现时的注意点当材料在相变区间内时它的实际比热不是恒定的而是等效变化的。建议整个相变区间的焓差显热加潜热用一个统一表达式计算避免分段点处的不连续导致残差跳变。2.3 参数实例给一个可以抄作业的案例我在一个储能项目里做过典型的石蜡相变模型参数直接拿来举例点地数据做了一定归一化工程上符合典型石蜡特性参数数值说明固相密度860 kg/m³固液相密度差不大时可近似常数固相比热2100 J/(kg·K)液相密度760 kg/m³如果考虑自然对流浮力这个差要体现液相比热2400 J/(kg·K)导热系数0.25 W/(m·K)固液导热接近暂不分段潜热180 kJ/kg典型石蜡级入能量密度固相线 T_s45℃液相线 T_l55℃实际石蜡熔化区间比这更宽我这里做的简化热膨胀系数0.001 /K用于自然对流浮力项动力粘度0.005 Pa·s液态石蜡典型粘度随温度变化略照这个参数你算一下相变区间的焓差h_s 2100 × 45 94500 J/kgh_l 94500 180000 2400 × (55 - 45) 94500 180000 24000 298500 J/kg也就是在 45℃ 到 55℃ 之间每公斤石蜡要吸收大约 298.5 kJ 的能量其中潜热占了 180 kJ比例非常高。这个数字直接决定了储能时间和加热速率你在后处理时看到某个区域内温度始终停在 45-55℃ 区间就是在相变储能温度曲线出现一个明显的平台段。有了这个参数表接下来搭建多物理场耦合就有了抓手。3. 温度场和流场的耦合方式从单向到双向的升级3.1 能量方程里的对流项别漏掉单纯导热模型只解能量方程温度场和流场互不影响。但实际相变过程里液相一旦出现自然对流就会显著改变热量输运。这时候能量方程需要写成包含对流项的形式ρ c_p ∂T/∂t ρ c_p u · ∇T ∇·(k ∇T) Q_latent其中 Q_latent 是潜热源项在焓法框架里可以做到温度方程里。也可以用焓形式直接写ρ ∂h/∂t ρ c_p u · ∇T ∇·(k ∇T)注意对流项里用的是 c_p 乘温度梯度不是焓梯度这在数学上更符合能量守恒的物理分解也是焓法实现时最容易出错的表达式。我在实操中的一个体会是如果材料固相和液相的比热差别较大对流项里的 c_p 到底取固相值还是液相值会导致明显的温差。通常做法是使用与温度相关的分段 c_p(T)借助液相分数对固液两相的值做插值c_p,eff c_p,s (1 - f_L) c_p,l f_L这个表达式看着简单但它保证了能量守恒在固液界面处的一致性。直接把 c_p 设成常数会高估或低估融化速度跑出来的相变时间差个 10%-20% 很正常。3.2 自然对流的浮力驱动Boussinesq 近似液态相变材料里的流动绝大多数是自然对流由温度引起的密度差驱动。如果你直接写流体密度随温度变化的全可压缩 N-S 方程求解器会处理大量和相变无关的声学尺度物理又慢又不稳定。工程上最常用的处理是 Boussinesq 近似。思路是在动量方程里只保留浮力项的密度变化把动量方程写成ρ0 ∂u/∂t ρ0 (u · ∇)u -∇p μ ∇²u ρ0 β (T - T_ref) g其中 ρ0 是参考密度β 是热膨胀系数T_ref 通常取初始温度或平均温度g 是重力矢量方向。浮力项只在液相区有物理意义固相区因为速度被阻尼冻结该项自然不起作用。这里有一个先决条件密度差不能太大一般要求 ρ_l / ρ_s 在 0.9-1.1 之间。石蜡这种材料固液相密度差在 10%-15%Boussinesq 近似还能凑合用如果要更严格就得走变密度流模型。金属相变时固液相密度差也不大同样适用。3.3 糊状区阻尼保证固相区速度为 0 的物理手段这是相变耦合模型里最容易忽略、最影响收敛的关键项。试想液相分数 f_L 0 的区域材料还是固体不应该有任何流动速度。但如果只靠前面那些方程速度场会把整个计算域都当流体去算固相区会出现泄漏速度物理上完全错误。解决思路是在动量方程里加一个“糊状区阻力项”模拟固相区对流动的巨大阻力。最常用的是基于 Darcy 定律的 Carman-Kozeny 表达式F -C (1 - f_L)² / (f_L³ ε) · u这里的 C 是糊状区常数工程上取 10^4 到 10^7 甚至更大 ε 是个很小的正数防止分母为零通常取 10^-3 量级。当 f_L 0全固相时这个阻力项趋向无穷大物理上被 C 限制在一个超大值速度被强推为 0。当 f_L 1全液相时阻力项为 0不影响流动。有个实用经验C 值不能无条件取最大。C 取 10^7 时动量方程里这个阻力项的数量级非常大会破坏矩阵的条件数导致非线性求解器每步都很难收敛。我的习惯是先从 10^4 开始试算确认固相区的速度泄漏量级已经低于 10^-6 m/s 之后再逐步加加到 10^6 为止。大部分情况下 10^5 到 10^6 足够。这个阻尼项还顺便解决了相变界面的数值稳定性问题。没有它液相分数在 0-1 之间的界面区域会产生虚假的对流严重时会在界面两侧形成振荡的温度涟漪画面非常难看。3.4 耦合方式是单向还是双向什么时候必须双向如果你的相变材料是纯导热主导比如热源功率很小、液相区流动很弱那可以只做单向耦合先解温度场根据温度分布算出液相分数再把液相分数映射到材料参数粘度、热容、导热系数更新流场速度。这种解法和解耦求解速度快很多。但如果自然对流的 Rayleigh 数高热源功率大流场反过来显著影响温度分布尤其你想准确预测“融化前锋的形状”时单向耦合完全不够。必须做双向耦合也就是每个时间步内同步求解动量方程、连续性方程、能量方程和焓反解所有变量联动更新。焓法的非线性本身已经够强双向耦合会让非线性进一步增强这时候求解策略就成了关键。我个人的取舍标准是看无量纲数 Ra 是否超过临界值Ra g β (T_hot - T_cold) H³ / (ν α)H 是特征尺度ν 是运动粘度α 是热扩散率。Ra 超过 10^5 时自然对流已经是主要传热机制必须双向耦合。低于 10^4 时单向耦合误差可以接受。这个判断标准我用了很多年基本靠谱。4. 实操流程从几何搭建到后处理的完整执行方案4.1 物理场选择与接口组合我不强调具体软件品牌工程圈都知道这类多物理场工具的通用流程但无论你用哪套平台典型的相变-流-热耦合模型物理场组合都是流体流动接口求解 N-S 方程 连续性方程固体传热/流体传热接口用于能量方程自定义变量/方程用于计算液相分数 f_L、焓 h、温度反解 T(h)这里有个自由度的选择问题。用哪个温度变量去求解决定了表达式写法。我用过两种方案。方案 A主变量直接用温度 T焓法以“表观热容”加“潜热源项”的形式实现。具体做法是能量方程里把热容替换成热焓对时间的偏导数项同时在瞬间添加一个潜热释放源项。优点是与传热接口兼容度最高后处理直接看温度。缺点是实质上回到等效热容路线对网格敏感。方案 B主变量改成焓 h 或使用热容形式材料但这种写法通常要到更高阶的自定义弱形式才能实现通用接口下反而难操作。对于绝大多数工程场景我建议你使用表观热容形式的写法但要把热容表达式精心构造避免直接用尖峰函数。具体做法定义等效热容 C_eff(T)在整个温区内连续C_eff c_p,s (1 - f_L) c_p,l f_L L / (T_l - T_s) · df_L/dT由于 f_L 是线性插值时df_L/dT 在两个相变边界点是阶跃的。为了让数值稳定我给 df_L/dT 用一个小平滑实际效果比理想阶跃好得多。4.2 几何建模与网格分区建立二维或轴对称模型时我通常把计算域画成矩形或圆柱腔体。关键网格处理原则相变前沿可能经过的区域加密加密区的网格尺寸不超过总特征长度的 1/30壁面处添加边界层网格至少 5 层第一层厚度根据预期对流强度设定确保 y 满足低 Re 涡模拟要求如果使用高度函数初始化液相分数 0初始时刻全固相逐步融化焓法和 Level Set 不同不需要界面附近高精度贴体网格因为界面是隐式捕捉的只需要确保相变区间温度变化对应的焓值能被网格积分覆盖。一个可量化的经验在最大局部热流密度下网格间距需要满足使得一个时间步内温度变化小于 0.5×(T_l - T_s)。4.3 边界条件设置实例以一个带内热源的方形储能单元为例左侧壁面加热温度恒定 80℃右侧壁面自然对流散热对流系数 10 W/(m²·K)环境温度 25℃上下壁面绝热初始温度 25℃全固相速度场初始为 0重力方向沿-y轴浮力驱动液相向上流动这些条件写起来都不复杂但耦合设置里有两个坑第一个坑是压力参考点必须设置。纯自然对流问题速度场是由浮力驱动的压力没有唯一绝对参考不设参考点的话压力方程奇异求解器会报错或长时间不收敛。我在所有涉及自然对流的模型里计算域角落固定一个压力零点。第二个坑是壁面附近的浮力方向。如果加热壁在左侧液体受热后沿左壁上升在顶部转向右在右壁冷却下降形成一个逆时针的大涡环。你不能把所有壁面都设为无滑移就完事还要检查浮力项的方向设定是否正确。我第一次做的时候把重力项写错方向算出来的涡环方向完全反了温度场倒是没太大差异但流场后处理一眼就看出来有问题。4.4 瞬态求解器参数与时间步策略瞬态求解设置里最影响结果是时间步进方式用 BDF 2二阶向后差分或等价的隐式时间格式避免显示格式的时间步限制初始时间步设得非常小取 0.001 秒量级让求解器从全固相平稳启动开启自适应时间步进步长上限控制在 1 秒左右下限 1e-5 秒非线性残差容差设置为默认值的 1/10比如 0.001保证界面推移的精度焓法模型的非线性非常强尤其是温度穿越相变区间的瞬间物理特性突变。这种情况下直接全体方程全耦合求解求解器每步的雅可比矩阵更新成本很高但值得。我常用的策略是“先导热后对流”两步走第一步关闭流体求解只算纯导热的相变过程 0.5 秒得到初始温度场第二步再开启流体耦合求解。这个方法能快速跨过最困难的启动段节省大量时间并且显著提高稳定性。4.5 后处理重点指标与可视化模型算完后处理要看哪些东西第一个是液相分数云图。它是最直观的相变过程体现你能够在不同时间切片上看到固相逐渐缩小的全过程轮廓。可以用 0.5 等值线代表界面位置看它的推进路径是否合理。第二个是温度历史曲线。选取几个监测点分别位于靠近热源、中间区域、远离热源位置画出温度随时间变化。典型的相变曲线应该出现斜率先升、中间变缓相变平台、随后斜率恢复的过程。平台段的长度对应潜热吸收时间。第三个是速度场。看液相区的涡环结构是否形成涡心位置大概在哪个区域最大流速数量级是否符合自然对流的常识通常几毫米每秒到几十毫米每秒。第四个是能量守恒校验。这个方法非常必要把整个计算域的焓随时间的变化率与边界热流积分总量去做对比。我每次都要算这个总量差的曲线如果能量不平衡超过 5%说明模型或者求解参数有问题。曾经有一次我漏掉了相变区间的热量计算能量曲线直接缺口排查了半天才发现是液相分数插值表达式写错了一个边界。5. 常见问题与排查技巧实录5.1 温度场锯齿振荡界面像锯齿怎么办锯齿振荡基本都出在液相分数和潜热的处理上。排查顺序第一步检查相变区间的宽度是不是太小。如果 T_l - T_s 只有 0.5℃网格稍粗就被温度跳跃跨过了必然振荡。把区间放大到 2-5℃如果振荡消失就是这个问题。第二步检查等效热容是不是存在尖峰。理想线性插值产生的 df_L/dT 是常数不会尖峰但如果你用了非线性光滑函数相位变化集中在极窄区间内等效热容就有尖峰隐患。第三步检查网格在相变前沿区域是否足够细。界面处温度梯度大网格太粗时界面推进是一格一格跳的看起来像锯齿。5.2 计算发散提示非线性求解器不收敛最典型的场景发生在材料刚开始熔化那几步。这时候潜热突然释放温度对焓的非线性映射变化剧烈牛顿迭代容易振荡。排查技巧第一调小初始时间步长。把这个参数缩小一个数量级通常就能跨过启动难关。第二给液相分数加“松动”。英语里叫 under-relaxation也就是液分更新不用全量每一步只更新一部分f_L,new f_L_old α (f_L_calc - f_L_old)α 取 0.5-0.8。这个做法牺牲一点点时间精度但换来了极其稳定的收敛过程。第三如果时间步长已经非常小而仍然发散查看是不是糊状区常数 C 太大。把 C 从 10^6 降到 10^4 试试。有的模型里 C 过强会让动量方程局部条件数剧增。第四检查初始条件是否自洽。初始温度恰好卡在相变区间边界上有时候会触发初始残差极大。把初始温度稍微调整 0.1℃或者初始给一个极小的液相分数 1e-4也能绕过。5.3 速度场在固相区出现非零速度如果固相区速度不是零先检查糊状区阻尼项是不是写进了方程。这个大概率是漏写了这个源项。不是漏写的话检查 C 值是不是太小增大一个数量级看看。另外注意 f_L 的数值精度。如果 f_L 在固相区因为数值误差变成 1e-5 而不是精确的 0那个小值被代入阻尼公式1-f_L²/(f_L³ε)虽然分母有 ε 保护但 C 值如果不够大仍旧会出现微小泄漏。5.4 相变时间过短或过长和理论值明显不符相变持续总时长对不上理论估算通常在两个地方第一潜热值是不是被正确“激活”了。有时候材料定义里 L 被填成 0潜热完全不起作用融化时间比理论值少 40% 以上。排查方式是在能量守恒曲线里看潜热贡献段。第二相变区间宽度被设置得过宽。假设实际材料相变区间只有 1℃但你在模型里为了数值稳定把它拉成 20℃潜热在这个过程中均匀释放虽然消耗的总能量一致但热流分配方式和实际情况不一样会导致融化过程持续时间偏长。这里要做一个平衡数值稳定和物理保真之间的折中。5.5 快速排查清单症状首要怀疑点推荐措施温度锯齿振荡相变区间太窄/网格太粗增大 T_l-T_s加密网格不收敛初始步长过大初始步长缩到 1e-4 再试流场冻结/泄漏糊状区阻尼缺失检查 C(1-f_L)²/(f_L³ε)项相变时间偏短潜热没接上检查 L 的表达式链接温度平台不明显等效热容无尖峰但区间宽适当压缩相变区间5.6 一个非常隐蔽的坑热物性参数插值固相和液相的热物性不同尤其导热系数可能差好几倍。如果你只在材料属性里用分段函数切换而没有通过液相分数做光滑过渡界面处的导热系数突变会带来严重的数值不连续性。建议所有依赖于相态的物性参数都统一写成k_eff k_s (1 - f_L) k_l f_L这样的形式。这不是额外的物理假设而是数值处理上的必需操作。我当年做金属相变时固态的导热系数比液态高大约 60%最初没做过渡处理温度场界面处出现了无物理意义的高频扰动。改成这个直接插值后一切都稳了。关于参数选择我最后再啰嗦几句相变仿真做到最后瓶颈往往不在软件操作而在你对材料相变行为的物理理解有多深。焓法这个数学框架本身非常优雅它把一个强非线性移动边界问题转化为单域固定网格问题只要掌控了液相分数、焓反解、糊状区阻尼这三个核心环节绝大多数工程相变场景都能稳定算出来。最后分享一个我个人的实践经验拿到任何相变模型第一件事不是急着调求解器而是先验算焓差和能量总量。在 Excel 里花十分钟把固相升至熔点、潜热吸收、液相继续升温这三段能量相加和总加热功率乘以预期时间做个对比。如果这个粗算和仿真结果的量级对不上后面调什么都白费。把这一步养成习惯之后你做相变模型的速度和一次成功率高得惊人。
