多次元非定常熱伝導
理论与物理
概述 — 为何需要多维分析
老师，多维非稳态热传导问题只能用有限元法求解吗？我在教科书上看过一维的Heisler图，但实际零件都是三维的吧…。
问得好。实际上，如果正交坐标系中各方向的边界条件相互独立，就可以使用“乘积解（product solution）”。例如，短角柱（类似砖块的形状）的温度可以这样求得：
即，只需将各方向的一维解相乘，就能得到三维非稳态温度场。即使没有有限元法，也能通过手算获得三维非稳态温度场。
诶，只用乘法就能知道三维温度分布吗！？这是什么原理呢？
关键在于变量分离法。如果多维热传导方程在各轴方向的边界条件是可分离的，那么它就可以分解为各方向的一维问题。从数学上讲，这意味着偏微分方程可以归结为各变量常微分方程的乘积。例如，考虑汽车锻件（短圆柱形状）的淬火冷却。可以将侧面视为圆柱，上下表面视为平板，然后将各自的一维解相乘即可。
原来如此！那么，这是否也可以用于验证有限元法的结果呢？
正是如此。乘积解作为有限元法的基准（验证解）非常有用。可以将网格或时间步长的合理性与精确的解析解进行定量比较。在实际工作中，是否了解这一点，对分析结果的可信度影响巨大。
控制方程
那么，首先请告诉我多维非稳态热传导的基本方程。
起点是傅里叶热传导方程的多维版本。对于各向同性、均匀且无内热源的介质：
其中 $\alpha = k/(\rho c_p)$ 是热扩散率 [m²/s]。这个方程是抛物型偏微分方程，给定初始条件 $T(\mathbf{x}, 0) = T_i$ 和边界条件后，解就唯一确定了。
在圆柱坐标系中会变成什么样呢？像发动机气缸或圆筒形零件有很多吧。
在圆柱坐标系 $(r, \phi, z)$ 中：
对于轴对称（$\phi$ 方向均匀）且长度有限的情况，就变成了 $r$ 方向和 $z$ 方向的二维问题。这正是乘积解大显身手的地方。
我们来整理一下无量纲温度 $\theta^*$、毕渥数 $\text{Bi}$、傅里叶数 $\text{Fo}$ 的定义：
其中 $h$ 是对流换热系数，$L_c$ 是特征长度（平板为半厚，圆柱为半径），$T_\infty$ 是周围流体温度，$T_i$ 是初始温度。
乘积解（Product Solution）
请再详细解释一下乘积解的数学依据。为什么可以相乘呢？
核心在于变量分离的成立条件。假设解的形式为 $\theta^*(\mathbf{x},t) = X(x,t) \cdot Y(y,t) \cdot Z(z,t)$，代入控制方程后得到：
如果各项仅依赖于各自的变量，那么就可以分离为各方向的一维问题。这样得到的 $X$, $Y$, $Z$ 就是各自一维非稳态热传导的解（Heisler解）。
对于具体的几何形状会怎样呢？
举三个典型例子：
1. 短角柱（Short Rectangular Bar） — 三个方向平板解的乘积
2. 短圆柱（Short Cylinder） — 圆柱解 × 平板解
3. 半无限大固体的角部（Semi-infinite Corner） — 三个方向半无限大固体解的乘积
例如，考虑电子元件树脂封装（长方体）的回流加热，就可以直接使用案例1。
假设 $x$ 方向冷却了50%，$y$ 方向也独立地冷却了50%。那么整体上就是 $0.5 \times 0.5 = 0.25$，即冷却了75%。各方向的冷却进程是独立的，它们通过乘法组合在一起，这就是乘积解的本质。
乘积解不能使用的情况是什么样的呢？
这是个重要的点。以下情况不能应用乘积解：
- 存在内热源的情况（$\dot{q} \neq 0$）— 方程变为非齐次，无法变量分离
- 物性值随温度变化的情况 — 方程变为非线性
- 边界条件在方向间耦合的情况（倾斜边界、接触热阻的空间分布等）
- 坐标系不正交的形状（L形、T形等）
这些情况需要数值解法（有限元法、有限差分法、有限体积法）。不过，许多实际工作中的验证问题都设定为可以用乘积解处理的形状，所以了解这个方法的价值很大。
叠加原理（Superposition）
乘积解是“乘法”，那么可以用“加法”来组合解吗？
是的，这就是叠加原理（superposition principle）。对于线性的偏微分方程，可以将各个解相加来构造新的解。这是与乘积解不同的另一种强大方法。
例如，即使只有某个面被快速加热而其他面绝热这种非对称边界条件，也可以通过分解为对称问题和反对称问题再相加来求解。钢材的单面淬火就属于这种模式。
乘积解和叠加原理有什么区别呢？感觉容易混淆…。
整理如下：
| 方法 | 运算 | 使用区分 |
|---|---|---|
| 乘积解 | 乘法（$\theta^* = \theta^*_1 \cdot \theta^*_2$） | 组合不同空间方向的一维解以构成多维解 |
| 叠加原理 | 加法（$T = T_1 + T_2 - T_{\text{ref}}$） | 在同一空间组合不同边界条件的解 |
可以这样记：乘积解是“空间方向的分解”，叠加原理是“边界条件的分解”。将两者结合，就能解析处理相当复杂的问题。
形状系数法（Shape Factor Method）
教科书上也提到了“形状系数”这个概念，它和非稳态有什么关系呢？
形状系数法主要用于稳态多维热传导。定义导热形状系数（conduction shape factor）$S$：
其中 $Q$ 是单位时间的热流量 [W]。例如，埋地管道向地面的散热、半导体芯片到基板的热流路径等，只需从表中查找已有解析解的标准形状的 $S$ 值，就能计算热流量。
具体有哪些形状的 $S$ 值在表里呢？
举几个代表性的：
| 形状 | 形状系数 $S$ | 条件 |
|---|---|---|
| 同心圆柱（长度 $L$） | $\dfrac{2\pi L}{\ln(r_2/r_1)}$ | $L \gg r_2$ |
| 同心球壳 | $\dfrac{4\pi r_1 r_2}{r_2 - r_1}$ | — |
| 埋地管道 | $\dfrac{2\pi L}{\cosh^{-1}(z/r)}$ | $L \gg z$, $z > r$ |
| 立方体角部 | $0.15 L$ | 三个面为等温面 |
Incropera的传热学教科书里记载了数十种形状系数。实际工作中，通常在估算阶段使用此法，然后在详细设计阶段转向有限元法。
形状系数和热阻的概念有关系吗？
正是如此。多维导热热阻可以表示为 $R_{\text{cond}} = 1/(kS)$。将其与对流热阻 $R_{\text{conv}} = 1/(hA)$ 串联或并联组合，就可以用电路类比法手算评估复杂系统的总热阻。这是电子设备散热设计中日常使用的方法。
数值解法与实现
2D/3D Heisler图的使用方法
要将Heisler图扩展到多维，具体需要遵循什么步骤呢？
以短圆柱（半径 $r_0$、半长 $L$）为例说明步骤。假设所有表面都进行对流冷却（$h$, $T_\infty$）：
- 步骤1： 对于 $r$ 方向，计算 Bi$_{\text{cyl}} = hr_0/k$、Fo$_{\text{cyl}} = \alpha t/r_0^2$，然后从无限长圆柱的Heisler图读取中心温度 $\theta^*_{\text{cyl}}(0,t)$
- 步骤2： 对于 $z$ 方向，计算 Bi$_{\text{wall}} = hL/k$、Fo$_{\text{wall}} = \alpha t/L^2$，然后从无限大平板的Heisler图读取中心面温度 $\theta^*_{\text{wall}}(0,t)$
- 步骤3： 短圆柱的中心温度由乘积求得
如果需要表面温度或任意位置的温度，则还需使用Heisler第二图（位置修正图）对每个位置进行修正。
请用数值例子让我确认一下！比如不锈钢短圆柱急冷的情况？
我们来具体算一下。SUS304短圆柱（$r_0 = 25$ mm、$L = 50$ mm = 半长）：
- $k = 14.9$ W/(m-K)、$\alpha = 3.95 \times 10^{-6}$ m²/s
- $h = 200$ W/(m²-K)、$T_i = 500$°C、$T_\infty = 25$°C
- 目标时刻：$t = 300$ s
$r$ 方向： Bi$_{\text{cyl}} = 200 \times 0.025 / 14.9 = 0.336$、Fo$_{\text{cyl}} = 3.95 \times 10^{-6} \times 300 / 0.025^2 = 1.896$
由图得 $\theta^*_{\text{cyl}}(0,t) \approx 0.23$
$z$ 方向： Bi$_{\text{wall}} = 200 \times 0.05 / 14.9 = 0.671$、Fo$_{\text{wall}} = 3.95 \times 10^{-6} \times 300 / 0.05^2 = 0.474$
由图得 $\theta^*_{\text{wall}}(0,t) \approx 0.61$
中心温度： $\theta^* = 0.23 \times 0.61 = 0.140$
$T_{\text{center}} = T_\infty + \theta^*(T_i - T_\infty) = 25 + 0.140 \times 475 = 91.6$°C
300秒后的中心温度约为92°C。如果用有限元法求解相同条件并与该值比较，就可以验证网格、时间步长的合理性。
基于FEM的多维离散化
当乘积解无法使用时，就得依赖有限元法了吧。用有限元法求解多维非稳态热传导时，关键点是什么？
有限元法公式化的关键在于空间离散化和时间离散化的分离。对弱形式（伽辽金法）应用空间离散，会得到半离散化的常微分方程组：
其中 $[\mathbf{C}]$ 是容量矩阵（热学版的质量矩阵），$[\mathbf{K}]$ 是传导矩阵（热学版的刚度矩阵），$\{F\}$ 是外部热负荷向量。
对此，在时间方向上应用欧拉法（前向/后向）或克兰克-尼科尔森法等时间积分格式。二维问题使用三角形或四边形单元，三维问题使用四面体或六面体单元。
时间积分格式的选择
时间积分格式该怎么选呢？有前向欧拉、后向欧拉等等很多种吧。
可以用广义的 $\theta$ 法统一表示（这里的 $\theta$ 是参数，与无量纲温度不同）：
| $\theta$ 值 | 格式名称 | 特点 |
|---|---|---|
| $0$ | 前向欧拉（显式法） | 条件稳定。受限于 $\Delta t \le \Delta x^2 / (2\alpha)$ |
| $1/2$ | 克兰克-尼科尔森 | 无条件稳定、二阶精度。但有振荡解的风险 |
| $2/3$ | 伽辽金法 | 无条件稳定、抑制振荡。商用软件中常用 |
| $1$ | 后向欧拉（隐式法） | 无条件稳定、一阶精度。数值扩散较大 |
实际工作中，通常选择<
相关主题
なった
詳しく
報告