斯蒂芬-玻尔兹曼定律 — 辐射传热的基础与CAE实现

斯蒂芬-玻尔兹曼定律的理论基础

概述 — 为什么是4次幂

普朗克分光辐射亮度为

对全波长 $\lambda = 0 \sim \infty$ 积分得到全辐射能 $E_b$。使用变量替换 $x = hc/(\lambda k_B T)$，积分变为

这个定积分是 $\Gamma(4)\,\zeta(4) = 6 \times \pi^4/90 = \pi^4/15$ 的确定值。因此 $E_b \propto T^4$。4次幂来自于"三维空间中光子状态密度与 $\nu^2$ 成正比"和"玻色-爱因斯坦分布的能量平均"——这是空间维度和量子统计的必然结果。

自2019年SI重新定义以来，$k_B$、$h$、$c$ 都是确定的定义值，所以 $\sigma$ 也有了精确确定的值。CAE教科书常用 $5.67 \times 10^{-8}$ 四舍五入，但知道这个精确值在精度讨论中很有帮助。

支配方程

斯蒂芬-玻尔兹曼定律的基本形式如下。

实际物体不是黑体，因此引入放射率（emissivity） $\varepsilon$：

下表列出了代表性材料的放射率。

灰体的辐射交换

这是灰体的基本方程。重要的是是 $T^4$ 的差，不是 $(T_1 - T_2)^4$。这种非线性使CAE求解变得复杂。

例如，铸坯（$T_1 = 1200°C = 1473\,\text{K}$）和周围环境（$T_2 = 30°C = 303\,\text{K}$）的情况：

同样条件下的自然对流，取 $h \approx 10$ W/(m²K) 得 $q_{\text{conv}} = 10 \times 1170 \approx 12$ kW/m²。即辐射约为对流的18倍。在1200°C的环境中，辐射绝对主导。

多面间的辐射交换与形态系数

封闭空间（Enclosure）内多面间的辐射交换用形态系数（View Factor） $F_{ij}$ 来描述。

灰体和扩散面采用辐射度法（Radiosity Method）：

其中 $J_i$ 是面 $i$ 的辐射度（辐射+反射能的总和）。这是一个N元联立方程，CAE求解器通过矩阵运算求解。形态系数有以下重要性质：

• 互易性： $A_i F_{ij} = A_j F_{ji}$
• 求和性： $\sum_{j=1}^{N} F_{ij} = 1$（封闭空间）
• 凸面： $F_{ii} = 0$（凸面看不到自己）

辐射主导条件 — 什么时候辐射成为主角

但在真空中，没有对流，所以航天器热设计中，即使常温，辐射也是唯一的散热手段。国际空间站（ISS）的散热器工作在约280K，完全依赖 $\sigma T^4$ 散热。温度低意味着需要更大的散热面积。

斯蒂芬-玻尔兹曼定律的数值计算方法

T⁴项的线性化

其中辐射热传导系数为

使用前次迭代的温度计算 $h_{\text{rad}}$，通过牛顿-拉夫逊法迭代更新。$h_{\text{rad}}$ 可以与对流的 $h_{\text{conv}}$ 相加，集成到现有的热传导系数框架中。

这个被放入NR法的切线刚度矩阵。Abaqus的*RADIATION和Ansys的SRDOPT内部使用这种方法。收敛快，但初始温度估计不好容易发散。

有限元法中的辐射边界条件

热传导方程的弱形式中加入辐射边界条件。能量方程为：

在辐射面 $\Gamma_r$ 上的边界条件：

用伽勒金法处理弱形式后，单元级别的非线性残差向量为：

其中 $\mathbf{N}$ 是形状函数矢量。切线矩阵为：

耦合求解器策略

1. 弱耦合（Sequential）: CFD求出对流系数 → 在FEM中求解热传导+辐射 → 温度反馈给CFD。实现简单但收敛慢。
2. 强耦合（Monolithic）: 流体、固体、辐射在一个方程系统中同时求解。Ansys Fluent的S2S(Surface-to-Surface)模型就是这种。收敛快但内存消耗大。
3. DO/MC模型： 有关联介质（燃烧气体等）时，用离散纵坐标法或蒙特卡洛法解辐射输运方程(RTE)。斯蒂芬-玻尔兹曼定律用作壁面边界条件。

实务经验是"先收敛不含辐射的流场，再启用辐射模型"最稳定。全部一起打开容易发散。

斯蒂芬-玻尔兹曼定律的实务应用

分析流程

1. 确定辐射面： 高温面和它"能看见"的面。检查遮挡情况。
2. 设置放射率： 检查是否有温度依赖性。氧化状态会大幅改变放射率，不能盲目参照文献值。
3. 计算形态系数： 用求解器自动计算（半立方体法等）或用解析解验证。
4. 网格考虑： 辐射面由于T⁴的非线性，温度梯度陡峭的地方（边角）要细化网格。
5. 初始温度： 给定接近实际工况的初始温度。从300K开始达到1500K，迭代次数会爆炸。
6. 分步启用： 先启用对流 → 再加辐射。不要一次全开。

案例1: 连续铸造的辐射冷却

1. 喷雾冷却区： 水喷雾的强制对流（$h = 500\sim2000$ W/(m²K)）主导。这里斯蒂芬-玻尔兹曼辐射占总冷却的10〜20%。
2. 辐射冷却区： 喷雾停止后，铸坯温度仍在900〜1100°C。这一区间 $q_{\text{rad}} = \varepsilon\sigma(T_s^4 - T_\infty^4)$ 占冷却的80%以上。
3. 辊接触区： 导向辊接触的传导。接触面积小但局部重要。

CAE模型通常用沿铸造方向的拉格朗日追踪（分层法）。切下铸坯的二维横截面，按拉拔速度让它通过各冷却区。凝固潜热 $L = 270$ kJ/kg 用焓法处理。

案例2: 航天器的热设计

太空中没有对流，排热完全依靠斯蒂芬-玻尔兹曼辐射。设计基础方程为：

太空背景辐射温度 $T_{\text{space}} \approx 3\,\text{K}$（宇宙微波背景）可以忽略，实际上只用 $\varepsilon\sigma A T_{\text{rad}}^4$。ISS散热器用 $\varepsilon \approx 0.9$、$T_{\text{rad}} \approx 280\,\text{K}$ 运行。$280^4 = 6.15 \times 10^9$，乘以 $\sigma$ 得约 $350\,\text{W/m}^2$ ——这就是每平方米散热器的散热功率。

太阳入射热（$\alpha_s \cdot S_0 \approx 0.3 \times 1361 = 408\,\text{W/m}^2$）的热平衡决定了散热器面积和工作温度。这里太阳吸收率与红外放射率的比值（$\alpha_s / \varepsilon$）是热设计的最关键参数。

案例3: 工业炉的炉壁设计

特别是燃气炉，燃烧气体（CO₂、H₂O）自身参与辐射，单靠斯蒂芬-玻尔兹曼定律不足。需要用加权灰气体(WSGG)模型计算气体放射率，求解辐射输运方程。Ansys Fluent的DO模型、STAR-CCM+的DO/S2S都支持这种功能。

斯蒂芬-玻尔兹曼定律的软件对比

主流CAE工具的辐射功能

开源软件的实现

radiationModel  viewFactor;
viewFactorCoeffs
{
    nBands    1;          // 灰体=1波段
    smoothing true;       // 形态系数平滑化
}
emissivityModel lookup;   // 表查询

ElmerFEM的 Heat Equation 求解器也能用 Radiation = Diffuse Gray 启用灰体面间辐射。但开源软件在大规模形态系数计算（百万面级）的并行效率上有短板，比商用工具慢10〜50倍。

斯蒂芬-玻尔兹曼定律的前沿研究

分光辐射模型

灰体近似（放射率不随波长变化、为常数）在多数工程问题中足够，但以下情况需要分光（光谱）辐射模型：

• 波长选择性材料： 玻璃、陶瓷涂层、半导体硅片等，波长相关放射率变化大
• 太阳光+红外共存： 太阳光峰值~0.5μm，物体红外放射峰值~10μm，两者放射率差异大
• 薄膜干涉： 光学镀膜的热设计

波段模型将波长域分成数个波段，各波段用不同的 $\varepsilon_i$。Ansys Fluent的Non-Gray DO模型、COMSOL的Multi-band radiation都是这种实现。

参与介质的辐射

高温气体（CO₂、H₂O、CO、CH₄）是参与介质（Participating Media），在特定波长范围吸收和放射。此时斯蒂芬-玻尔兹曼定律只作为壁面边界条件，气体内部的辐射传输用辐射输运方程(RTE)描述：

其中 $I$ 是辐射强度，$\kappa$ 是吸收系数，$\sigma_s$ 是散射系数，$\Phi$ 是散射相位函数。壁面处 $I_b = \sigma T^4/\pi$ 就是斯蒂芬-玻尔兹曼定律的贡献。

GPU加速蒙特卡洛法

学术界还在探索用机器学习（PINN: Physics-Informed Neural Networks）近似辐射输运方程，学好的网络能高速解RTE，缩短燃烧模拟的周期。但实务应用还比较有限。

斯蒂芬-玻尔兹曼定律的故障排除

常见错误与解决方案

收敛困难处理

1. 欠松弛因子： 温度更新用 $T^{n+1} = \omega\,T^{n+1}_{\text{calc}} + (1-\omega)\,T^n$ 的方式，取 $\omega = 0.3\sim0.7$。
2. 温度夹紧： 限制每迭代温度变化最大100K，特别是初始几步。
3. 分步启用放射率： 先用 $\varepsilon = 0.1$ 收敛，再逐步上升至实际值（分层）。
4. 改进初值： 即使稳态，给定接近预期温度分布的初值，迭代次数能减半。
5. 网格质量： 辐射面长宽比>10的单元会导致温度振荡。重新划分网格。

Abaqus在 *STEP, NLGEOM 的 *CONTROLS 中可限制 FIELD=TEMPERATURE 的许可增量。Ansys用 NROPT,UNSYM（非对称NR）应对辐射雅可比非对称的问题。