蓄热材料（PCM）的热模拟
理论与物理
什么是PCM
老师，PCM就是蓄热材料对吧？在仿真中具体要解什么呢？
问得好。PCM（相变材料）是一种在从固体变为液体时——也就是熔化时——能吸收大量热量的材料。利用这种“潜热”来控制温度，正是蓄热的本质。
用潜热控制温度…？能再具体点吗！
例如，冰在0°C熔化时，温度不会上升，却能吸收334 kJ/kg的热量，对吧？PCM就是利用这个原理。石蜡系材料在熔点28°C附近潜热约200 kJ/kg，水合盐在熔点32°C附近潜热约260 kJ/kg。为了将电池包温度维持在40°C以下，近年来在电芯之间填充PCM的设计急剧增加。
原来如此，熔化期间温度不上升，所以就像一道“温度墙”对吧！那仿真要解的是…？
预测固体与液体的边界——熔化前沿如何移动、何时何处会熔化、完全熔化后温度如何上升。这就是被称为“斯蒂芬问题”的经典移动边界问题。在CAE中，主流方法是使用焓法，无需显式追踪边界即可求解。
斯蒂芬问题与移动边界
斯蒂芬问题这名字我听说过。就是熔化位置不断移动的问题对吧？
没错。我们来考虑一维的斯蒂芬问题。假设有一块单面被高温加热的PCM平板。设固液界面的位置为 $s(t)$，则在界面处满足“吸收潜热的同时界面前进”这一条件：
这里 $L$ 是单位质量潜热 [J/kg]，$k_s, k_l$ 分别是固相和液相的热导率，$T_s, T_l$ 是固相和液相的温度场。界面两侧热通量的差值，成为熔化的能量来源。
追踪界面的话，网格也得跟着动，看起来好麻烦…
说得对。所以在实际工作中，绝大多数情况下都使用回避界面追踪的焓法。其思路是“无需知道界面在哪里，只要能量收支平衡，温度场就能正确得出”。
焓法的控制方程
请告诉我焓法的公式！
焓法中，求解以单位体积焓 $H(T)$ 为变量的能量守恒方程：
这里焓 $H(T)$ 作为温度的函数定义如下：
这里 $f_l(T)$ 是液相率（liquid fraction），是一个从 $0$（完全固体）变化到 $1$（完全液体）的函数。$\epsilon$ 是糊状区的半宽，对于纯物质 $\epsilon \to 0$，对于合金或石蜡混合物则为几°C左右。
想象一下 $H(T)$ 的图，是在熔点处突然跳跃的感觉吗？
答对了！对于纯物质，在熔点处会垂直跳跃。对于像石蜡或水合盐这样有熔化温度范围的材料，则会呈现倾斜的阶梯状。这个跳跃的高度就对应着潜热 $L$。在仿真中，这种急剧变化是非线性的根源，因此需要注意时间步长和网格。
各项的物理意义
- 蓄热项 $\rho \partial H / \partial t$：单位体积的焓变化率。不仅包含显热，也包含潜热。熔化过程中即使温度几乎不变，$H$ 也会发生很大变化。【具体示例】石蜡（RT28HC）在28°C熔化时，$H$ 增加200 kJ/kg，但温度几乎不变。
- 扩散项 $\nabla \cdot (k \nabla T)$：基于傅里叶定律的热传导。PCM的热导率通常较低（石蜡系 $k \approx 0.2$ W/(m·K)）。这是PCM最大的弱点，随着熔化进行，液相层会起到隔热层的作用，导致熔化速度急剧下降。
- 热源项 $Q$：电池电芯发热或加热器输入等。对于电池包，$Q = I^2 R_{internal} / V_{cell}$ 用于评估电芯的内阻发热。
量纲分析与单位制
| 变量 | SI单位 | 注意点 |
|---|---|---|
| 焓 $H$ | J/kg | 从参考温度开始的累积值。根据求解器不同，定义可能有所不同 |
| 潜热 $L$ | J/kg | 石蜡: 150〜250，水合盐: 200〜300，金属: 200〜400 |
| 液相率 $f_l$ | 无量纲 [0, 1] | 0=完全固体，1=完全液体。在糊状区内连续变化 |
| 热导率 $k$ | W/(m·K) | 石蜡固相: 0.24，液相: 0.15。水合盐: 0.5〜1.1 |
| 斯蒂芬数 $\mathrm{Ste}$ | 无量纲 | $c_p \Delta T / L$。Ste < 1 为潜热主导，Ste > 1 为显热主导 |
等效比热法（Apparent Heat Capacity）
除了焓法还有其他方法吗？
另一个著名的方法是等效比热法（Apparent Heat Capacity Method）。它将潜热建模为熔化温度范围内的表观比热 $c_{p,eff}$：
在实现上，用高斯分布或帽函数近似狄拉克δ函数，表现为宽度为 $2\epsilon$ 的峰值：
焓法和等效比热法，哪个更好呢？
从实用性来说，焓法具有压倒性优势。等效比热法在时间步长较粗时，可能会跳过熔化温度区间而漏掉潜热——即发生所谓的“潜热跳过”。Ansys Fluent和COMSOL的相变模块在内部都使用焓法，等效比热法是在复用FEM标准热分析求解器时的简易方法。
| 特性 | 焓法 | 等效比热法 |
|---|---|---|
| 能量守恒性 | ◎ 严格守恒 | △ 依赖时间步长 |
| 潜热跳过风险 | 无 | 有（$\Delta t$ 大时发生） |
| 实现容易度 | ○ 需要专用模块 | ◎ 可用于标准热分析 |
| 糊状区的表现 | 可自然表现 | $\epsilon$ 的选择很关键 |
| 主要采用的求解器 | Fluent, STAR-CCM+, COMSOL | Abaqus, Nastran（通过UDL） |
斯蒂芬数 Ste 的物理意义
斯蒂芬数是什么？怎么用呢？
斯蒂芬数是表示显热与潜热之比的无量纲数。它是量化PCM仿真“难度”的指标：
这里 $\Delta T$ 是过热度（例如壁面温度与熔点的差值 $T_{wall} - T_m$），$L$ 是潜热。
Ste 大了会有什么变化呢？
粗略来说：
- $\mathrm{Ste} \ll 1$（潜热主导）：熔化缓慢进行。电池包的PCM冷却属于此区域。若 $\Delta T = 12$°C，$c_p = 2000$ J/(kg·K)，$L = 200{,}000$ J/kg，则 $\mathrm{Ste} = 0.12$。时间步长可以取得相对较大。
- $\mathrm{Ste} \gg 1$（显热主导）：几乎接近普通的热传导问题。铸造凝固（$\Delta T$ 数百°C）属于此类。时间步长若不细化，容易产生数值振荡。
- $\mathrm{Ste} \sim 0.1$：可以使用纽曼近似解的区域。熔化前沿位置可估计为 $s(t) \approx 2\lambda\sqrt{\alpha t}$，便于验证分析结果。
PCM材料的种类与物性值
PCM材料除了石蜡还有其他种类吗？
大致可分为三类。连同仿真所需的物性值一起整理如下：
| PCM类别 | 代表材料 | 熔点 [°C] | 潜热 [kJ/kg] | $k$ [W/(m·K)] | 用途 |
|---|---|---|---|---|---|
| 有机系（石蜡） | RT28HC, n-十八烷 | 24〜32 | 180〜250 | 0.2 | 电池冷却、建筑 |
| 无机系（水合盐） | CaCl₂·6H₂O, Na₂SO₄·10H₂O | 29〜32 | 170〜265 | 0.5〜1.1 | 地暖、蓄冷 |
| 共晶系 | 有机+有机、有机+无机混合 | 设计自由 | 120〜200 | 0.2〜0.8 | 特定温度区域设计 |
| 金属系 | Ga, Bi-Sn合金 | 30〜140 | 60〜80 | 20〜40 | 高热通量电子设备 |
石蜡的热导率0.2，这也太低了吧…
是的，这正是PCM仿真最大的实际课题。热导率低会导致熔化速度极慢。所以在实际产品中，会加入铝翅片或分散碳纳米管来提高有效热导率。在仿真中，是将其作为“PCM+翅片”的复合体处理，还是使用均质化的有效物性值，这里正是设计者展现本领的地方。
NASA的PCM温度控制
NASA从1960年代的阿波罗计划就开始将PCM用于航天器的温度控制。月球表面有向阳+127°C / 背阴-173°C 的严酷温度环境，但搭载了石蜡系PCM的月球车电子设备壳体，成功将内部温度维持在-20〜+50°C的范围内。2020年代的火星探测器Perseverance也搭载了改进型PCM模块。其仿真必须包含真空中的辐射散热耦合分析，这正是CAE大显身手之处。
数值解法与实现
空间离散化与时间积分
如何在计算机上求解焓法呢？
如果是FVM（有限体积法），可以将每个单元的焓作为守恒量直接处理，所以兼容性好。Ansys Fluent和STAR-CCM+采用这种方式。FEM（有限元法）则需要先转换为弱形式再进行离散化。
展示基于FVM的离散化方程。对于单元 $i$ 的能量守恒方程进行隐式离散化得到：
这里 $V_i$ 是单元体积，$A_f$ 是面面积，$\delta_f$ 是单元中心间距离，下标 $nb$ 表示相邻单元。由于 $H_i^{n+1}$ 和 $T_i^{n+1}$ 之间存在 $H(T)$ 关系，因此需要用牛顿法或基于源项的线性化方法迭代求解。
时间积分用隐式格式和显式格式哪个好呢？
对于PCM，隐式格式（后向欧拉）是铁律。原因有二：
- 石蜡热导率低，要满足傅里叶数 $\mathrm{Fo} = \alpha \Delta t / \Delta x^2$ 的稳定性条件 $\mathrm{Fo} \leq 0.5$，需要非常小的时间步长
- 焓的急剧变化（潜热跳跃）会产生强非线性，而隐式格式可以在每个时间步内迭代收敛
糊状区的处理
经常听到“糊状区”，那是什么？
既不是固体也不是液体的中间状态——固液共存的区域。对于石蜡混合物，熔点范围在28〜32°C之间，所以整个这个温度带都是糊状区。
Ansys Fluent的 Solidification/Melting 模型中，糊状区的液相率 $f_l$ 按如下线性插值：
此外，为了抑制糊状区内的流动，在运动方程中添加Darcy型源项（焓-孔隙度法）：
这里 $C_{mush}$ 是糊状区常数（典型值 $10^5$ 〜 $10^8$），$\varepsilon_{mush}$ 是防止除零的微小常数（通常 $10^{-3}$）。当 $f_l \to 0$（固体）时 $S \to \infty$（完全阻止流动），$f_l \to 1$（液体）时 $S \to 0$（自由流动）。
$C_{mush}$ 的值怎么确定呢？ $10^5$ 和 $10^8$ 差太多了吧…
实际上，$C_{mush}$ 的选择是PCM仿真的“手艺活”部分。纯石蜡多用 $10^5$ 〜 $10^6$，并根据实验的熔化前沿位置进行调整。$C_{mush}$ 太大，糊状区会变得太硬，熔化变慢；太小，则会在固相区域产生非物理的流动。论文中一定会进行“$C_{mush}$ 灵敏度分析”，建议尝试三个以上水平的值。
液相自然对流建模
熔化后的液体会流动吗？只用热传导不行吗？
流动得非常厉害。这正是PCM仿真的难点所在。熔化的石蜡因温差产生密度差，从而引发浮力驱动的自然对流。当瑞利数超过 $10^5$ 时，对流效应将占主导地位，熔化速度可能比纯传导预测的快2〜5倍。
使用Boussinesq近似…
なった
詳しく
報告