地暖系统模拟 — 埋设管道热传导分析与地板表面温度分布预测
理论与物理
概要 — 地板采暖CAE的整体概览
老师，地板采暖的模拟是分析什么的呢？不就是在地板下面埋管子吗？
简单来说，就是用FEM求解埋设管道的热传导，来预测地板表面的温度分布。虽说“只是埋管子”，但其实管道布置方式会极大地影响舒适度。
诶，是这样吗？具体会怎么变化呢？
例如，管道间距200mm和300mm，地板表面的温度不均匀性（最高温度与最低温度之差 $\Delta T$）会有很大不同。200mm间距时 $\Delta T \approx 1\text{--}2°\text{C}$ 就够了，但300mm时 $\Delta T \approx 3\text{--}5°\text{C}$。脚底会感觉到“这里冷，那里热”的程度。
原来如此…那日本住宅的目标温度大概是多少呢？
日本的建筑标准推荐地板表面温度在29°C以下。WHO的指南也建议居室地板表面温度在19〜29°C为宜。温度过高有低温烫伤的风险，过低则失去了采暖的意义。为了瞄准这个“恰到好处的区间”，CAE被用于管道间距、埋设深度、隔热材料厚度的组合优化。
三种热传递机制
地板采暖，地板会变暖我明白，但整个房间是怎么变暖的呢？不像空调那样有风吹出来吧？
问得好。地板采暖的热传递由三种机制构成：
- 传导：热水管道 → 砂浆 → 地板材料，热量在固体内部传递
- 对流：从温暖的地板表面到室内空气的自然对流。暖空气上升
- 辐射：地板表面以红外线形式向墙壁、天花板、家具辐射热量。实际上总散热量的50〜70%是这种辐射
辐射占比那么大吗！这和篝火前感觉温暖是同一个原理吗？
正是如此。空调通过对流加热空气，所以容易导致天花板附近很热而脚下寒冷——即所谓的“头热足寒”。地板采暖则相反，实现“头寒足热”。Fanger的舒适性研究（1972年）也表明，足部24°C、头部19°C对人体来说最为舒适。因此，在CAE中辐射分量的建模极为重要。
控制方程
基本是三维非稳态热传导方程。求解地板板（混凝土+砂浆+饰面材料）内部的温度分布：
是普通的热传导方程呢。但是管道中流动的热水怎么处理呢？
管道内的热水通常用一维的对流方程来处理。如果将所有管道的弯曲和分支都用3D CFD求解，计算成本会非常庞大，因此在实际应用中，“1D管道+3D固体传导”的耦合模型是主流：
这里 $T_w$ 是水温，$v$ 是流速，$s$ 是管道长度方向坐标，$A_p$ 是管道截面积，$d_i$ 是管道内径，$q_p$ 是管道壁面单位长度的散热量。这样就能追踪热水在管道中流动时逐渐冷却的过程。
啊，所以管道入口附近和出口附近的地板温度会不一样对吧！
没错。在热水式地板采暖中，入口附近的热水温度较高（通常40〜60°C），出口附近会下降5〜10°C。这成为“管道布置模式”设计中的重要点。如果采用螺旋型（涡旋状），入口管道和出口管道会交替排列，从而使温度不均匀性得到平均化。
管道布置的热阻模型
从热水到地板表面，热量会受到多大的阻力呢？
从热水到室内空气的总热阻可以整理为“串联电路”。和电路中的欧姆定律是同样的思路：
| 热阻 | 物理意义 | 典型值 [m²·K/W] |
|---|---|---|
| $R_{conv,w}$ | 热水→管道内壁的对流 | 0.001〜0.003 |
| $R_{pipe}$ | PEX管壁的传导 | 0.003〜0.01 |
| $R_{mortar}$ | 砂浆埋设层的传导 | 0.02〜0.06 |
| $R_{finish}$ | 饰面材料（地板等）的传导 | 0.04〜0.15 |
| $R_{conv,r}$ | 地板表面→室内空气的对流 | 0.10〜0.15 |
| $R_{rad}$ | 地板表面→墙壁·天花板的辐射 | 0.05〜0.10 |
饰面材料的阻力意外地大呢。地板材料的厚度和材质不同，差别会很大吧…
很敏锐。实际上，地板材料的选择是地板采暖设计中最容易被忽视的点。实木（$k \approx 0.12$ W/(m·K)）和瓷砖（$k \approx 1.0$ W/(m·K)）的热导率相差约8倍。瓷砖更容易将管道的热量传递到地板表面，因此即使热水温度相同，地板表面温度也会更高。欧洲浴室多采用瓷砖+地板采暖的组合就是这个原因。
这里 $h_{conv} \approx 5\text{--}7$ W/(m²·K)（向上表面的自然对流），$h_{rad} \approx 5\text{--}6$ W/(m²·K)（室温20°C附近的线性化辐射系数）。合计 $h_{total} \approx 10\text{--}13$ W/(m²·K) 是实际应用中的标准值。
舒适性指标与CAE的结合点
“舒适”我觉得是主观的东西，用模拟怎么评价呢？
ISO 7730定义的PMV（预测平均投票）指标是标准。它从温度、湿度、气流速度、辐射温度、衣着量、代谢量这6个要素出发，将舒适性量化为“热（+3）〜冷（-3）”的尺度。将地板采暖的分析结果（地板表面温度分布、室内气流场）输入到PMV计算中，就可以对房间哪里舒适、哪里不舒适进行空间映射。
CAE的温度计算结果可以直接用于人的舒适性评价呢。好厉害。
特别是地板采暖，ASHRAE 55对地板表面温度的不均匀性（$\Delta T_{floor}$）有限制。坐着的人脚底如果有3°C以上的温差就会感到不适。通过CAE进行这项检查，可以在施工前判断管道布置模式的好坏。
古罗马的地坑式供暖
地板采暖的历史可以追溯到公元前80年左右的古罗马。被称为“地坑式供暖（hypocaust）”的系统，是在地板下设置空洞，让炉灶的热气循环的构造。广泛应用于浴场（thermae），通过改变地板下支柱（pilae）的高度来控制温度分布。韩国的“温突”也是同样的原理，拥有超过1000年的历史。现代CAE工程师解决的本质问题——“如何才能均匀地温暖地板表面”——早在2000多年前的人们就已经凭经验解决了。
各项的物理意义
- 蓄热项 $\rho c_p \partial T/\partial t$：表示混凝土板的蓄热性。混凝土比热约880 J/(kg·K)，密度约2,300 kg/m³，数值较大，因此加热需要时间，但一旦变热就能长时间保温。地板采暖“升温时间长（30分钟〜2小时）”问题的根源就在这里。
- 热传导项 $\nabla \cdot (k \nabla T)$：地板材料内部的热扩散。砂浆（$k \approx 1.5$ W/(m·K)）与混凝土（$k \approx 1.6$ W/(m·K)）数值相近，但如果夹有木质地板（$k \approx 0.12$ W/(m·K)），则会产生较大的温度降。隔热材料（XPS: $k \approx 0.035$ W/(m·K)）则作为有意不让热量通过的层发挥作用。
- 热源项 $Q$：对于电热式地板采暖，直接建模为加热线的焦耳发热 $Q = I^2 R / V_{heater}$（W/m³）。对于热水式，则作为热水管道管壁的边界条件间接表达。
量纲分析与单位制
| 变量 | SI单位 | 地板采暖中的典型值 |
|---|---|---|
| 热水温度 $T_w$ | °C | 35〜60°C（低温水式：35〜45°C） |
| 地板表面温度 $T_s$ | °C | 24〜29°C（推荐上限29°C） |
| 热导率 $k$（砂浆） | W/(m·K) | 1.3〜1.6 |
| 热导率 $k$（地板） | W/(m·K) | 0.10〜0.17 |
| 管道间距 | mm | 150〜300（标准200mm） |
| 管道埋设深度 | mm | 30〜60 |
| 散热量 $q$ | W/m² | 50〜100（标准住宅采暖负荷） |
数值解法与实现
基于FEM的离散化
刚才的热传导方程，实际上是怎么用计算机求解的呢？
用有限元法（FEM）对空间进行离散化。将地板板分割成许多小单元，近似求解每个单元内的温度。转换成弱形式（变分形式）并用Galerkin法进行公式化，最终得到以下联立方程：
这里 $[C]$ 是热容矩阵，$[K]$ 是热传导矩阵，$\{F\}$ 是外部热输入向量。和结构分析的 $[K]\{u\} = \{F\}$ 形式很像吧？热分析中温度是未知数，只是用热传导矩阵代替了刚度矩阵。
结构分析的知识可以直接用上呢。用什么单元呢？
地板采暖分析中，8节点六面体单元（线性）是标配。地板板是平板状，与六面体单元很匹配。只在管道周边细化，远离地板表面的部分则粗化——这就是所谓的不均匀网格的基本做法。
| 单元类型 | 节点数 | 地板采暖中的用途 | 注意事项 |
|---|---|---|---|
| 8节点六面体（线性） | 8 | 地板板主体 | 最通用。最适合层状结构的扫掠网格 |
| 20节点六面体（二次） | 20 | 管道周边的高精度区域 | 在温度梯度陡峭处有效 |
| 4/10节点四面体 | 4/10 | 复杂形状的管道分支部分 | 用于无法使用六面体部分的补充 |
| 1D管道单元 | 2 | 热水配管 | COMSOL: Pipe Flow, Fluent: 1D pipe network |
管道建模方法
把管道按3D圆形切出来划分网格，看起来超级麻烦啊…
实际应用中有三种方法：
- 全3D解析：忠实网格化管道的圆形截面。最精确但计算成本巨大。用于研究
- 等效热源法：从3D模型中省略管道，在管道位置的节点上施加等效的热输入。计算快但温度梯度精度下降
- 1D-3D耦合法：将1D管道单元嵌入到3D固体单元中。COMSOL的“Pipe Flow”模块和Ansys Fluent的“embedded pipe”就是这种。实际应用的主流
1D-3D耦合是最平衡的呢。管道的拐角之类的也能追踪吗？
当然可以。1D管道单元沿着管道路径配置节点，因此无论是蛇行型还是螺旋型都能自由追踪。想象一下，每个节点处管道外壁与周围的3D固体单元在热学上结合。管道内的水温在哪里下降了多少也能实时计算出来。
CFD+传导耦合分析
刚才提到“室内空气的自然对流”对吧。那个也要模拟吗？
需要高精度时会做。用CFD（计算流体力学）求解室内空气的流动和温度场，在地板表面的边界上与热传导分析进行耦合。这就是所谓的共轭传热。
什么时候需要做到CFD呢？总是需要吗？
实际上，在大多数情况下不需要CFD。只需将地板表面对流换热系数 $h_{conv}$ 设为常数（5〜7 W/(m²·K)）就能获得足够的精度。需要CFD耦合的情况是：
- 天花板非常高的空间（挑高空间、体育馆）且气流模式复杂时
- 需要评估窗户冷风渗透与地板采暖的相互作用时
- 需要考虑家具布置对气流的阻碍时
- 需要高精度空间映射PMV分布时
原来如此，先只用传导求解，必要时再追加CFD，是分阶段的方法呢。
时间步长与非稳态分析
听说地板采暖“变暖需要时间”，那个升温时间也能模拟吗？
当然可以。通过非稳态（瞬态）分析追踪时间演变。时间方向的离散化通常使用隐式解法（后向欧拉法或Crank-Nicolson法）。隐式解法即使取较大的时间步长也能稳定，因此数小时的模拟只需几十个时间步就能完成。
时间步长的基准，是根据地板板热扩散的特征时间 $\tau$ 来确定的：
なった
詳しく
報告