建筑用PCM（相变蓄热材料）的热模拟
理论与物理
PCM是什么
老师，我听说建筑中会用到PCM，墙壁会融化吗？感觉有点可怕...
墙壁本身不会融化。PCM（相变储能材料）是一种在从固体变为液体时能吸收大量热量的材料。在建筑中，会将一种石蜡进行微胶囊化，制成直径5~50μm的颗粒，然后混入石膏板或砂浆中。
微胶囊？也就是说石蜡被包裹在壳里，只在内部融化或凝固，是吗？
没错。例如，巴斯夫的Micronal DS系列就是将熔点约23°C的石蜡包裹在三聚氰胺树脂壳中。白天，当室温超过23°C时，墙体内的PCM开始融化，将太阳热能作为潜热储存起来。夜晚室温下降时，它又会凝固并释放热量。通过这种"储热→放热"循环，可以将室温波动抑制2~4°C。
哦，墙壁变成了"热缓冲材料"啊。但是普通混凝土也能储热吧？PCM的优点是什么呢？
问得好。混凝土的储热是显热——即通过温度升高来储存。这样储热密度低。水的显热是每升高1°C储存4.18 kJ/kg，但石蜡C18的熔化潜热约为244 kJ/kg。以相同质量比较，PCM在仅1~2°C的温度变化范围内，就能储存相当于混凝土温度升高50°C所储存的热量。而且它是在室温"舒适区"附近工作，因此能直接减少空调负荷。
PCM在储放热过程中的总焓变，由固相的显热+潜热+液相的显热表示。
其中 $m$ 是质量，$c_s, c_l$ 是固相/液相比热，$T_m$ 是熔点，$L_f$ 是熔化潜热，$T_1, T_2$ 是初始/最终温度。建筑用石蜡系PCM的代表性物性值如下所示。
| 物性 | 值 | 单位 |
|---|---|---|
| 熔点 $T_m$ | 21~26 | °C |
| 熔化潜热 $L_f$ | 170~250 | kJ/kg |
| 密度 $\rho$（固体） | 850~950 | kg/m³ |
| 热导率 $k$ | 0.15~0.25 | W/(m·K) |
| 比热 $c_p$（固/液） | 1.8~2.4 | kJ/(kg·K) |
Stefan问题与移动边界
PCM在融化或凝固时，固相和液相的边界会移动对吧？这个怎么公式化呢？
那就是Stefan问题。由奥地利物理学家Jozef Stefan于1891年公式化，是伴随相变的移动边界问题。在固相区域和液相区域分别求解热传导方程，并满足界面 $x = s(t)$ 处的能量守恒条件（Stefan条件）。一维情况下如下所示。
固相区域 ($0 < x < s(t)$):
液相区域 ($s(t) < x < L$):
Stefan条件（界面 $x = s(t)$ 处的能量平衡）:
界面位置 $s(t)$ 是未知数，而且会随时间移动呢。这个能解析求解吗？
在单侧加热半无限固体等特殊条件下，存在 $s(t) = 2\lambda\sqrt{\alpha t}$ 这样的Neumann解析解。但对于建筑围护结构这样的有限厚度、两侧温度变化、多层材料的问题，解析解是不可能的。必须使用数值解法，这时焓法就登场了。
焓法的控制方程
焓法是求解Stefan问题的便捷方法吗？思路是怎样的呢？
Stefan问题的棘手之处在于"界面位置的追踪"。焓法是避开这个问题的天才想法，它以体积焓 $H$ 代替温度 $T$ 作为因变量。由于相变前后 $H(T)$ 函数会因潜热而产生跳跃，因此可以用单一方程统一处理固相、液相和混合相（糊状区）。无需显式地追踪界面。
焓法的控制方程可以写成如下形式。
其中体积焓 $H(T)$ 定义为温度的函数。
对于建筑用PCM，相变并非发生在单一温度，而是在 $T_s$ 到 $T_l$ 的温度范围（通常为2~4°C宽）内发生。这个范围就是糊状区，是固相和液相共存的区域。
原来如此！用温度写方程在界面处会不连续而难以处理，但用焓写就能平滑处理，是这个意思吧。
没错。当糊状区宽度 $\Delta T = T_l - T_s$ 趋近于零时，就归结为理想的"尖锐界面"Stefan问题。实际的建筑用PCM具有 $\Delta T = 2\sim4$°C的有限宽度，因此焓法反而是更自然的公式化方法。
PCM围护结构的传热机制
实际的建筑围护结构，不仅有PCM层，还有隔热材料、混凝土等各种层，对吧？整体上如何建模呢？
用一个例子来说明典型的PCM围护结构构成吧。从外侧开始，通常是"外饰材料（10mm）→ 隔热材料（50mm）→ PCM石膏板（12.5mm）→ 空气层→ 内饰面"这样的多层结构。需要联立求解各层的一维非稳态热传导。
对于多层围护结构的每一层 $i$:
层间连接条件（温度连续性和热流密度守恒）:
外表面/内表面的边界条件包括对流（$h_\mathrm{ext}$, $h_\mathrm{int}$）、太阳辐射吸收（$\alpha_s I_\mathrm{sol}$）和长波辐射（$\varepsilon \sigma (T^4 - T_\mathrm{sky}^4)$）。
也就是说，只有PCM层的 $c_p(T)$ 是温度依赖的，变成非线性了。其他层保持通常的热传导方程就行，对吧。
是的。所以计算上的要点就集中在"如何细致地离散化PCM层"。隔热材料和混凝土即使划分得粗一些也足够，但PCM层在相变温度范围 $\Delta T$ 内焓值变化剧烈，因此需要精细划分。
Stefan问题的历史
Stefan问题的名称来源于1891年提出该问题的奥地利物理学家Jozef Stefan（1835-1893），他也是玻尔兹曼的老师。他当时研究的问题是"北极海的冰能生长到多厚"。在给定冰表面温度和海水温度的情况下，冰的厚度 $s(t)$ 与 $\sqrt{t}$ 成比例增长的Neumann解，至今仍被用作PCM模拟的验证基准问题。建筑领域的应用真正开始于1980年代，起因是美国能源部（DOE）的Doris Boehm等人的研究。
各项的物理意义
- 焓积累项 $\partial H / \partial t$：单位体积的焓变化率。包含显热（温度变化）和潜热（相变）。在PCM熔点附近，$H$ 会随微小的温度变化而急剧增加——这就是"吸热"效果。用日常例子来说，冰水能在0°C保持很长时间，就是因为冰在持续吸收熔化潜热。
- 热传导项 $\nabla \cdot (k \nabla T)$：基于傅里叶定律的热扩散。PCM的热导率很低，为0.15~0.25 W/(m·K)，大约是混凝土（1.4 W/(m·K)）的1/6。因此，PCM层具有"擅长储热，但传热慢"的特性。在围护结构设计中，重要的是在发挥PCM储热能力的同时，确保向内侧的放热路径，取得平衡。
- 热源项 $Q$：内部发热（在PCM围护结构中通常 $Q=0$）。不过，对于内置电加热器的主动型PCM面板，有时需要考虑焦耳发热 $Q = I^2R/V$。
假设条件与适用范围
- 一维传热假设：当围护结构面积相对于厚度足够大时，可以忽略面内方向的热传导，作为仅厚度方向的一维问题处理。窗框周围和角落部分会形成热桥，因此需要2D/3D分析。
- 自然对流忽略：在微胶囊型PCM中，液相被封闭在微小胶囊内，因此不会产生宏观的自然对流。但是，对于块状PCM（如储热罐等），液相的对流占主导地位，需要与Navier-Stokes方程耦合求解。
- 体积变化忽略：石蜡在相变时伴随约10%的体积变化，但通常被微胶囊的弹性变形吸收，因此结构影响通常可以忽略。
- 滞后忽略：实际的PCM在熔化和凝固时表现出不同的温度-焓曲线（滞后现象），但在大多数BES工具中被忽略。精确评估需要滞后模型。
量纲分析与无量纲数
| 无量纲数 | 定义 | 物理意义 | PCM围护结构中的典型值 |
|---|---|---|---|
| Stefan数 $\mathrm{Ste}$ | $c_p \Delta T / L_f$ | 显热与潜热之比 | 0.02~0.1（潜热主导） |
| Biot数 $\mathrm{Bi}$ | $h L / k$ | 表面对流与内部传导之比 | 0.5~5 |
| Fourier数 $\mathrm{Fo}$ | $\alpha t / L^2$ | 无量纲时间（扩散进行度） | 年周期下 $\gg 1$ |
数值解法与实现
焓法的离散化
要实际用计算机求解焓法，如何进行离散化呢？
在建筑的BES（建筑能耗模拟）工具中，主流做法是沿厚度方向用有限差分法（FDM）对围护结构进行离散化。以EnergyPlus的CondFD算法为例，将PCM层分割为 $N$ 个节点，在每个节点 $j$ 处求解以下隐式差分方程。
$H$ 和 $T$ 都是未知数，所以需要用 $H(T)$ 的关系式联系起来，对吧？但这是非线性的吧？
很敏锐。每个时间步都需要通过反函数 $T = H^{-1}(H)$ 从 $H^{n+1}$ 求出 $T^{n+1}$。具体来说，首先以上一个时间步的温度作为初始估计值，在 $H(T)$ 曲线上进行更新温度的反迭代计算。在糊状区内，可以用 $T = T_s + (H - H_s)\Delta T / (\rho L_f)$ 显式求解，因此迭代通常1~3次就能收敛。
与表观比热法的比较
我也听说过"表观比热法"，它和焓法有什么区别呢？
表观比热法（Apparent Heat Capacity法）是将潜热表现为温度依赖的等效比热 $c_\mathrm{eff}(T)$ 的方法。
使用表观比热法，就可以作为通常的热传导方程 $\rho c_\mathrm{eff}(T) \partial T / \partial t = \nabla \cdot (k \nabla T)$ 来求解，因此易于在现有求解器中实现。但是有一个很大的陷阱——如果时间步长太粗，可能会"跳过"相变温度范围，导致潜热完全从计算中漏掉。
诶，那太可怕了！也就是说，如果一步就从 $T_s$ 跳到 $T_l$，潜热部分的能量就消失了，是吗？
没错。这就是表观比热法最大的弱点。另一方面，焓法直接追踪 $H$，因此即使时间步长稍粗，也能保证潜热守恒。实际工作中按以下方式区分使用。
| 方法 | 优点 | 缺点 | 推荐场景 |
|---|---|---|---|
| 焓法 | 能量守恒严格，鲁棒性强 | 需要 $H \to T$ 的反变换 | BES工具（如EnergyPlus等） |
| 表观比热法 | 易于在现有求解器中实现 | 依赖时间步长，有漏掉潜热的风险 | CFD/FEM（如COMSOL、Fluent等） |
时间步长与网格划分
对于PCM围护结构的分析，合适的时间步长和网格细密程度大概是多少呢？
这是相当关键的一点。PCM的热扩散率非常小，$\alpha = k/(\rho c_p) \approx 8 \times 10^{-8}$ m²/s。对于12.5mm厚的PCM石膏板，特征时间为 $\tau = L^2/\alpha \approx 2400$ 秒（40分钟）。以此为标准：
- 时间步长: 1~3分钟（年模拟中3分钟较实用，精密评估则用1分钟）
- 空间划分: PCM层至少划分4份（$\Delta x \le 3$ mm），推荐6~8份
- Fourier数约束: $\mathrm{Fo} = \alpha \Delta t / \Delta x^2 \le 0.5$（显式法的情况。隐式法则无约束）
年模拟的话，365天 $\times$ 1440分钟 = 525,600步... 即使3分钟一步也有175,200步呢。计算量怎么样？
因为是1D围护结构的有限差分，所以单步计算成本很小。在EnergyPlus中，对于包含PCM围护结构的单区域模型年模拟，也只需要几十秒到几分钟。即使是多区域整体建筑，也只需10分钟左右。如果用COMSOL进行2D/3D详细分析，计算量会大得多，但除非是想观察围护结构内部的自然对流或界面附近的细节，否则1D-BES就足够了。
非线性收敛的处理
PCM层的非线性不会引起收敛问题吗？
问得好。PCM的 $H(T)$ 曲线在糊状区内斜率很大，这会导致热容 $C = dH/dT$ 急剧变化。如果使用牛顿-拉夫逊法，雅可比矩阵的条件数会变差，可能引起振荡或发散。因此，BES工具通常采用以下策略：
なった
詳しく
報告