参与性介质辐射传热
理论与物理
什么是参与介质
老师，我第一次听到“参与介质”这个词，普通的空气也会影响辐射吗？
问得好。像氮气 (N₂) 或氧气 (O₂) 这样的对称双原子分子，实际上几乎不吸收也不发射红外线。所以常温下的空气对辐射几乎是“透明”的。
但是，像 CO₂ 或 H₂O 这样的多原子分子就不同了。它们会在与分子振动固有频率相对应的特定波长带（2.7μm, 4.3μm, 15μm 带等）强烈吸收红外线，同时也发射红外线。这种“与辐射能量相互作用的介质”就称为参与介质（participating media）。
原来如此…但是在日常生活中，感觉不到空气在吸收红外线啊？在什么情况下它会变得重要呢？
关键在于温度和尺度。在常温下1米的距离，气体辐射几乎可以忽略。但是，在以下情况下就完全不同了：
- 燃气轮机燃烧室（1500°C以上）：燃烧气体中的CO₂和H₂O产生的辐射占壁面热流密度的30〜50%
- 工业炉/锅炉（1000〜1600°C）：炉内高温气体向被加热物体的主要传热模式变为辐射
- 大规模火灾：火焰的辐射热是向周围蔓延的原因。烟尘颗粒的贡献很大
- 大气辐射：地球的温室效应本身就是CO₂・H₂O的红外吸收现象
粗略地说，在存在高温气体的系统中，忽略辐射会导致致命的设计错误。
啊，壁面热流密度的30〜50%是气体辐射吗！那确实不能忽视啊…只用表面间的角系数不行吗？
没错。表面间的辐射只需要用角系数 $F_{ij}$ 计算表面之间的交换就可以了，但是有参与介质时，介质本身既是发射体同时也是吸收体。光线每次穿过介质时，能量都会增加或减少。所以必须求解辐射传递方程（RTE: Radiative Transfer Equation）这个积分微分方程。
辐射传递方程（RTE）
光是听到RTE就觉得很难…是什么样的方程呢？
RTE是描述“光线沿某一方向前进时，其能量强度如何变化”的方程。在稳态下形式如下：
我们来整理一下各项的含义：
- $I(\mathbf{r}, \hat{\mathbf{s}})$：位置 $\mathbf{r}$、方向 $\hat{\mathbf{s}}$ 的辐射强度 [W/(m²·sr)]
- $\kappa$：吸收系数 [1/m] — 介质吸收辐射的程度
- $\sigma_s$：散射系数 [1/m] — 介质散射辐射的程度
- $I_b = \frac{\sigma T^4}{\pi}$：黑体辐射强度（普朗克函数的方向积分）
- $\Phi(\hat{\mathbf{s}}' \to \hat{\mathbf{s}})$：散射相函数 — 从方向 $\hat{\mathbf{s}}'$ 散射到 $\hat{\mathbf{s}}$ 的概率分布
左边是“光线前进 $ds$ 距离时的强度变化”，右边的四项分别对应“因发射而增加”、“因吸收而减少”、“因散射到其他方向而减少”、“因从其他方向散射过来而增加”。
原来这四项分别代表增减啊。话说这个方程和普通的微分方程有什么不同，难点在哪里呢？
很好的着眼点。RTE之所以难，有三个原因：
- 7维问题：位置3维 × 方向2维 × 波长1维 × 时间1维。考虑全光谱和时间依赖性的话，规模非常庞大
- 积分微分方程：散射项包含对整个立体角 $4\pi$ 的积分，因此所有方向的解是联立的
- 非线性耦合：通过发射项 $I_b \propto T^4$ 与能量方程非线性耦合
所以它的计算成本与Navier-Stokes方程相当——有时甚至更高。
从RTE导出的一个重要关系是辐射源项（加入能量方程的项）：
其中 $G = \int_{4\pi} I \, d\Omega$ 是入射辐射（irradiation），$\sigma$ 是Stefan-Boltzmann常数（$5.67 \times 10^{-8}$ W/(m²·K⁴)）。这一项作为热源被纳入能量方程，从而使辐射与传导、对流耦合。
光学性质与光学厚度
我经常看到“光学厚度”这个术语，具体是什么意思呢？
光学厚度（optical thickness） $\tau$ 是一个无量纲数，表示介质使辐射衰减的程度：
其中 $\beta = \kappa + \sigma_s$ 是消光系数（extinction coefficient），$L$ 是特征长度。
根据光学厚度的值，介质的特性会发生很大变化：
| 光学厚度 $\tau$ | 分类 | 物理行为 | 解法指南 |
|---|---|---|---|
| $\tau < 0.1$ | 光学薄 | 辐射几乎直接穿过介质。吸收/发射的贡献很小 | 可用表面间辐射+薄气体修正来应对 |
| $0.1 < \tau < 3$ | 中间区域 | 吸收、发射、散射都很重要。最难分析 | 需要DOM（S4以上）或蒙特卡洛法 |
| $\tau > 3$ | 光学厚 | 辐射在短距离内被吸收/再发射。表现为扩散行为 | P1近似或Rosseland扩散近似有效 |
例如，燃气轮机燃烧器内富含CO₂的燃烧气体 $\tau \approx 5\sim10$，属于光学厚。另一方面，大气中1米左右的距离下 $\tau \ll 0.1$，属于光学薄。
原来可用的近似方法会随光学厚度而变化啊。中间区域最棘手，这很直观。
气体辐射物性模型
老师，CO₂和H₂O的吸收系数随波长变化很大吧？这个要怎么处理呢？
这正是气体辐射中最棘手的问题。实际的CO₂・H₂O吸收光谱具有极其复杂的结构，由数十万条谱线组成。吸收系数在波长0.01cm⁻¹ 的间隔内就会变化好几个数量级，因此直接计算全光谱的逐线（Line-by-Line, LBL） 法需要数万到数十万个波长带，对于工程计算来说不现实。
因此在实际应用中，会使用以下带模型或全局模型：
WSGGM（Weighted Sum of Gray Gases Model，灰气体加权和模型）
这是应用最广泛的模型。它将真实气体的辐射特性近似为 $N$ 个假想灰气体（gray gas）的加权和：
其中 $a_{\varepsilon,i}(T)$ 是温度依赖的权重系数，$\kappa_i$ 是第 $i$ 个灰气体的吸收系数，$pL$ 是压力×光程长。$i=0$ 对应“透明窗口”，设 $\kappa_0 = 0$。
WSGGM的优点在于，通常只需3〜5个灰气体就能获得工程上足够的精度（与LBL相比误差在5%以内）。代表性的权重系数以Smith et al. (1982) 的最为著名，已默认集成到大多数CFD求解器中。
但也有一些注意事项：
- 前提是均匀温度、均匀组成的气体柱 → 应用于非均匀场时需要修正
- 对于CO₂/H₂O摩尔比变化的系统（煤燃烧 vs. 气体燃烧），应使用不同的系数集
- 高压下（10 atm以上）需要考虑压力展宽效应，需要相应的系数
除了WSGGM还有其他方法吗？
当然有。根据精度和计算成本的权衡来选择使用：
| 模型 | 精度 | 计算成本 | 适用场景 |
|---|---|---|---|
| 逐线法（Line-by-Line, LBL） | 最高（基准解） | 极高 | 基准/验证用 |
| SNB（Statistical Narrow Band，统计窄带模型） | 高 | 高 | 研究目的的详细分析 |
| WSGGM | 工程上足够 | 低 | 工业CFD的标准 |
| 灰气体（Gray Gas） | 低 | 最低 | 概念设计阶段 |
烟尘颗粒的辐射
燃烧分析中经常提到“烟尘（soot）”吧。有烟尘的话辐射会怎么变化呢？
烟尘从辐射的角度来看是非常特殊的存在。气体分子只在特定波长带吸收/发射，而烟尘颗粒则在连续光谱的全范围与辐射相互作用。也就是说，烟尘表现为“近乎灰色”的参与介质。
烟尘引起的吸收系数可用下式近似：
其中 $f_v$ 是烟尘的体积分数（典型值为 $10^{-7} \sim 10^{-5}$），$T$ 是温度 [K]，$C$ 是常数（约 1860 m⁻¹K⁻¹）。总的吸收系数为 $\kappa_{total} = \kappa_{gas} + \kappa_{soot}$。
例如，对于体积分数 $f_v = 10^{-6}$、温度1500K的含烟尘火焰，$\kappa_{soot} \approx 2.8$ m⁻¹。这与CO₂的吸收系数相当，表明即使少量烟尘对辐射的影响也很大。在柴油发动机燃烧室或池火中，烟尘的辐射常常占主导地位。
Hottel的气体辐射图 — 辐射工程的起源
参与介质辐射的研究始于1920年代MIT的Horace Hottel。他通过实验将CO₂・H₂O的总发射率（total emissivity）整理为温度和压力×光程长（pL）的函数图表。在没有数字计算机的时代，工程师仅凭这张图和手算就能设计工业炉。1982年Smith, Shen, Friedman提出的WSGGM，正是对这张Hottel图进行公式拟合，使其成为可在CFD求解器中实现的划时代工作。即使在今天，WSGGM的精度也常常通过与Hottel图的一致性来验证。
RTE各项的物理含义
- 发射项 $\kappa I_b$：根据基尔霍夫定律，好的吸收体也是好的发射体。介质根据局部温度 $T$ 各向同性地发射黑体辐射 $I_b = \sigma T^4/\pi$。气体温度越高，发射的辐射能量越强。
- 吸收项 $-\kappa I$：光线穿过介质时能量被吸收并转化为热。这是比尔定律 $I = I_0 e^{-\kappa s}$ 的微分形式。吸收越强，光传播得越近。
- 外散射项 $-\sigma_s I$：关注方向的光线因散射到其他方向而损失。在存在烟尘、粉尘、液滴等颗粒的系统中很重要。纯气体中散射几乎为零。
- 内散射项 $\frac{\sigma_s}{4\pi}\int I \Phi \, d\Omega'$：从其他方向散射过来的辐射能量加入到关注方向。散射相函数 $\Phi$ 是前向散射占优（米氏散射）还是各向同性散射，会影响分布。
量纲分析与单位制
| 变量 | SI单位 | 典型值/注意事项 |
|---|---|---|
| 辐射强度 $I$ | W/(m²·sr) | 方向依赖量。对整个立体角积分得到辐射通量密度 |
| 吸收系数 $\kappa$ | 1/m | CO₂(1atm, 1000K): 0.1〜10 m⁻¹（波长依赖） |
| 散射系数 $\sigma_s$ | 1/m | 纯气体: ≈0，含烟尘: 0.01〜1 m⁻¹ |
| 光学厚度 $\tau$ | 无量纲 | $\tau = \beta L$。燃烧炉: 5〜10，大气1m: ≪0.01 |
| 黑体辐射强度 $I_b$ | W/(m²·sr) | $I_b = \sigma T^4/\pi$。1500K: 约94 kW/(m²·sr) |
| 入射辐射 $G$ | W/m² | $G = \int_{4\pi} I \, d\Omega$ |
数值解法与实现
离散坐标法（DOM / Discrete Ordinates Method）
实际在计算机上求解RTE要怎么做呢？方向的积分看起来很困难…
最广泛使用的方法是离散坐标法（DOM: Discrete Ordinates Method）。思路很简单，就是用有限个离散方向 $\hat{\mathbf{s}}_m$（$m=1, \ldots, M$）来代表连续的方向 $\hat{\mathbf{s}}$，然后对每个方向分别求解RTE。
在DOM中，RTE 对每个离散方向 $\hat{\mathbf{s}}_m$ 近似为如下形式：
其中 $I_m = I(\mathbf{r}, \hat{\mathbf{s}}_m)$，$w_{m'}$ 是方向 $m'$ 的求积权重。方向的选择广泛使用 S$_N$ 近似，$N$ 的值决定了精度：
| $S_N$ 阶数 | 方向数（3D） | 精度 | 计算成本 | 应用指南 |
|---|---|---|---|---|
| S2 | 8 | 低 | 低 | 初步计算/趋势把握 |
| S4 | 24 | 中 | 中 | 一般工业计算的最低要求 |
| S6 | 48 | 高 | 高 | 标准的详细分析 |
| S8 | 80 | 非常高 | 非常高 | 光学薄区域的精密分析 |
增加方向数精度会提高，但计算成本也会上升啊。实际工作中S4是标准吗？
是的，实际工作中S4〜S6用得最多。S2的方向分辨率太低，容易出现称为射线效应（ray effect）的人工条纹。相反，S8以上计算成本会急剧上升，所以只在精度要求非常高的情况下使用。
DOM的优势在于：
- 不受光学厚度限制（从光学薄到光学厚都能应对）
- 也能应对散射强的系统
- 与非结构网格兼容性好
- 几乎所有的CFD求解器都已实现
弱点是会形成“方向数 × 空间网格数”的联立方程，因此内存消耗大。
P1近似（球谐函数法）
除了DOM还有其他方法吗？我听说过P1近似。
P1近似是另一种方法，它假设辐射强度 $I(\mathbf{r}, \hat{\mathbf{s}})$ 的方向依赖性可以用球谐函数展开，并且只保留一阶项。这相当于假设辐射强度在方向上是近似各向同性的，或者其方向变化可以用一个简单的线性函数来描述。
在P1近似下，RTE可以简化为一个关于入射辐射 $G$ 的椭圆型偏微分方程：
这个方程看起来很像扩散方程，因此P1近似有时也被称为辐射扩散近似。它的优点是计算成本极低，因为只需要求解一个标量方程（对于 $G$），而不是像DOM那样求解多个方向上的方程。
然而，P1近似有严格的适用条件：
- 介质必须是光学厚的（$\tau \gg 1$），这样辐射在局部区域才能达到近似平衡，方向分布才接近各向同性。
- 散射必须是各向同性的，或者至少是近似各向同性的。
- 对于光学薄的区域或存在强方向性光源（如激光、小开口）的情况，P1近似误差很大。
因此，P1近似常用于核反应堆、恒星内部等光学极厚的系统，或者作为复杂问题中光学厚区域的快速估算方法。在一般的工业燃烧计算中，如果光学厚度足够大，也可以考虑使用P1近似来节省计算资源。
明白了，就是用一个扩散方程来近似辐射传递。那有没有介于DOM和P1之间的方法呢？
有的，比如简化球谐函数法（SPN）。它保留了球谐函数展开的更高阶项（如P3, P5, P7），从而能在更宽的光学厚度范围内保持精度，同时计算成本仍低于全方向的DOM。SPN方法在医学光学、半导体工艺等领域的辐射传输问题中应用较多，在燃烧CFD中也逐渐受到关注。
另外，对于光学厚度变化剧烈的复杂系统，还有一种策略是混合方法：在光学厚的区域使用P1近似，在光学薄的区域使用DOM或射线追踪法，然后在交界区域进行耦合。这需要复杂的算法实现，但能有效平衡精度和效率。
蒙特卡洛法（Monte Carlo Method）
我还听说过蒙特卡洛法，它和DOM、P1有什么不同？
蒙特卡洛辐射法（MCRT）是一种基于概率统计的方法。它的核心思想不是直接求解RTE，而是追踪大量“光子束”或“能量束”的随机行走过程，通过统计这些束的归宿（如被吸收、到达某个表面）来估算辐射热流、温度分布等。
基本步骤是：
- 发射：从辐射源（高温表面或高温气体体积元）按照其发射特性随机发射能量束。
- 传输：束在介质中沿直线传播，传播距离根据介质的消光系数按概率分布随机确定。
- 相互作用：在传播终点，根据吸收系数和散射系数的比例，随机决定束是被吸收还是被散射。
- 如果被吸收，能量沉积到当地网格，该束历史结束。
- 如果被散射，则根据散射相函数随机确定新的传播方向，然后返回步骤2继续传输。
- 如果束到达边界表面，则根据表面的反射特性（漫反射、镜面反射、吸收）决定其下一步是反射还是被吸收。
重复以上过程追踪大量束的历史，最后对所有束的能量沉积和边界热流进行统计平均，得到结果。
听起来很直观，就像模拟实际光子的行为。它的优缺点是什么？
是的，它的物理图像非常清晰。主要优点包括：
- 精度高，无方向离散误差：方向是随机采样的，不存在DOM的射线效应。
- 几何适应性极强：可以处理极其复杂的几何形状，包括非结构网格、移动边界等，因为束的传输是沿着直线进行的，不依赖于网格拓扑。
- 易于处理复杂物性：可以相对容易地纳入全光谱模型、各向异性散射、复杂表面反射特性等。
- 天然并行：每个束的历史是独立的，非常适合大规模并行计算。
主要缺点：
- 计算成本高，收敛慢：为了获得统计上平滑的结果，需要追踪海量的束（数百万甚至数十亿）。计算成本与精度（统计噪声水平）的平方成反比。要减少一半的噪声，需要增加四倍的束数量。
- 统计噪声：结果中存在固有的统计波动，对于梯度很大的物理量（如壁面热流峰值）可能需要非常多的样本才能准确捕捉。
- 与流场耦合困难：MCRT通常是“后处理”式的，即给定一个温度场和物性场，计算辐射场。要实现与流动、传热的强耦合迭代，需要在每次迭代中都进行昂贵的MCRT计算，成本很高。
因此，MCRT通常用于：1) 为其他方法（如DOM）提供基准解；2) 几何或物性极其复杂，其他方法难以应用的场合；3) 作为一次性的辐射热流计算工具，用于设计验证。
在CFD中的实现与耦合
在实际的CFD计算中，辐射模型是怎么和流动、传热方程耦合在一起的呢？
这是一个关键问题。辐射与流动、传热的耦合是双向的：
- 辐射 → 能量方程：辐射计算得出的辐射源项 $\nabla \cdot \mathbf{q}_r$ 作为源项加入到流体的能量方程中（或固体导热方程中）。这个源项代表由于辐射的吸收和发射导致的当地净能量增益或损失。
- 温度/组分 → 辐射：流场计算出的温度 $T$ 和气体组分浓度（决定 $\kappa, \sigma_s$）是辐射计算所需的输入。
在瞬态计算中，这种耦合需要在每个时间步进行。在稳态计算中，则通过迭代来实现耦合。
典型的耦合求解策略有：
- 松散耦合（Loosely Coupled / Segregated）：在一个迭代步内，先假设辐射源项不变求解流动和能量方程，更新温度场；然后用新的温度场求解辐射方程，更新辐射源项；如此交替进行直到收敛。
この記事の評価ご回答ありがとうございます！参考に
なったもっと
詳しく誤りを
報告