一维非定常热传导的变量分离法
一维非定常热传导变量分离法的理论基础
概述 — 变量分离法是什么
变量分离法用于热传导的哪些问题呢？
这是求解均匀截面平板、圆柱、球体非定常温度分布的精确方法。通过将偏微分方程分解为"仅关于空间的函数"与"仅关于时间的函数"的乘积，利用特征函数展开获得无限级数的精确解。
精确解是指像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。相差两个数量级以上。这就是为什么铜块瞬间升温，但混凝土墙需要数小时的原因。
为了统一处理，我们引入无量纲温度和无量纲数：
其中 $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在里面好像很难……
这是超越方程，解析方法无法求解。需要用Newton法或Brent法数值求根。教科书附表中列出了按毕渥数计算的 $\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%以内。
就只加一项？高次项那么快就消失了吗？
让我们来具体算一下。第一类BC时，$\zeta_1 = \pi/2 \approx 1.571$，$\zeta_2 = 3\pi/2 \approx 4.712$。当Fo=0.2时：
- 第1项：$\exp(-\zeta_1^2 \times 0.2) = \exp(-0.493) \approx 0.611$ — 还很大
- 第2项：$\exp(-\zeta_2^2 \times 0.2) = \exp(-4.44) \approx 0.012$ — 只有1%
- 第3项：$\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) 一行代码就搞定。
傅里叶的执着与拿破仑
确立变量分离法的约瑟夫·傅里叶（1768-1830）曾随拿破仑远征埃及。他亲身体验了撒哈拉沙漠白天的酷热和夜晚的极度寒冷，对温度变化产生了浓厚兴趣。1822年发表的名著《热的解析理论》首次提出变量分离法和傅里叶级数，但遭到拉格朗日和拉普拉斯等同时代大师的质疑——"任意函数能表示为三角函数之和？不可能！"最后证明傅里叶正确无误，这场论争成为现代泛函分析的起点。
一维非定常热传导变量分离法的数值计算方法
解析解的数值计算步骤
在计算机上实际算出变量分离法的解，步骤是什么？
分4个步骤：
- 计算特征值：用Newton法求解 $\zeta_n \tan(\zeta_n) = \text{Bi}$。第 $n$ 个根必在区间 $((n-1)\pi, \; (n-0.5)\pi)$ 内，用Brent法二分最保险
- 计算展开系数：对每个 $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}$）时停止求和
用Python怎么写？
用SciPy大概是这样：
from scipy.optimize import brentq
import numpy as np
def eigenvalues_plate(Bi, N=50):
"""计算平板、第三类BC的特征值共N个"""
zeta = []
for n in range(1, N+1):
a = (n - 1) * np.pi + 1e-12
b = (n - 0.5) * np.pi - 1e-12
f = lambda z: z * np.tan(z) - Bi
zeta.append(brentq(f, a, b))
return np.array(zeta)
def theta_plate(x_L, Fo, Bi, N=50):
"""计算无量纲温度 theta(x/L, Fo, Bi)"""
zn = eigenvalues_plate(Bi, N)
Cn = 4 * np.sin(zn) / (2*zn + np.sin(2*zn))
return np.sum(Cn * np.cos(zn * x_L)
* np.exp(-zn**2 * Fo))与FEM的比较验证
用这个解析解怎么验证FEM的结果？
标志性的例子是NAFEMS T2基准。0.1m厚钢板（$\alpha = 3.5 \times 10^{-6}$ m²/s），一面固定100°C，另一面绝热。32秒后 $x = 0.08$ m 处的温度参考值是 $T = 36.60$ °C。
FEM的结果与36.60°C差多少算不合格？
通常要求误差在1%以内（±0.37°C）。但这是网格和时间步长都足够细致时的标准。FEM的验收流程是：
- 粗网格×大时间步做初次计算
- 网格加倍细分再算——检查温度变化
- 时间步减半再算——检查温度变化
- 重复直到两者的变化都在1%以内
- 最终结果与变量分离法的精确解对比
时间积分方案的影响
FEM的时间积分方法对精度有影响吗？
影响很大。看一下主要的三种方案：
| 方案 | $\theta_{\text{method}}$ | 精度 | 稳定性 | 备注 |
|---|---|---|---|---|
| 前向欧拉（显式） | 0 | $O(\Delta t)$ | 条件稳定：$\Delta t < \Delta x^2 / (2\alpha)$ | 时间步长约束很严 |
| Crank-Nicolson | 0.5 | $O(\Delta t^2)$ | 无条件稳定 | 可能出现数值振荡 |
| 后向欧拉（隐式） | 1 | $O(\Delta t)$ | 无条件稳定 | 安全但扩散可能过度 |
Crank-Nicolson精度最好，但有"数值振荡"？很危险啊……
比如淬火最初瞬间温度剧变，Crank-Nicolson可能会让温度出现不物理的负值。为避免这点，Ansys等商用软件常把 $\theta = 2/3$（Galerkin法）作为默认值。与变量分离法精确解的对比，可以快速发现这类数值人工误差。
混凝土墙的日射响应手算
厚200mm的混凝土墙（$\alpha = 7 \times 10^{-7}$ m²/s）突然受到日射升温20°C时，1小时后的傅里叶数 Fo = $7 \times 10^{-7} \times 3600 / 0.1^2 = 0.252$。由于 Fo > 0.2，可以用单项近似，快速估算内侧温升只有几摄氏度。这样的手算检验在JIS A 1470隔热性能测试中广泛应用，是"跑FEM前的合理性检查"。
一维非定常热传导变量分离法的实务应用
淬火（急冷）解析的应用
淬火处理中，变量分离法怎么具体用？
典型案例："直径50mm的钢棒从850°C投入水中急冷"。用无限圆柱模型，水的热传系数 $h = 5000$ W/(m²·K)，钢的 $k = 50$ W/(m·K)，则：
- $\text{Bi} = h \cdot r_0 / k = 5000 \times 0.025 / 50 = 2.5$
- 求解特征值方程 $\zeta_1 J_1(\zeta_1)/J_0(\zeta_1) = 2.5$，得 $\zeta_1 \approx 1.814$
- 中心温度降到300°C时，$\theta_0 = (300 - 20)/(850 - 20) = 0.337$
- 反推傅里叶数：$0.337 = C_1 \exp(-\zeta_1^2 \text{Fo})$ → $\text{Fo} \approx 0.152$ → $t = \text{Fo} \cdot r_0^2 / \alpha \approx 4.8$ 秒
不到5秒中心就冷到300°C。将此冷却速率与连续冷却转变（CCT）图结合，就能预测淬火硬度。
但水中急冷时，沸腾会导致 $h$ 变化，对吧？
说得好。实际淬火经历膜沸腾→核沸腾→对流冷却三个阶段，$h$ 剧烈变化。变量分离法只能处理 $h$ 恒定的情况。所以它的角色是"粗估"——"大概5秒冷下来"——可信的。精细分析必须用FEM，让 $h(T)$ 作为温度函数输入。
隔热材·多层壁的设计检讨
多层壁（隔热材夹在中间那种）也能用吗？
严格来说不行。多层材料物性不同，简单的变量分离法失效。但实务上有些对策：
- 等效物性法：多层壁等效为单层（计算等效 $\alpha$）再用变量分离法
- 各层分别分离：每层用自己的特征函数，层间施加温度和热流连续条件——数学很复杂
- 实际推荐：层数≥3时直接用1D-FEM更省力
单层的问题——比如"50mm玻璃棉的一面突然加热，另一面升温需多久？"——可以直接用变量分离法。
商用工具中的实现步骤
在CAE软件中设置变量分离法的验证问题时，要注意什么？
有几个常见陷阱：
| 项目 | 注意事项 |
|---|---|
| 网格 | 一维问题也得用2D/3D单元的工具很多。平板就沿厚度剪裁、其他方向用对称BC固定 |
| 初始条件 | $T_i$ 要均匀赋给全部节点。Ansys例：IC,ALL,TEMP,850 |
| 时间步 | Fo=0.2对应的时刻前后要重点监测。Abaqus例：*HEAT TRANSFER, DELTMX=5 限制最大温变 |
| 物性温度相关性 | 验证目的一定要关闭温度相关物性。COMSOL用"Materials: custom"确保常数 |
| 输出 | 导出节点温度时间序列。事先确定要对比的位置（中心、表面、特定 $x$） |
一维非定常热传导变量分离法的软件比较
工具别非定常热传导响应
不同CAE软件中，变量分离法的验证问题有差异吗？
看看各工具的非定常热传导设置和与变量分离法验证的相容性：
| 工具 | 分析指令 | 时间积分 | 验证相容性 |
|---|---|---|---|
| Ansys Mechanical | ANTYPE,TRANS | 广义梯形（默认$\theta=2/3$） | APDL宏可自动与解析解对比 |
| Abaqus | *HEAT TRANSFER | 后向欧拉（默认） | Python脚本+odb后处理自动化 |
| COMSOL | Heat Transfer → Time Dependent | BDF（变步长） | 可在"Analytic"功能中直接输入解析式并绘制差值 |
| OpenFOAM | laplacianFoam | Euler / Crank-Nicolson | Python后处理脚本。纯固体只需laplacianFoam |
COMSOL的"Analytic"功能听起来很方便。
是的。COMSOL可以在GUI里直接输入解析式，非常便捷。不过要注意特征值计算的精度——最好自己也验算一遍。
开源工具
商用软件以外还有什么工具能算变量分离法的解？
- Python + SciPy：前面展示的脚本，从特征值求解到温度分布绘图全自动。最灵活
- Wolfram Mathematica：用
DSolve命令可符号求解，学术论文输出很漂亮 - Octave / MATLAB：用
fzero求特征值，besselj处理圆柱问题 - hteqn（Python库）：一行函数调用就得到平板、圆柱、球的解。内置与海斯勒图对比功能
10行Python代码写出"变量分离法计算器"
装好NumPy和SciPy，变量分离法的一维非定常热传导解只需10行搞定。特征值用Brent法求50个，级数求和。代码写好后，"毕渥数3的圆柱，中心温度衰到一半需多久？"这样的问题一瞬间就有答案。工作台常驻这个脚本，比启动FEM快得多。
一维非定常热传导变量分离法的前沿研究
PINN与变量分离法的融合
200年前的方法跟最新AI技术有关联？
有！物理信息神经网络（PINN）的最新研究表明，把变量分离法的"时间×空间"结构硬编码到网络架构中，效果大幅改善。
- 空间子网络：学习接近特征函数的基函数
- 时间子网络：学习指数衰减规律
- 结果：比传统PINN收敛快得多，数据需求少一个数量级，且解的物理正确性有保证
两个世纪前傅里叶发现的数学结构，今天在深度学习中作为"归纳偏置"重生。
非傅里叶传导的扩展
"非傅里叶传导"是什么？变量分离法用不了？
傅里叶定律假设"热流对温度梯度瞬时响应"。但在超短脉冲激光加工（皮秒-飞秒）或极低温（<1K）下，热传播有有限速度。这种情况用Cattaneo-Vernotte方程描述：
时间二阶导数……像波动方程。
完全一致。确实出现了"热波（thermal wave）"。惊人的是，这个方程也能用变量分离法处理——特征值方程形式改变但求解思路相同。不过在实际工程中，这是个非常小众的应用（只有半导体飞秒激光加工等才用）。常规热分析还是傅里叶定律。
一维非定常热传导变量分离法的故障排除
级数不收敛的情况
级数一直在振荡，不收敛是怎回事？
这是Fo极小（极初期）的典型症状。解决办法：
| 现象 | 原因 | 对策 |
|---|---|---|
| 100项都不收敛 | Fo < 0.01（极初期） | 改用半无限体的误差函数解 |
| 表面温度振荡（吉布斯现象） | 初始条件不连续（阶梯型） | 取500项以上，或避开表面区域（用 $x/L \lesssim 0.95$） |
| 特征值求解失败 | 搜索区间设置错 | 严格遵守 $\zeta_n \in ((n-1)\pi, (n-0.5)\pi)$。先用二分法，再切换Newton法 |
| 中心温度对但表面错 | 毕渥数输入错（分不清半厚全厚） | 变量分离法中 $L$ = 半厚。与FEM的输入对标 |
FEM与解析解不一致的情况
FEM结果与变量分离法差数%。谁对谁错？
冷静排查。对照单清单：
- 变量分离法侧检查：
- 特征值计算够20个以上吗？
- $L$ 定义是半厚还是全厚？（教科书有差异）
- 展开系数公式对应该边界条件吗？
- FEM侧检查：
- 物性设定为温度无关吗？温度函数表没混进去？
- 网格收敛验证过了吗？特别是表面附近（高温度梯度区）
- 时间步够小吗？自动时间步的最大值限制合理吗？
- 对称面设为绝热吗？（绝热≠自由面）
- 如果还不一致：从最简单的情况（第一类BC、均匀初温）开始验证，复杂BC留后面
明白了。最简单的条件先试，复杂的后来？
正是。调试的铁律是"一次改一个变量"。变量分离法有"确定的正确答案"这个最强武器，要充分利用。
更详细
错误