MCM520 ← 资料站首页 2018A 高温作业专用服装（一）：四层非稳态热传导建模与参数辨识 打开交互阅读器 →

2018A 高温作业专用服装（一）：四层非稳态热传导建模与参数辨识

1 问题重述

2018 年国赛 A 题要求为高温作业专用服装建立传热模型，并优化各层厚度，使体表（皮肤）温度在作业过程中不超过安全阈值。本题给出四层织物（外→内依次为阻燃外壳、隔热层、舒适层、空气隙），其中第Ⅱ层（隔热层）与第Ⅳ层（空气隙）为设计变量，而第Ⅰ、Ⅲ层固定。题目要求：①在高温外环境（如 65°C、75°C、80°C）下，建立多层一维非稳态热传导模型，预测皮肤温度随时间的变化；②以"皮肤温度 ≤ 47°C（安全）/ ≤ 44°C（灼痛预警）"为约束，求解满足约束的最小总厚度；③考虑材料与环境的不确定性，给出稳健设计方案。

本篇聚焦第一与第三环节：先把物理模型建对、验证对，再用一组 75°C 假人实验数据把未知参数（皮肤换热系数 h_skin 与隔热层导热系数 k2）辨识出来。只有正演模型可信，后续的优化与稳健设计才站得住。许多参赛队伍在这一题上失分，正是因为一上来就调用现成偏微分方程求解器，却没验证求解器在多层突变界面、瞬态短时工况下的正确性——结果优化建立在错误的温度场上。

2 基本假设与符号

为把问题简化为可计算的一维模型，作如下假设：

  • (H1) 热量沿垂直于服装表面的方向一维传导，忽略织物平面内扩散与四肢弯曲带来的几何复杂性；
  • (H2) 每层材料均匀、各向同性，物性参数（导热系数 k、密度 ρ、比热 c）取定值，不随温度变化；
  • (H3) 外表面与高温环境直接接触，温度被环境强制拉到 T_env（Dirichlet 边界）；
  • (H4) 皮肤面与体核通过血流与组织对流换热，用等效换热系数 h_skin 描述为 Robin 边界；
  • (H5) 初始时刻整层温度等于体核温度 T0 = 37°C，即穿着前服装已与人体热平衡。

记第 j 层厚度为 LjL_j、导热系数为 kjk_j、体积热容 ρjcj\rho_j c_j。四层物性（基线设计）为：

层 名称 厚度(mm) k (W/(m·K)) ρc (J/(m³·K))
Ⅰ 阻燃外壳(固定) 0.6 0.082 300×1377
Ⅱ 隔热层(设计) L2 0.370 862×2100
Ⅲ 舒适层(固定) 3.6 0.045 74.2×1726
Ⅳ 空气隙(设计) L4 0.028 1.18×1005

皮肤等效换热系数取 hskin=14 W/(m2⋅K)h_{skin}=14\ \mathrm{W/(m^2\cdot K)}，体核温度 Tcore=37∘CT_{core}=37^\circ\mathrm{C}，网格节点数 N=60N=60，时间步长 Δt=30 s\Delta t=30\ \mathrm{s}。需要强调，h_skin 并非单纯的导热系数，它把皮肤内部的血流对流散热一并折算进来，因此是个体间差异最大的参数之一——这也正是第三篇蒙特卡洛中它成为首要不确定源的原因。

3 传热控制方程与数值方法

3.1 控制方程

沿厚度方向 x∈[0,L]x\in[0,L]，温度场 T(x,t)T(x,t) 满足一维非稳态热传导方程

ρc ∂T∂t=∂∂x ⁣(k ∂T∂x)\rho c\,\frac{\partial T}{\partial t}=\frac{\partial}{\partial x}\!\left(k\,\frac{\partial T}{\partial x}\right)

边界条件为：x=0x=0（外表面）给定 T(0,t)=TenvT(0,t)=T_{env}；x=Lx=L（皮肤面）采用 Robin 条件

−k∂T∂x∣x=L=hskin (Tskin−Tcore)-k\frac{\partial T}{\partial x}\Big|_{x=L}=h_{skin}\,(T_{skin}-T_{core})

这是一个双曲型边界 + 抛物型内域的初边值问题。由于各层物性在界面处突变，且第Ⅳ层存在辐射换热，解析求解极其困难，因此采用数值解。

3.2 有限体积离散与 Thomas 求解

把控制方程在区间 [xi,xi+1][x_i,x_{i+1}] 上积分，用界面热流守恒得到半离散格式。对空间导数采用中心差分，对每个节点 i 有

ρiciTin+1−TinΔt=ki+1/2(Ti+1n+1−Tin+1)Δxi2−ki−1/2(Tin+1−Ti−1n+1)Δxi−12\rho_i c_i \frac{T_i^{n+1}-T_i^n}{\Delta t}=\frac{k_{i+1/2}(T_{i+1}^{n+1}-T_i^{n+1})}{\Delta x_i^2}-\frac{k_{i-1/2}(T_i^{n+1}-T_{i-1}^{n+1})}{\Delta x_{i-1}^2}

其中界面导热系数 ki+1/2k_{i+1/2} 在层贴合网格下直接取该层材料值，无需调和平均。整理后每个时间步得到三对角线性方程组 ATn+1=bA\mathbf{T}^{n+1}=\mathbf{b}，用 Thomas 算法（追赶法） 以 O(N) 复杂度求解。

时间推进采用隐式后向欧拉格式：未知数全取第 n+1 步。隐式格式是无条件稳定的，因此可以用较大的 Δt=30 s\Delta t=30\ \mathrm{s} 而不发散——这比显式格式（受傅里叶数限制、步长需小到秒级）快两三个数量级，使大规模参数扫描与蒙特卡洛成为可能。代价是每个时间步需解一次三对角系统，但 N=60 时成本可忽略。

3.3 层贴合有限体积网格

早期尝试用均匀网格剖分整层，但发现层界面常落在单元内部，导致界面导热率需要调和平均、且皮肤温度出现 ±0.1°C 的非单调"抖动"——这正是二分法寻优锁到错误凹坑、辨识结果发散的根源。本篇改用层贴合（layer-conforming）有限体积网格：每层内部均匀剖分，使每个层界面恒为节点。这样每个控制体单元完全属于一种材料，导热率无需调和平均，界面条件自动精确满足，稳态一致性误差降至 4.69\times10^{-7}\,^\circ\mathrm{C}。

3.4 空气隙的辐射修正

第Ⅳ层是空气隙，厚度仅几毫米，实际换热包含分子导热与两壁面间的辐射。把辐射项线性化（围绕平均温度展开）后为等效导热系数

keff=k4+HRAD⋅L4,HRAD=6.5 W/(m2⋅K)k_{eff}=k_4+H_{RAD}\cdot L_4,\qquad H_{RAD}=6.5\ \mathrm{W/(m^2\cdot K)}

其有效热阻 Rair=L4/keffR_{air}=L_4/k_{eff} 随 L4L_4 增大而增大，但存在上界 1/HRAD=0.1538 m2⋅K/W1/H_{RAD}=0.1538\ \mathrm{m^2\cdot K/W}——这意味着单纯加厚空气隙的增益会饱和，这是后续"为什么不能只靠空气隙"的物理基础。线性化带来的近似误差在本题温度范围内可忽略，因为辐射热流本身仅占空气隙换热的一小部分。

图1

4 模型验证

数值模型必须经得住检验，才能用于设计。本篇用两种独立方式验证。

4.1 半无限体 erfc 解析解对照

对单一均质半无限体（表面突升到 TsurfT_{surf}、初始 T0T_0），有精确解

T(x,t)=T0+(Tsurf−T0) erfc ⁣(x2αt),α=kρcT(x,t)=T_0+(T_{surf}-T_0)\,\mathrm{erfc}\!\left(\frac{x}{2\sqrt{\alpha t}}\right),\quad \alpha=\frac{k}{\rho c}

用同样的隐式 FD + Thomas 求解器复算该问题，在多个 (x,t)(x,t) 检查点上与 erfc 解对比，得到最大绝对误差 0.3833°C。该误差主要来自有限域截断（半无限体被截断为 0.2 m 有限域）与离散步长，量级合理，说明求解器核心正确。这一验证独立于多层结构，专门检验"瞬态热传导 + 隐式时间积分"这一最基础的数值内核。

图2

4.2 稳态串联热阻闭式解对照

当时间足够长，系统趋于稳态，四层等价于串联热阻

Rtot=∑jLjkj,Tskinss=Tenv−Tenv−Tcore1+hskinRtot (hskinRtot)R_{tot}=\sum_j \frac{L_j}{k_j},\qquad T_{skin}^{ss}=T_{env}-\frac{T_{env}-T_{core}}{1+h_{skin}R_{tot}}\,(h_{skin}R_{tot})

对 6 组不同 (L2,L4,Tenv)(L_2,L_4,T_{env}) 工况，把数值稳态皮肤温与上式闭式解对照，最大误差仅 4.69\times10^{-7}\,^\circ\mathrm{C}——层贴合网格消除了界面错位后，稳态一致性近乎完美。这条验证检验的是"多层串联 + 界面热流连续"在稳态极限下是否正确。

图3
图4

5 各层热阻与热容分解

把四层的热阻 Rj=Lj/kjR_j=L_j/k_j 与热容 Cj=ρjcjLjC_j=\rho_j c_j L_j 分别求和，并计算每层占比，得到关键结构认知：

  • 热阻：层Ⅲ占 42.15%、层Ⅳ占 45.45%，二者合计近 88%——它们是阻挡热流的主要屏障；层Ⅱ虽为设计变量，但其固定基准厚度下热阻占比仅 8.5%；
  • 热容：层Ⅱ占 93.82%，远高于其余层（层Ⅰ 2.1%、层Ⅲ 4.0%、层Ⅳ 0.06%）——隔热层几乎独占了整衣的"热惯性"。

这一分解预示了后续优化中一个深刻现象：在长时间工况（接近稳态）下，层Ⅱ靠"热阻"起作用；在短时间工况（远未达稳态）下，层Ⅱ靠"热容"起作用。集总参数法给出时间常数 τ≈RtotCtot=2197 s\tau\approx R_{tot}C_{tot}=2197\ \mathrm{s}，但数值实测的 63% 温升时间 t63≈200 st_{63}\approx 200\ \mathrm{s}——集总法严重高估了响应速度（差一个数量级），进一步说明必须做分布式的瞬态求解，不能套用简单的一阶惯性环节。

图5

6 空气隙辐射与增益饱和

空气隙的有效热阻 Rair=L4/keffR_{air}=L_4/k_{eff} 随间隙增厚单调上升，但上界 1/HRAD=0.15381/H_{RAD}=0.1538。在可行域 L4∈[0.6,6.4] mmL_4\in[0.6,6.4]\ \mathrm{mm} 内，每增厚 1mm 的热阻增量从 0.0275 一路降到 0.00557（单位 10−3 m2⋅K/W10^{-3}\,\mathrm{m^2\cdot K/W} 每 mm），衰减明显。这说明空气隙"越厚越不划算"，为问题三的权衡提供物理依据：当空气隙已经较厚时，继续加厚它的边际增益远低于加厚隔热层，因此最优设计必然落在两层之间的某个折中。

图6

7 基于 75°C 假人实验的参数辨识

题设之外，本篇构造了一组合成的 75°C 假人实验：固定基线设计 (L2,L4)=(6.0,5.5) mm(L_2,L_4)=(6.0,5.5)\ \mathrm{mm}，用正演模型生成皮肤温度时间序列，并叠加高斯噪声（标准差约 0.5°C，模拟探测误差与个体波动），得到 181 个观测点。我们用这组"数据"反演两个未知参数：皮肤等效换热系数 hskinh_{skin} 与隔热层导热系数 k2k_2。

辨识方法：在 (h,k2)(h,k_2) 平面上做网格扫描，对每个候选参数用正演模型算出皮肤温曲线，与观测序列比较残差平方和，取最小者。为避免陷入局部极小，网格覆盖 h∈[10,18]h\in[10,18]、k2∈[0.25,0.49]k_2\in[0.25,0.49] 的较宽范围，并对噪声实现做固定随机种子以保证可复现。结果为

h^skin=14.000 W/(m2⋅K),k^2=0.3640 W/(m⋅K)\hat h_{skin}=14.000\ \mathrm{W/(m^2\cdot K)},\qquad \hat k_2=0.3640\ \mathrm{W/(m\cdot K)}

而真值分别为 14.0 与 0.370，残差 RMSE = 0.350°C，决定系数 R2=0.9520R^2=0.9520。模型曲线与含噪实测高度吻合（图7），说明正演模型既能"正算预测"也能"反演参数"，可信度得到闭环验证。值得注意，k2 的估计值 0.3640 略低于真值 0.370，这是噪声与模型近似带来的微小偏差，但在工程容许范围内——辨识的目的不是完美还原，而是确认模型对关键参数"敏感且可识别"。

图7
图8

8 结论

本篇建立了四层服装的一维非稳态热传导模型，核心工程选择是：层贴合有限体积网格 + 隐式后向欧拉 + Thomas 求解 + 空气隙辐射修正。模型经半无限体 erfc 解（误差 0.3833°C）与稳态闭式解（误差 < 5×10⁻⁷°C）双重验证；热阻/热容分解揭示"层Ⅲ+Ⅳ主导热阻、层Ⅱ主导热容"的结构特征，并指出集总参数法会严重高估响应速度；最后用合成 75°C 假人实验把 hskinh_{skin}、k2 辨识到真值附近（R2=0.952R^2=0.952）。模型已具备支撑后续优化与稳健设计的可信度，下一篇将在此基础上完成确定性厚度优化。

附录：核心 Python 实现（可复现上述数字）

import sys, os
sys.path.insert(0, os.path.join(os.path.dirname(__file__), "..", "..", "tools"))
import gen_data as GD

R = GD.gen_2018a()                       # 确定性生成全部数据（唯一真源）
print("验证(erfc)最大误差 =", round(R["val_err"], 4), "°C")
print("稳态一致性最大误差 =", R["steady_err"], "°C")
print("层热阻占比 R_frac =", [round(x, 4) for x in R["R_frac"]])
print("层热容占比 C_frac =", [round(x, 4) for x in R["C_frac"]])
print("集总 tau =", round(R["tau_base"], 1), "s | 实测 t63 =", round(R["ss_base"], 1), "s 量级")
print("空气隙热阻上界 1/H_RAD =", round(1.0 / R["H_RAD"], 4))
print("辨识 h_hat=%.3f (真14.0)  k2_hat=%.4f (真0.370)" % (R["h_hat"], R["k2_hat"]))
print("辨识 RMSE=%.3f  R2=%.4f" % (R["fit_rmse"], R["fit_r2"]))
print("75°C 实验 CSV 已写出：", R["csv"])