一维非稳态热传导变量分离法
理论与物理
概述 — 什么是分离变量法
分离变量法用于解决热传导中的哪些问题？
这是一种精确求解均匀截面平板、圆柱、球体非稳态温度分布的方法。它将偏微分方程分离为“仅空间函数”和“仅时间函数”的乘积，通过特征函数展开得到无穷级数解。
您说的精确解，是指不像FEM那样近似，而是能得到完全准确的答案吗？
正是如此。但前提是“物性值不依赖于温度”、“形状简单（平板、圆柱、球体）”、“边界条件是线性的”。只要满足这些条件，得到的解在数学上是精确的。因此它作为FEM的验证基准极为重要。
具体在什么场景下使用呢？
举三个典型用途：
- 淬火（急冷）处理的内部温度预测 — 钢材投入水中后，中心温度何时通过马氏体相变点。这可以用分离变量法以秒为单位进行预测。
- FEM解的验证（V&V） — NAFEMS T2基准问题就是以分离变量法的精确解作为参考值的。
- 设计初期的概算 — 在运行FEM之前，用于工计算估算混凝土墙的日照响应、电子基板突然通电时的温升。
哦，就像是用来确认FEM是否正确的“答案核对表”一样的东西呢。
正是如此。不了解分离变量法的解析解而只依赖FEM，就像没有计算器却相信心算结果一样。
控制方程与无量纲化
请先从基本方程开始讲解。
一维热传导方程如下。无内热源、物性值恒定的情况：
其中 $\alpha = k / (\rho c_p)$ 是热扩散率 [m²/s]，$k$ 是热导率，$\rho$ 是密度，$c_p$ 是比热容。
$\alpha$ 的值在不同材料间大概有多大差异？
粗略来说，铜的 $\alpha \approx 1.1 \times 10^{-4}$ m²/s，钢的 $\alpha \approx 1.2 \times 10^{-5}$ m²/s，混凝土的 $\alpha \approx 7 \times 10^{-7}$ m²/s。存在两个数量级以上的差异。铜块能瞬间变热，而混凝土墙则需要好几个小时，对吧。
为了统一处理分析，引入无量纲温度 $\theta$ 和无量纲数：
其中 $T_i$ 是初始温度，$T_\infty$ 是周围流体温度，$L$ 是特征长度（平板半厚度），$h$ 是对流传热系数。$\theta=1$ 对应初始状态，$\theta=0$ 对应稳态（与周围温度一致）。
分离变量的步骤
终于到“分离变量”的核心内容了。具体怎么做呢？
将温度分离为“仅空间函数 $X(x)$”和“仅时间函数 $\Gamma(t)$”的乘积：
将其代入热传导方程得：
左边只是 $t$ 的函数，右边只是 $x$ 的函数，却相等……这意味着两边必须都是常数，否则说不通，是吗？
没错！这就是分离变量法的关键。这个常数 $-\lambda^2$ 称为分离常数。加上负号是为了反映温度随时间衰减（不发散）的物理现象。
由此得到两个常微分方程：
时间部分就是简单的指数衰减呢。空间部分是三角函数……有点像波动方程？
很好的观察。数学上确实是相同形式的特征值问题。区别在于波动会持续振动，而热传导由于 $e^{-\alpha\lambda^2 t}$ 的作用会随时间衰减并趋于稳态。
边界条件与特征值方程
$\lambda$ 的值如何确定？
施加边界条件后，$\lambda$ 可取的值被限制为离散的特征值 $\lambda_n$。令 $\zeta_n = \lambda_n L$ 进行无量纲化，则根据边界条件的类型，特征值方程就确定了。
对于厚度 $2L$ 的平板，在中心 $x=0$ 施加对称条件，在两面 $x=\pm L$ 施加各种边界条件的情况：
| 边界条件 | 特征值方程 | 特征函数 $X_n$ |
|---|---|---|
| 第一类（表面温度恒定） | $\cos(\zeta_n) = 0$ 即 $\zeta_n = (2n-1)\pi/2$ | $\cos(\zeta_n x/L)$ |
| 第二类（表面热流恒定） | $\sin(\zeta_n) = 0$ 即 $\zeta_n = n\pi$ | $\cos(\zeta_n x/L)$ |
| 第三类（对流: $-k \partial T/\partial x = h(T-T_\infty)$） | $\zeta_n \tan(\zeta_n) = \text{Bi}$ | $\cos(\zeta_n x/L)$ |
第三类的 $\zeta_n \tan(\zeta_n) = \text{Bi}$，这个能解吗？有 tan 在里面，手算似乎不可能……
这是超越方程，无法解析求解。需要用牛顿法或布伦特法数值求根。你见过教科书附录里按毕渥数列出的 $\zeta_1$ 值表吧？那就是事先计算好列成表的。例如 Bi=1 时，$\zeta_1 \approx 0.8603$。
利用特征函数的正交性，展开系数 $C_n$ 可由初始条件唯一确定。平板初始温度均匀为 $T_i$ 时：
最终的解由以下无穷级数表示：
傅里叶数与单项近似
无穷级数，要计算多少项才够呢？
这正是此方法美妙之处。看 $\exp(-\zeta_n^2 \text{Fo})$ 这一项。$n$ 越大，$\zeta_n$ 也越大，因此指数衰减得越快。当 Fo > 0.2 时，仅第一项（$n=1$）就能达到 0.1% 以内的精度。
诶，只用一项？高次项消失得那么快吗？
我们来具体计算一下。设第一类边界条件下 $\zeta_1 = \pi/2 \approx 1.571$，$\zeta_2 = 3\pi/2 \approx 4.712$。Fo=0.2 时：
- 第一项: $\exp(-\zeta_1^2 \times 0.2) = \exp(-0.493) \approx 0.611$ — 仍然较大
- 第二项: $\exp(-\zeta_2^2 \times 0.2) = \exp(-4.44) \approx 0.012$ — 只剩 1% 左右
- 第三项: $\exp(-\zeta_3^2 \times 0.2) = \exp(-12.3) \approx 4.5 \times 10^{-6}$ — 实际上为零
所以“单项近似（one-term approximation）”就足够了。海斯勒图正是将这个单项近似图表化的结果。
应用单项近似后的中心温度 $\theta_0$（$x=0$）公式：
反过来，Fo 很小的时候，也就是加热/冷却的初始阶段，就需要很多项了，是吗？
是的。Fo < 0.05 左右的极初期，有时需要100项以上。这种情况下，使用半无限大物体的解（误差函数解）比用分离变量法更高效。这将在另一个专题中详细讨论。
向圆柱和球的扩展
除了平板，还能用于其他形状吗？
无限长圆柱和球体也可以使用相同的方法。只是坐标系和特征函数会改变。
| 形状 | 坐标 | 特征函数 | 特征值方程（第三类BC） |
|---|---|---|---|
| 平板 | 直角坐标 $x$ | $\cos(\zeta_n x/L)$ | $\zeta_n \tan(\zeta_n) = \text{Bi}$ |
| 无限长圆柱 | 圆柱坐标 $r$ | $J_0(\zeta_n r/r_0)$ | $\zeta_n J_1(\zeta_n)/J_0(\zeta_n) = \text{Bi}$ |
| 球 | 球坐标 $r$ | $\sin(\zeta_n r/r_0)/(\zeta_n r/r_0)$ | $1 - \zeta_n \cot(\zeta_n) = \text{Bi}$ |
$J_0$ 是什么？没见过这个符号。
是第一类零阶贝塞尔函数。在圆柱坐标下写出拉普拉斯算子后，代替平板的 $\cos$ 或 $\sin$ 出现的就是贝塞尔函数。在 Python 中，用 scipy.special.j0(x) 就能直接计算，所以在实际工作中不用怕。
各项的物理意义
- 蓄热项 $\rho c_p \partial T/\partial t$：单位体积的热能储存率。铁锅难热难冷（$\rho c_p$ 大），而铝锅易热易冷。这个乘积越大，温度变化越缓慢。
- 热传导项 $k \partial^2 T/\partial x^2$：基于傅里叶定律的热输运。与温度分布的曲率（凹凸）成比例地移动热量。金属勺子放入热锅后连把手都会变热，就是因为 $k$ 大。
- 热扩散率 $\alpha = k/(\rho c_p)$：表示“温度均匀化的速度有多快”。越大则温度扩散越快。铜 ($1.1 \times 10^{-4}$) 和木材 ($1.3 \times 10^{-7}$) 有3个数量级的差异。
- 傅里叶数 $\text{Fo} = \alpha t / L^2$：无量纲时间。热扩散距离与物体尺寸之比的平方。Fo=1 大致对应“热已遍布整个物体”的状态。
- 毕渥数 $\text{Bi} = hL/k$：表面对流热阻与内部传导热阻之比。Bi < 0.1 则内部温度几乎均匀（可使用集总热容法）。Bi > 100 则表面温度立即变为 $T_\infty$（第一类边界条件近似）。
假设条件与适用范围
- 物性值不依赖于温度：$k$, $\rho$, $c_p$ 恒定。大温差（数百度以上）时不能忽略温度依赖性，分离变量法不适用。需要FEM。
- 各向同性材料：热导率不依赖于方向。对于CFRP等复合材料，需要考虑各向异性的其他方法。
- 无内热源：存在焦耳热或化学反应热时，有时可以通过与稳态解叠加来处理。
- 线性边界条件：无法直接应用于包含辐射（与 $T^4$ 成比例）的边界条件。需要线性化近似。
- 一维问题：对于二维以上的问题，如果是直角坐标系，可以扩展为各方向解的乘积（product solution）。
量纲分析与单位制
| 变量 | SI单位 | 典型值・注意事项 |
|---|---|---|
| 温度 $T$ | K 或 °C | 温差计算用K或°C结果相同。辐射计算必须用绝对温度 |
| 热导率 $k$ | W/(m·K) | 铜: 400、钢: 50、混凝土: 1.4、空气: 0.026 |
| 热扩散率 $\alpha$ | m²/s | 铜: $1.1\times10^{-4}$、钢: $1.2\times10^{-5}$、水: $1.4\times10^{-7}$ |
| 对流传热系数 $h$ | W/(m²·K) | 自然对流: 5〜25、强制对流: 25〜250、水中急冷: 1,000〜10,000 |
| 傅里叶数 Fo | 无量纲 | Fo > 0.2 时可应用单项近似。实际工作中最频繁使用的判定标准 |
| 毕渥数 Bi | 无量纲 | Bi < 0.1 则用集总热容法。0.1 < Bi < 100 是分离变量法的用武之地 |
傅里叶的热情与拿破仑
确立了分离变量法的约瑟夫·傅里叶（1768-1830）曾随拿破仑远征埃及。据说他亲身体验了沙漠的酷热与夜晚的极端寒冷，从而痴迷于“温度为何如此变化”这一问题。他在1822年的名著《热的解析理论（Theorie analytique de la chaleur）》中发表了分离变量法和傅里叶级数，但当时的拉格朗日和拉普拉斯并不相信“任意函数都能用三角函数的和表示”这一主张。最终证明傅里叶是正确的，而这场争论成为了现代函数分析学的起点。
数值解法与实现
解析解的数值计算步骤
要实际用计算机计算分离变量法的解，应该按什么步骤进行？
步骤分为4步：
- 特征值的计算: 用牛顿法解 $\zeta_n \tan(\zeta_n) = \text{Bi}$。第 $n$ 个根必定存在于区间 $((n-1)\pi, \; (n-0.5)\pi)$ 内，所以用布伦特法夹逼更可靠。
- 展开系数的计算: 对每个 $n$ 计算 $C_n = 4\sin(\zeta_n) / (2\zeta_n + \sin(2\zeta_n))$。
- 级数求和: 对于所需的 $x$ 和 $t$（Fo）计算 $\theta = \sum C_n \cos(\zeta_n x/L) \exp(-\zeta_n^2 \text{Fo})$。
- 截断判定: 当 $|C_n \exp(-\zeta_n^2 \text{Fo})| < \epsilon$（例如 $10^{-10}$）时，截断级数。
相关主题
なった
詳しく
報告