焓法相变传热分析

固体熔化或液体凝固时，有个特别的现象：温度不变，但热量还在进出——这就是潜热。普通热传导方程 $\rho c_p \partial T/\partial t = \nabla\cdot(k\nabla T)$ 没法处理这个问题，因为此时 $\partial T/\partial t = 0$ 但热量还在流动，热收支就不平衡了。焓法把总焓 $H$（包括显热和潜热）作为主变量，自然地把潜热纳入能量方程。铸造凝固分析、PCM（相变材料）蓄热设计、电池热管理的相变冷却，都离不开这套理论。

想象一锅铝合金在675℃开始凝固。在凝固过程中，铸型持续从熔液抽走热量，但铸件表面温度"卡"在675℃，直到潜热 $L_f = 397\,\text{kJ/kg}$ 全部释放完毕后，温度才继续下降。标准的 $\rho c_p \partial T/\partial t$ 方程在 $\partial T/\partial t = 0$ 这段时间内会认为热流量为零——但实际上还有 $\rho L_f \partial f_s/\partial t$（凝固速率×潜热）的热量在流动。漏掉这部分能量，凝固过程的温度场完全错误。

用总焓 $H$ 表示的统一形式控制方程：

$H$ 同时包含显热和潜热两部分。引入液相分率 $f_l \in [0,1]$，总焓展开为：

$L_f$ 是熔解潜热（J/kg）。对于合金和树脂的糊状区（Mushy Zone），$f_l$ 在固相线 $T_s$ 和液相线 $T_l$ 之间线性变化（Lever Rule近似）：

纯物质（如纯铝、纯水）在单一熔点 $T_m = T_s = T_l$ 处阶跃，是糊状区宽度为零的极端情况。

糊状区是固液共存的"浆状"区域，是固体树枝晶（Dendrite）在液体基体中生长的过渡带。铸造时可以观察到：切开刚凝固一半的铸件，会看到外层是固体壳，中间有一段像烂泥一样的糊状区，内部还是液态——这就是糊状区的宏观表现。糊状区的宽度取决于合金成分（相图决定）：共晶成分几乎无糊状区（两相直接转变），宽成分区间的过共晶铝合金糊状区可达几十毫米。糊状区内液体流动（补缩流）不充分时就会在此形成缩孔（Shrinkage Porosity）——这是铸造凝固仿真最重要的预测目标之一。

通过 $f_l$ 对温度的导数，把潜热加到比热上：

在糊状区内 $\partial f_l/\partial T = 1/(T_l - T_s)$，所以：

这样现有的热传导求解器只需更换材料属性（将 $c_{eff}$ 替换 $c_p$），无需修改任何求解逻辑——这是商用CAE软件大多采用有效比热法的原因。以铝合金A356（$T_s=557°C, T_l=615°C, L_f=397\,\text{kJ/kg}$）为例，糊状区内等效比热 $c_{eff} = 0.96 + 397000/58 \approx 6840\,\text{J/(kg·K)}$，是显热比热的7倍多。

三种方法各有适用场景：

工程推荐：合金材料（有糊状区）→有效比热法（展宽糊状区3～5℃防振荡）；纯物质（尖锐熔点）→潜热源项法或焓法；CFD+凝固耦合→潜热源项法（自然集成对流项）。

Voller & Swaminathan（1991年）的经典方案，将潜热项移到方程右侧作为源项：

每个时间步的迭代流程：

1. 用当前温度 $T^k$ 更新液相分率 $f_l^k = f_l(T^k)$
2. 计算源项 $S^k = -\rho L_f (f_l^k - f_l^{k-1})/\Delta t$
3. 求解热传导方程：$\rho c_p (T^{k+1} - T^k)/\Delta t = \nabla\cdot(k\nabla T^{k+1}) + S^k$
4. 判断 $f_l^{k+1}(T^{k+1})$ 是否在 $[0,1]$ 范围内（防止过冲）
5. 重复步骤1～4直到收敛（$|f_l^{k+1}-f_l^k| < 10^{-4}$）

关键：步骤4的防过冲处理是稳定性的保障。如果 $f_l^{k+1} < 0$ 则固定为0并调整温度；如果 $f_l^{k+1} > 1$ 则固定为1。

几个实用稳定化技巧：

• 糊状区展宽（Smearing）： 把糊状区宽度适当展宽3～5℃，降低比热尖峰幅度。例如将A356从58℃展宽至65℃，$c_{eff}$ 峰值从6840降至6100 J/(kg·K)，振荡显著减小
• 液相分率亚松弛： $f_l^{new} = \omega f_l^{calc} + (1-\omega)f_l^{old}$（$\omega = 0.7 \sim 0.9$），防止液相分率在迭代间剧烈变化
• 自适应时间步： 检测当前 $\partial f_l/\partial t$，相变激烈时自动缩小 $\Delta t$（0.5～0.1倍），稳定后恢复
• 时间步Fourier约束： $\Delta t \le \Delta x^2 / (2\alpha)$（Fo ≤ 1），在糊状区内用 $\alpha = k/(\rho c_{eff})$，比显热区的 $\alpha$ 小很多
• 初始化策略： 从全液态（均匀高温）出发，用小时间步让凝固前沿稳定建立，再放大时间步

典型铸造凝固分析流程：

1. 材料数据收集： 固相线 $T_s$、液相线 $T_l$、熔解潜热 $L_f$、温度相关热导率 $k(T)$ 和比热 $c_p(T)$。从相图软件（Thermo-Calc、Pandat）提取或查文献手册。
2. 界面热阻（HTC）设置： 铸型-熔液界面气隙随时间增大，HTC从初始高值（金属型~2000 W/(m²K)，砂型~400 W/(m²K)）逐渐降低。这是影响凝固速度最大的参数。
3. 网格生成： 糊状区移动范围内加密（单元尺寸≤1～2 mm），铸型域也要包含在耦合网格内（导热耦合）。
4. 时间步设置： 初始 $\Delta t = 0.01 \sim 0.1\,\text{s}$（取Fo~0.1），启用自适应时间步，最大步长 $\le 1\,\text{s}$。
5. 结果提取： 凝固前沿速度 $v_s = -d(凝固分率)/dt$、凝固时间分布、温度梯度 $G$（影响晶粒细化）、局部凝固时间 $t_f$（预测缩孔）。
6. 验证： 与热电偶实测冷却曲线对比，确认凝固温度和凝固时间的误差在5～10%以内。

PCM蓄热分析（太阳能蓄热、建筑节能、电池热管理）与铸造凝固的主要区别：

• 循环充放热： PCM需要反复熔化-凝固循环，不像铸造是单次凝固。需要多周期仿真确认热响应的周期稳定性（通常3～5个周期后温度波形趋于稳定）
• 对流影响： 液态PCM中自然对流显著影响熔化速率（特别是大容量装置）。纯导热模型（焓法+Fourier方程）会低估熔化速度20～50%。液相区要用CFD（Navier-Stokes+Boussinesq近似）或至少加入等效对流导热率 $k_{eff} = Nu \cdot k_{liquid}$
• 体积变化： PCM熔化时体积膨胀（水凝固膨胀，有机PCM熔化膨胀），需要设计膨胀空间，否则容器会破裂
• 过冷问题： 部分PCM（如Na₂SO₄·10H₂O）凝固时容易过冷，实际凝固温度比熔点低5～20℃，需要在 $f_l(T)$ 中加入过冷度修正

主要工具对比：

糊状区里液体在固体树枝晶骨架间流动，类似多孔介质中的渗流，Carman-Kozeny方程描述这种流动阻力：

$A_{mush}$ 是糊状区常数（典型值 $10^4 \sim 10^7\,\text{kg/(m³s)}$），$\epsilon$ 是防除零小量（$\sim 10^{-3}$）。当 $f_l \to 0$（完全凝固），阻力项 $\to \infty$，速度被强制为零——等效于固体不流动。当 $f_l = 1$（完全液态），阻力为零，恢复正常流体。$A_{mush}$ 越大，糊状区流动越受限，适合细密树枝晶合金；纯金属可取较小值。

几个有趣的前沿方向：

• 相场法（Phase-Field）与焓法融合： 焓法处理宏观凝固（cm级），相场法处理枝晶生长（μm级）。将两者通过均一化连接，可以建立"枝晶尺度到铸件尺度"的多尺度相变模型，预测组织分布和力学性能。
• PINN辅助糊状区预测： 用物理信息神经网络（PINN）从多元合金相图数据学习非线性的 $f_l(T,\vec{c})$（$\vec{c}$ 是合金成分向量），比Scheil模型或Lever Rule更准确。
• 增材制造（AM）熔池分析： 激光粉末床熔融（LPBF）的快速凝固（$10^6$ K/s量级），远超传统铸造（$1\sim100$ K/s）。需要考虑非平衡凝固效应（成分分配系数变化、亚稳相形成），焓法结合快速凝固动力学是研究热点。
• 格子Boltzmann法（LBM）相变： LBM自然处理界面动力学，与焓法结合的LBM-Enthalpy方案在气泡/液滴相变（沸腾、冷凝）中精度优于传统VOF+焓法。
• 数字孪生中的实时相变预测： 将降阶焓法（ROM）嵌入生产线数字孪生，实时预测铸件内部凝固状态，用于自动化工艺优化。