主题012：辐射换热的有限差分法

摘要

有限差分法是求解辐射换热问题最经典、最基础的数值方法之一。本主题将系统介绍辐射换热有限差分法的理论基础、离散格式、边界条件处理以及实际应用。通过详细的数学推导和丰富的Python代码示例，帮助读者掌握从一维到多维、从稳态到瞬态、从灰体到非灰体的各种辐射换热问题的有限差分求解技术。

关键词： 有限差分、空间离散、角度离散、截断误差、数值稳定性、辐射传递方程

---

1. 引言

1.1 有限差分法的发展历史

有限差分法（Finite Difference Method, FDM）是数值分析中最古老的方法之一，其历史可以追溯到18世纪欧拉和拉格朗日的时代。在辐射换热领域，有限差分法的应用始于20世纪50年代，随着计算机技术的发展而逐渐成熟。

早期的辐射换热研究主要依赖解析解和近似方法，如Hottel的 zone method 和 Gebhart 的吸收因子法。然而，这些方法在处理复杂几何形状和非均匀介质时面临巨大挑战。有限差分法的引入为辐射换热问题的数值求解提供了统一而灵活的框架。

20世纪70-80年代，随着计算流体力学（CFD）的兴起，有限差分法在辐射换热中的应用得到了快速发展。研究人员开发了各种离散格式，如迎风格式、中心差分格式、QUICK格式等，并研究了它们的数值特性。

进入21世纪，有限差分法与其他数值方法（如有限体积法、有限元法）的结合，以及高性能计算技术的应用，使得复杂辐射换热问题的求解成为可能。

1.2 有限差分法的优势与局限

优势：

1. 概念简单直观：有限差分法基于泰勒展开，用差商代替微商，物理意义清晰，易于理解和实现。

2. 编程实现容易：相比有限元法需要复杂的网格生成和刚度矩阵组装，有限差分法的程序结构相对简单。

3. 计算效率高：在规则网格上，有限差分法的计算效率通常高于有限元法，特别适合大规模并行计算。

4. 理论分析成熟：有限差分法的稳定性、收敛性和精度分析有完善的数学理论支持。

局限：

1. 几何适应性差：有限差分法对复杂几何形状的处理能力较弱，难以处理曲线边界和不规则网格。

2. 守恒性不足：传统的有限差分法难以严格保证物理量的守恒，这在辐射换热问题中尤为重要。

3. 边界条件处理复杂：在曲线边界上，边界条件的离散需要特殊处理，如虚拟节点法或坐标变换法。

4. 网格限制：有限差分法通常要求结构化网格，对自适应网格和非结构化网格的支持不如有限元法。

1.3 本主题的组织结构

本主题按照由浅入深的原则组织内容：

• 第2节介绍有限差分法的基本原理，包括泰勒展开、差商近似和截断误差分析。
• 第3节详细讲解空间离散方法，从一维到多维的各种差分格式。
• 第4节讨论角度离散方法，这是辐射换热特有的离散维度。
• 第5节分析波谱离散方法，处理非灰体辐射问题。
• 第6节研究数值稳定性和收敛性。
• 第7-10节通过四个完整的Python案例展示有限差分法的应用。

---

2. 有限差分法的基本原理

2.1 泰勒展开与差商近似

有限差分法的核心思想是用离散的差商代替连续的微商。这基于函数的泰勒展开：

对于光滑函数 $f(x)$，在点 $x_i$ 附近有：

$$f(x_{i+1})=f(x_i)+f^{\prime }(x_i)\Delta x+\frac{f^{\prime \prime }(x_i)}{2!}(\Delta x{)}^2+\frac{f^{\prime \prime \prime }(x_i)}{3!}(\Delta x{)}^3+O((\Delta x{)}^4)$$

$$f(x_{i-1})=f(x_i)-f^{\prime }(x_i)\Delta x+\frac{f^{\prime \prime }(x_i)}{2!}(\Delta x{)}^2-\frac{f^{\prime \prime \prime }(x_i)}{3!}(\Delta x{)}^3+O((\Delta x{)}^4)$$

其中 $\Delta x=x_{i+1}-x_i$ 是网格间距。

一阶导数的差分近似：

从前向泰勒展开，可以得到前向差分：

$$f^{\prime }(x_i)\approx \frac{f(x_{i+1})-f(x_i)}{\Delta x}+O(\Delta x)$$

从后向泰勒展开，可以得到后向差分：

$$f^{\prime }(x_i)\approx \frac{f(x_i)-f(x_{i-1})}{\Delta x}+O(\Delta x)$$

将两个展开式相减，可以得到中心差分：

$$f^{\prime }(x_i)\approx \frac{f(x_{i+1})-f(x_{i-1})}{2\Delta x}+O((\Delta x{)}^2)$$

二阶导数的差分近似：

将两个泰勒展开式相加，消去一阶导数项：

$$f^{\prime \prime }(x_i)\approx \frac{f(x_{i+1})-2f(x_i)+f(x_{i-1})}{(\Delta x{)}^2}+O((\Delta x{)}^2)$$

这就是二阶中心差分，具有二阶精度。

2.2 截断误差与精度阶数

截断误差是差分近似与精确导数之间的差异。例如，前向差分的截断误差为：

$$R=-\frac{\Delta x}{2}{f}^{\prime \prime }(x_i)-\frac{(\Delta x{)}^2}{6}{f}^{\prime \prime \prime }(x_i)+O((\Delta x{)}^3)$$

当 $\Delta x\to 0$ 时，如果 $R=O((\Delta x{)}^p)$，则称该差分格式具有p阶精度。

• 前向差分和后向差分：一阶精度，$O(\Delta x)$
• 中心差分：二阶精度，$O((\Delta x{)}^2)$

高阶差分格式：

通过使用更多的网格点，可以构造更高阶的差分格式。例如，四阶精度的中心差分：

$$f^{\prime }(x_i)\approx \frac{-f(x_{i+2})+8f(x_{i+1})-8f(x_{i-1})+f(x_{i-2})}{12\Delta x}+O((\Delta x{)}^4)$$

2.3 差分格式的分类

根据所使用的网格点位置和数量，差分格式可以分为：

1. 按方向分类：

• 前向差分（Forward Difference）：使用 $x_i$ 和 $x_{i+1}$
• 后向差分（Backward Difference）：使用 $x_{i-1}$ 和 $x_i$
• 中心差分（Central Difference）：使用 $x_{i-1}$ 和 $x_{i+1}$

2. 按精度分类：

• 一阶精度：$O(\Delta x)$
• 二阶精度：$O((\Delta x{)}^2)$
• 高阶精度：$O((\Delta x{)}^p)$，$p\ge 3$

3. 按性质分类：

• 显式格式：当前时刻的值直接由已知值计算
• 隐式格式：需要求解方程组
• 半隐式格式：结合显式和隐式的特点

---

3. 空间离散方法

3.1 一维空间离散

考虑一维稳态辐射传递方程：

$$\mu \frac{dI}{dx}=-\beta I+\kappa I_b+\frac{{\sigma }_s}{4\pi }{\int }_{4\pi }I({\Omega }^{\prime })\Phi (\Omega ,{\Omega }^{\prime })d{\Omega }^{\prime }$$

其中 $\mu =\cos \theta$ 是方向余弦，$\beta =\kappa +{\sigma }_s$ 是消光系数。

迎风格式（Upwind Scheme）：

根据传输方向选择差分格式：

• 当 $\mu >0$（正向传输）：使用后向差分
  $$\mu \frac{I_i-I_{i-1}}{\Delta x}=-{\beta }_i{I}_i+S_i$$

• 当 $\mu <0$（反向传输）：使用前向差分
  $$\mu \frac{I_{i+1}-I_i}{\Delta x}=-{\beta }_i{I}_i+S_i$$

其中 $S_i={\kappa }_i{I}_{b,i}+\frac{{\sigma }_{s,i}}{4\pi }\int I\Phi d{\Omega }^{\prime }$ 是源项。

迎风格式具有一阶精度，但具有较好的稳定性，特别适合对流主导的问题。

** diamond差分格式：**

为了提高精度，可以使用 diamond 差分格式：

$$I_{i+1/2}=\frac{1}{2}(I_i+I_{i+1})$$

将辐射传递方程在控制体上积分：

$$\mu (I_{i+1/2}-I_{i-1/2})=-{\beta }_i\Delta x\cdot I_i+S_i\Delta x$$

结合 diamond 关系，可以得到：

$$I_i=\frac{2\mu I_{i-1/2}+S_i\Delta x}{2\mu +{\beta }_i\Delta x},\mu >0$$

指数差分格式：

考虑局部解析解，可以得到指数差分格式：

$$I_{i+1/2}=I_{i-1/2}{e}^{-{\beta }_i\Delta x/\mu }+\frac{S_i}{{\beta }_i}(1-e^{-{\beta }_i\Delta x/\mu })$$

指数格式在光学厚极限下自动退化为扩散近似，具有很好的渐近特性。

3.2 二维空间离散

在二维情况下，辐射传递方程为：

$$\xi \frac{\partial I}{\partial x}+\eta \frac{\partial I}{\partial y}=-\beta I+S$$

其中 $\xi =\sin \theta \cos \varphi$，$\eta =\sin \theta \sin \varphi$。

方向分裂法：

将二维问题分解为两个一维问题：

$$\xi \frac{\partial I}{\partial x}=-\frac{\beta }{2}I+\frac{S}{2}$$

$$\eta \frac{\partial I}{\partial y}=-\frac{\beta }{2}I+\frac{S}{2}$$

分别用一维方法求解，然后进行迭代。

九点差分格式：

使用二维泰勒展开，可以构造九点差分格式：

$$\xi \frac{I_{i+1,j}-I_{i-1,j}}{2\Delta x}+\eta \frac{I_{i,j+1}-I_{i,j-1}}{2\Delta y}=-{\beta }_{i,j}{I}_{i,j}+S_{i,j}$$

这种格式具有二阶精度，但可能产生数值振荡。

3.3 三维空间离散

三维辐射传递方程：

$$\xi \frac{\partial I}{\partial x}+\eta \frac{\partial I}{\partial y}+\mu \frac{\partial I}{\partial z}=-\beta I+S$$

27点差分格式：

扩展到三维，可以使用27点差分格式，考虑所有方向的贡献。

有限体积离散：

在控制体 $V_{i,j,k}$ 上积分：

$${\int }_{V_{i,j,k}}\left(\xi \frac{\partial I}{\partial x}+\eta \frac{\partial I}{\partial y}+\mu \frac{\partial I}{\partial z}\right)dV={\int }_{V_{i,j,k}}(-\beta I+S)dV$$

应用散度定理：

$${\oint }_{\partial V}I(\mathbf{n}\cdot \mathbf{\Omega })dA={\int }_V(-\beta I+S)dV$$

其中 $\mathbf{n}$ 是控制体表面的外法向向量。

---

4. 角度离散方法

4.1 离散坐标法（DOM）

离散坐标法（Discrete Ordinates Method, DOM）是辐射换热中最常用的角度离散方法。基本思想是将连续的立体角空间离散为有限个方向（射线）。

S-N 方法：

S-N方法将单位球面上的方向离散为 $N(N+2)$ 个方向（对于三维问题）。常用的离散阶数包括：

• S2：8个方向
• S4：24个方向
• S8：80个方向
• S16：288个方向

方向余弦的选择：

对于S4方法，方向余弦和权重如下表所示：

方向	$\mu$	$\eta$	$\xi$	权重 $w$
1	0.2958759	0.2958759	0.9082483	0.5238095
2	0.2958759	0.9082483	0.2958759	0.5238095
…	…	…	…	…

对称性要求：

方向集合必须满足以下对称性条件：

1. 对于每个方向 $(\mu ,\eta ,\xi )$，必须有 $(-\mu ,\eta ,\xi )$、$(\mu ,-\eta ,\xi )$ 等对称方向
2. 所有方向的权重之和等于 $4\pi$
3. 方向余弦的加权平均为零

4.2 角度离散的控制方程

离散后的辐射传递方程为：

$${\mathbf{\Omega }}_m\cdot \nabla I_m=-\beta I_m+\kappa I_b+\frac{{\sigma }_s}{4\pi }\sum _{m^{\prime }=1}^M{w}_{m^{\prime }}{I}_{m^{\prime }}{\Phi }_{m,m^{\prime }}$$

其中 $m=1,2,...,M$ 是方向索引，$M$ 是总方向数，$w_{m^{\prime }}$ 是方向权重。

散射项的处理：

对于各向同性散射，${\Phi }_{m,m^{\prime }}=1$，散射项简化为：

$$\frac{{\sigma }_s}{4\pi }\sum _{m^{\prime }=1}^M{w}_{m^{\prime }}{I}_{m^{\prime }}={\sigma }_{s}G$$

其中 $G=\frac{1}{4\pi }{\int }_{4\pi }Id\Omega$ 是入射辐射。

对于各向异性散射，需要计算相函数：

$${\Phi }_{m,m^{\prime }}=1+\sum _{l=1}^L{a}_l{P}_l(\cos {\theta }_{m,m^{\prime }})$$

其中 $P_l$ 是勒让德多项式，$a_l$ 是展开系数。

4.3 射线效应与虚假散射

**射线效应（Ray Effect）：

射线效应是离散坐标法的固有缺陷。当源项或边界条件具有强烈的方向性时，离散的方向可能无法准确捕捉辐射的传播路径，导致数值解出现非物理的振荡。

缓解射线效应的方法：

1. 增加离散方向数（提高S-N阶数）
2. 使用有限体积法代替离散坐标法
3. 采用自适应角度加密
4. 使用修正的离散坐标法

**虚假散射（False Scattering）：

虚假散射是由于空间离散引起的数值扩散效应。即使物理上介质是透明的，数值解也会表现出类似散射的行为。

减少虚假散射的方法：

1. 使用高阶空间离散格式
2. 减小网格尺寸
3. 使用指数格式或解析格式
4. 采用自适应网格加密

---

5. 波谱离散方法

5.1 灰体假设与全波谱模型

灰体假设：

在灰体假设下，辐射物性（吸收系数、散射系数）与波长无关。辐射传递方程简化为：

$$\mathbf{\Omega }\cdot \nabla I=-\beta I+\kappa I_b+\frac{{\sigma }_s}{4\pi }{\int }_{4\pi }I\Phi d{\Omega }^{\prime }$$

其中 $I_b=\frac{\sigma T^4}{\pi }$ 是黑体辐射强度。

灰体假设大大简化了计算，但在处理选择性吸收介质（如燃烧产物中的CO₂和H₂O）时精度不足。

全波谱模型：

考虑波谱依赖性，辐射传递方程变为：

$$\mathbf{\Omega }\cdot \nabla I_{\lambda }=-{\beta }_{\lambda }{I}_{\lambda }+{\kappa }_{\lambda }{I}_{b,\lambda }+\frac{{\sigma }_{s,\lambda }}{4\pi }{\int }_{4\pi }{I}_{\lambda }{\Phi }_{\lambda }d{\Omega }^{\prime }$$

需要在大量波谱点（通常数千到数万）上求解方程，计算量巨大。

5.2 波谱带模型

**窄谱带模型（Narrow Band Model, NBM）：

窄谱带模型假设在一个小波谱范围内，辐射物性可以用统计方法描述。常用的模型包括：

1. Elsasser模型：假设谱线等间距排列
2. 统计模型：假设谱线随机分布
3. Goody模型：基于统计理论的简化模型

窄谱带模型的吸收率计算：

$${\bar{\kappa }}_{\Delta }=\frac{1}{\Delta \lambda }{\int }_{\Delta \lambda }{\kappa }_{\lambda }d\lambda$$

**宽谱带模型（Wide Band Model, WBM）：

宽谱带模型将整个吸收带作为一个整体处理，使用指数宽谱带模型：

$$\bar{\kappa }=C\cdot \omega \cdot {\left(\frac{\omega }{\delta }\right)}^{n-1}$$

其中 $C$ 是带强度参数，$\omega$ 是带宽参数，$\delta$ 是线间距，$n$ 是带宽指数。

5.3 加权求和灰气体模型（WSGGM）

WSGGM是工程中最常用的非灰体辐射模型。基本思想是将非灰体气体等效为若干灰气体的加权求和：

$$\epsilon =\sum _{k=1}^K{a}_k(T)(1-e^{-{\kappa }_{k}L})$$

其中 $K$ 是灰气体数（通常4-5个），$a_k(T)$ 是温度相关的权重系数，${\kappa }_k$ 是灰气体吸收系数，$L$ 是平均射线程长。

WSGGM的优点：

1. 计算效率高，只需在每个灰气体上求解一次辐射传递方程
2. 精度足够满足工程需求
3. 参数可以通过实验数据或高分辨率计算拟合得到

WSGGM的局限性：

1. 参数依赖于具体的压力和温度范围
2. 对于强非灰体效应（如低温下的CO₂）精度下降
3. 难以处理非均匀介质

---

6. 数值稳定性与收敛性分析

6.1 稳定性分析

冯·诺依曼稳定性分析：

对于线性差分方程，可以使用冯·诺依曼方法分析稳定性。假设解的形式为：

$$I_i^n=G^n{e}^{ik{x}_i}$$

其中 $G$ 是放大因子，$k$ 是波数。

稳定性要求 $|G|\le 1$ 对所有波数成立。

CFL条件：

对于显式格式，稳定性通常要求满足CFL（Courant-Friedrichs-Lewy）条件：

$$\frac{|\mu |\Delta t}{\Delta x}\le 1$$

在辐射传递问题中，由于光速极大，通常采用隐式格式以避免严格的时间步长限制。

6.2 收敛性分析

Lax等价定理：

对于适定的线性初值问题，如果差分格式是相容的，那么稳定性等价于收敛性。

相容性：

当 $\Delta x\to 0$ 和 $\Delta t\to 0$ 时，差分方程趋近于微分方程。

收敛阶数：

如果数值解与精确解的误差满足：

$$\Vert I_{numerical}-I_{exact}\Vert =O((\Delta x{)}^p+(\Delta t{)}^q)$$

则称格式具有空间 $p$ 阶和时间 $q$ 阶收敛精度。

6.3 误差估计与网格加密

后验误差估计：

通过比较不同网格上的数值解，可以估计离散误差：

$$\epsilon \approx \frac{\Vert I_h-I_{2h}\Vert }{2^p-1}$$

其中 $h$ 和 $2h$ 是粗细网格尺寸，$p$ 是格式精度阶数。

自适应网格加密：

基于误差估计，可以在误差大的区域自动加密网格：

1. 计算每个单元的误差指示器
2. 标记误差超过阈值的单元
3. 对标记单元进行加密（细分）
4. 重复直到满足精度要求或达到最大网格数

---

7. 案例一：一维灰体介质辐射传递

7.1 问题描述

考虑一维平行平板间的灰体介质辐射传递。两平板温度分别为 $T_1=1000$ K 和 $T_2=500$ K，板间距 $L=1$ m。介质为灰体，吸收系数 $\kappa =1.0$ m⁻¹，散射系数 ${\sigma }_s=0.5$ m⁻¹。

7.2 数学模型

稳态一维辐射传递方程：

$$\mu \frac{dI}{dx}=-\beta I+\kappa I_b+\frac{{\sigma }_s}{2}{\int }_{-1}^{1}I({\mu }^{\prime })d{\mu }^{\prime }$$

边界条件：

• $x=0$：$I(0,\mu >0)={\epsilon }_1{I}_{b1}+(1-{\epsilon }_1)I(0,\mu <0)$
• $x=L$：$I(L,\mu <0)={\epsilon }_2{I}_{b2}+(1-{\epsilon }_2)I(L,\mu >0)$

7.3 有限差分离散

使用S4离散坐标法和迎风格式：

对于 ${\mu }_m>0$：
$${\mu }_m\frac{I_{i,m}-I_{i-1,m}}{\Delta x}=-\beta I_{i,m}+\kappa I_{b,i}+\frac{{\sigma }_s}{2}\sum _{m^{\prime }}{w}_{m^{\prime }}{I}_{i,m^{\prime }}$$

对于 ${\mu }_m<0$：
$${\mu }_m\frac{I_{i+1,m}-I_{i,m}}{\Delta x}=-\beta I_{i,m}+\kappa I_{b,i}+\frac{{\sigma }_s}{2}\sum _{m^{\prime }}{w}_{m^{\prime }}{I}_{i,m^{\prime }}$$

7.4 Python实现

"""
案例一：一维灰体介质辐射传递
使用S4离散坐标法和迎风格式
"""

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation, PillowWriter
import matplotlib
matplotlib.use('Agg')

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

print("=" * 70)
print("案例一：一维灰体介质辐射传递")
print("=" * 70)

# 物理参数
T1 = 1000  # K, 左壁温度
T2 = 500   # K, 右壁温度
L = 1.0    # m, 板间距
kappa = 1.0  # m^-1, 吸收系数
sigma_s = 0.5  # m^-1, 散射系数
beta = kappa + sigma_s  # 消光系数
omega = sigma_s / beta  # 散射反照率
sigma = 5.670374e-8  # Stefan-Boltzmann常数

# 黑体辐射强度
Ib1 = sigma * T1**4 / np.pi
Ib2 = sigma * T2**4 / np.pi

# S4离散坐标 (方向余弦和权重)
mu_pos = np.array([0.2958759, 0.9082483])  # 正方向
w_pos = np.array([0.5238095, 0.5238095])   # 对应权重
mu = np.concatenate([-mu_pos[::-1], mu_pos])  # 全方向 [-0.908, -0.296, 0.296, 0.908]
w = np.concatenate([w_pos[::-1], w_pos])     # 对应权重

print(f"\n物理参数:")
print(f"  左壁温度: {T1} K")
print(f"  右壁温度: {T2} K")
print(f"  板间距: {L} m")
print(f"  吸收系数: {kappa} m^-1")
print(f"  散射系数: {sigma_s} m^-1")
print(f"  光学厚度: {beta * L:.2f}")
print(f"  散射反照率: {omega:.2f}")

# 空间离散
Nx = 50
dx = L / (Nx - 1)
x = np.linspace(0, L, Nx)

print(f"\n数值参数:")
print(f"  网格数: {Nx}")
print(f"  网格间距: {dx:.4f} m")
print(f"  方向数: {len(mu)}")

# 初始化辐射强度
I = np.zeros((Nx, len(mu)))

# 源项迭代求解
max_iter = 1000
tolerance = 1e-6

print("\n开始迭代求解...")

for iteration in range(max_iter):
    I_old = I.copy()
    
    # 计算入射辐射 G = sum(w_m * I_m)
    G = np.zeros(Nx)
    for i in range(Nx):
        G[i] = np.sum(w * I[i, :])
    
    # 计算源项 S = kappa * Ib + sigma_s / (4*pi) * G
    # 假设介质温度线性分布
    T_medium = T1 + (T2 - T1) * x / L
    Ib_medium = sigma * T_medium**4 / np.pi
    S = kappa * Ib_medium + sigma_s * G / (4 * np.pi)
    
    # 正方向扫描 (mu > 0)
    for m in range(len(mu)):
        if mu[m] > 0:
            # 左边界条件
            I[0, m] = Ib1  # 黑体边界
            
            # 空间扫描
            for i in range(1, Nx):
                # 迎风格式
                I[i, m] = (mu[m] * I[i-1, m] + dx * S[i]) / (mu[m] + beta * dx)
    
    # 负方向扫描 (mu < 0)
    for m in range(len(mu)):
        if mu[m] < 0:
            # 右边界条件
            I[-1, m] = Ib2  # 黑体边界
            
            # 空间扫描
            for i in range(Nx-2, -1, -1):
                # 迎风格式
                I[i, m] = (abs(mu[m]) * I[i+1, m] + dx * S[i]) / (abs(mu[m]) + beta * dx)
    
    # 检查收敛
    error = np.max(np.abs(I - I_old))
    if iteration % 100 == 0:
        print(f"  迭代 {iteration}: 最大误差 = {error:.2e}")
    
    if error < tolerance:
        print(f"\n收敛！总迭代次数: {iteration + 1}")
        break

# 计算辐射热流
q = np.zeros(Nx)
for i in range(Nx):
    q[i] = 2 * np.pi * np.sum(w * mu * I[i, :])

# 计算入射辐射
G_final = np.zeros(Nx)
for i in range(Nx):
    G_final[i] = np.sum(w * I[i, :])

print(f"\n计算结果:")
print(f"  左壁辐射热流: {q[0]:.2f} W/m²")
print(f"  右壁辐射热流: {q[-1]:.2f} W/m²")
print(f"  热流守恒误差: {abs(q[0] - q[-1]):.4f} W/m²")

# 可视化
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 辐射强度分布
ax1 = axes[0, 0]
for m in range(len(mu)):
    ax1.plot(x, I[:, m], linewidth=2, label=f'μ={mu[m]:.3f}')
ax1.set_xlabel('Position (m)')
ax1.set_ylabel('Radiation Intensity (W/(m²·sr))')
ax1.set_title('Radiation Intensity Distribution')
ax1.legend(fontsize=8)
ax1.grid(True, alpha=0.3)

# 辐射热流分布
ax2 = axes[0, 1]
ax2.plot(x, q, 'b-', linewidth=2)
ax2.axhline(y=q[0], color='r', linestyle='--', alpha=0.5, label=f'q={q[0]:.1f} W/m²')
ax2.set_xlabel('Position (m)')
ax2.set_ylabel('Radiative Heat Flux (W/m²)')
ax2.set_title('Radiative Heat Flux Distribution')
ax2.legend()
ax2.grid(True, alpha=0.3)

# 入射辐射分布
ax3 = axes[1, 0]
ax3.plot(x, G_final, 'g-', linewidth=2)
ax3.set_xlabel('Position (m)')
ax3.set_ylabel('Incident Radiation (W/m²)')
ax3.set_title('Incident Radiation Distribution')
ax3.grid(True, alpha=0.3)

# 角度分布 (在x=L/2处)
ax4 = axes[1, 1]
x_mid = Nx // 2
ax4.bar(range(len(mu)), I[x_mid, :], color='orange', alpha=0.7)
ax4.set_xlabel('Direction Index')
ax4.set_ylabel('Intensity (W/(m²·sr))')
ax4.set_title(f'Angular Distribution at x={x[x_mid]:.2f}m')
ax4.set_xticks(range(len(mu)))
ax4.set_xticklabels([f'{m:.3f}' for m in mu], rotation=45)
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('case1_1d_gray_medium.png', dpi=150)
plt.close()
print("\n已保存: case1_1d_gray_medium.png")

# 生成收敛过程动画
print("\n生成收敛过程动画...")

# 重新计算并保存中间结果
I_history = []
I_temp = np.zeros((Nx, len(mu)))

for iteration in range(min(200, max_iter)):
    # 计算入射辐射
    G_temp = np.zeros(Nx)
    for i in range(Nx):
        G_temp[i] = np.sum(w * I_temp[i, :])
    
    S_temp = kappa * Ib_medium + sigma_s * G_temp / (4 * np.pi)
    
    # 正方向扫描
    for m in range(len(mu)):
        if mu[m] > 0:
            I_temp[0, m] = Ib1
            for i in range(1, Nx):
                I_temp[i, m] = (mu[m] * I_temp[i-1, m] + dx * S_temp[i]) / (mu[m] + beta * dx)
    
    # 负方向扫描
    for m in range(len(mu)):
        if mu[m] < 0:
            I_temp[-1, m] = Ib2
            for i in range(Nx-2, -1, -1):
                I_temp[i, m] = (abs(mu[m]) * I_temp[i+1, m] + dx * S_temp[i]) / (abs(mu[m]) + beta * dx)
    
    if iteration % 5 == 0:
        I_history.append(I_temp.copy())

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# 左图：辐射强度
ax1 = axes[0]
lines1 = []
for m in range(len(mu)):
    line, = ax1.plot([], [], linewidth=2, label=f'μ={mu[m]:.3f}')
    lines1.append(line)
ax1.set_xlim([0, L])
ax1.set_ylim([0, np.max(I) * 1.1])
ax1.set_xlabel('Position (m)')
ax1.set_ylabel('Radiation Intensity (W/(m²·sr))')
ax1.set_title('Convergence of Radiation Intensity')
ax1.legend(fontsize=8)
ax1.grid(True, alpha=0.3)

# 右图：辐射热流
ax2 = axes[1]
line2, = ax2.plot([], [], 'b-', linewidth=2)
ax2.set_xlim([0, L])
ax2.set_ylim([min(q) * 1.2, max(q) * 1.2])
ax2.set_xlabel('Position (m)')
ax2.set_ylabel('Radiative Heat Flux (W/m²)')
ax2.set_title('Convergence of Heat Flux')
ax2.grid(True, alpha=0.3)

text = fig.text(0.5, 0.95, '', ha='center', fontsize=12)

def init():
    for line in lines1:
        line.set_data([], [])
    line2.set_data([], [])
    text.set_text('')
    return lines1 + [line2, text]

def update(frame):
    I_frame = I_history[frame]
    
    # 更新辐射强度曲线
    for m, line in enumerate(lines1):
        line.set_data(x, I_frame[:, m])
    
    # 计算并更新热流
    q_frame = np.zeros(Nx)
    for i in range(Nx):
        q_frame[i] = 2 * np.pi * np.sum(w * mu * I_frame[i, :])
    line2.set_data(x, q_frame)
    
    text.set_text(f'Iteration: {frame * 5}')
    return lines1 + [line2, text]

ani = FuncAnimation(fig, update, frames=len(I_history), 
                    init_func=init, interval=200, blit=False)
writer = PillowWriter(fps=5)
ani.save('case1_convergence_animation.gif', writer=writer)
plt.close()
print("已保存: case1_convergence_animation.gif")

print("\n" + "=" * 70)
print("案例一完成！")
print("=" * 70)

---

8. 案例二：二维方腔辐射换热

8.1 问题描述

考虑一个二维方腔，边长为 $L=1$ m。左壁温度为 $T_h=1000$ K，右壁温度为 $T_c=500$ K，上下壁绝热。腔内充满灰体介质，吸收系数 $\kappa =1.0$ m⁻¹。

8.2 数学模型

二维稳态辐射传递方程：

$$\xi \frac{\partial I}{\partial x}+\eta \frac{\partial I}{\partial y}=-\kappa I+\kappa I_b$$

其中 $(\xi ,\eta )=(\sin \theta \cos \varphi ,\sin \theta \sin \varphi )$。

8.3 Python实现

"""
案例二：二维方腔辐射换热
使用S4离散坐标法和迎风格式
"""

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation, PillowWriter
from matplotlib import cm
import matplotlib
matplotlib.use('Agg')

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

print("=" * 70)
print("案例二：二维方腔辐射换热")
print("=" * 70)

# 物理参数
L = 1.0  # m, 方腔边长
Th = 1000  # K, 热壁温度
Tc = 500   # K, 冷壁温度
kappa = 1.0  # m^-1, 吸收系数
sigma = 5.670374e-8  # Stefan-Boltzmann常数

# 黑体辐射强度
Ibh = sigma * Th**4 / np.pi
Ibc = sigma * Tc**4 / np.pi

# S4离散坐标 (二维)
# 方向: (xi, eta) 和权重 w
mu_1d = np.array([0.2958759, 0.9082483])
w_1d = np.array([0.5238095, 0.5238095])

# 构建二维方向
xi = []
eta = []
w = []
for i, m1 in enumerate(mu_1d):
    for j, m2 in enumerate(mu_1d):
        xi.append(m1)
        eta.append(m2)
        w.append(w_1d[i] * w_1d[j])
        xi.append(-m1)
        eta.append(m2)
        w.append(w_1d[i] * w_1d[j])
        xi.append(m1)
        eta.append(-m2)
        w.append(w_1d[i] * w_1d[j])
        xi.append(-m1)
        eta.append(-m2)
        w.append(w_1d[i] * w_1d[j])

xi = np.array(xi)
eta = np.array(eta)
w = np.array(w)
w = w / np.sum(w) * 4 * np.pi  # 归一化

print(f"\n物理参数:")
print(f"  方腔边长: {L} m")
print(f"  热壁温度: {Th} K")
print(f"  冷壁温度: {Tc} K")
print(f"  吸收系数: {kappa} m^-1")
print(f"  光学厚度: {kappa * L:.2f}")

# 空间离散
Nx = 30
Ny = 30
dx = L / (Nx - 1)
dy = L / (Ny - 1)
x = np.linspace(0, L, Nx)
y = np.linspace(0, L, Ny)
X, Y = np.meshgrid(x, y)

print(f"\n数值参数:")
print(f"  网格: {Nx} x {Ny}")
print(f"  方向数: {len(xi)}")

# 初始化辐射强度
I = np.zeros((Ny, Nx, len(xi)))

# 假设介质温度分布 (初始猜测)
T = Th - (Th - Tc) * X / L
Ib = sigma * T**4 / np.pi

# 源项迭代
max_iter = 500
tolerance = 1e-5

print("\n开始迭代求解...")

for iteration in range(max_iter):
    I_old = I.copy()
    
    # 计算入射辐射
    G = np.zeros((Ny, Nx))
    for j in range(Ny):
        for i in range(Nx):
            G[j, i] = np.sum(w * I[j, i, :])
    
    # 更新介质温度 (能量方程简化)
    # 假设局部热平衡: G = 4*pi*Ib
    Ib_new = G / (4 * np.pi)
    T_new = (np.pi * Ib_new / sigma)**0.25
    T = 0.5 * T + 0.5 * T_new  # 欠松弛
    Ib = sigma * T**4 / np.pi
    
    # 方向扫描
    for m in range(len(xi)):
        # 根据方向确定扫描顺序
        if xi[m] > 0 and eta[m] > 0:
            i_range = range(Nx)
            j_range = range(Ny)
        elif xi[m] < 0 and eta[m] > 0:
            i_range = range(Nx-1, -1, -1)
            j_range = range(Ny)
        elif xi[m] > 0 and eta[m] < 0:
            i_range = range(Nx)
            j_range = range(Ny-1, -1, -1)
        else:
            i_range = range(Nx-1, -1, -1)
            j_range = range(Ny-1, -1, -1)
        
        for j in j_range:
            for i in i_range:
                # 边界条件
                if i == 0 and xi[m] > 0:
                    I[j, i, m] = Ibh  # 左壁
                    continue
                if i == Nx-1 and xi[m] < 0:
                    I[j, i, m] = Ibc  # 右壁
                    continue
                if j == 0 and eta[m] > 0:
                    I[j, i, m] = I[j+1, i, m] if j+1 < Ny else Ibh  # 下壁 (绝热近似)
                    continue
                if j == Ny-1 and eta[m] < 0:
                    I[j, i, m] = I[j-1, i, m] if j-1 >= 0 else Ibh  # 上壁 (绝热近似)
                    continue
                
                # 迎风格式
                I_upwind_x = I[j, i-1, m] if xi[m] > 0 else I[j, i+1, m]
                I_upwind_y = I[j-1, i, m] if eta[m] > 0 else I[j+1, i, m]
                
                # 简化处理：使用一维近似
                if abs(xi[m]) > abs(eta[m]):
                    I[j, i, m] = (abs(xi[m]) * I_upwind_x + dx * kappa * Ib[j, i]) / \
                                 (abs(xi[m]) + kappa * dx)
                else:
                    I[j, i, m] = (abs(eta[m]) * I_upwind_y + dy * kappa * Ib[j, i]) / \
                                 (abs(eta[m]) + kappa * dy)
    
    # 检查收敛
    if iteration % 50 == 0:
        error = np.max(np.abs(I - I_old))
        print(f"  迭代 {iteration}: 最大误差 = {error:.2e}")
    
    if iteration > 10:
        error = np.max(np.abs(I - I_old))
        if error < tolerance:
            print(f"\n收敛！总迭代次数: {iteration + 1}")
            break

# 计算辐射热流
qx = np.zeros((Ny, Nx))
qy = np.zeros((Ny, Nx))

for j in range(Ny):
    for i in range(Nx):
        qx[j, i] = np.sum(w * xi * I[j, i, :])
        qy[j, i] = np.sum(w * eta * I[j, i, :])

q_magnitude = np.sqrt(qx**2 + qy**2)

print(f"\n计算结果:")
print(f"  最大热流: {np.max(q_magnitude):.2f} W/m²")
print(f"  平均温度: {np.mean(T):.1f} K")

# 可视化
fig, axes = plt.subplots(2, 2, figsize=(14, 12))

# 温度分布
ax1 = axes[0, 0]
im1 = ax1.contourf(X, Y, T, levels=20, cmap='hot')
ax1.set_xlabel('x (m)')
ax1.set_ylabel('y (m)')
ax1.set_title('Temperature Distribution (K)')
plt.colorbar(im1, ax=ax1)

# 辐射强度 (某个方向)
ax2 = axes[0, 1]
im2 = ax2.contourf(X, Y, I[:, :, 0], levels=20, cmap='viridis')
ax2.set_xlabel('x (m)')
ax2.set_ylabel('y (m)')
ax2.set_title(f'Radiation Intensity (dir: ξ={xi[0]:.3f}, η={eta[0]:.3f})')
plt.colorbar(im2, ax=ax2)

# 热流矢量
ax3 = axes[1, 0]
skip = 3
ax3.quiver(X[::skip, ::skip], Y[::skip, ::skip], 
           qx[::skip, ::skip], qy[::skip, ::skip],
           q_magnitude[::skip, ::skip], cmap='plasma', scale=50)
ax3.set_xlabel('x (m)')
ax3.set_ylabel('y (m)')
ax3.set_title('Radiative Heat Flux Vector')
ax3.set_aspect('equal')

# 入射辐射
ax4 = axes[1, 1]
G_final = np.zeros((Ny, Nx))
for j in range(Ny):
    for i in range(Nx):
        G_final[j, i] = np.sum(w * I[j, i, :])
im4 = ax4.contourf(X, Y, G_final, levels=20, cmap='coolwarm')
ax4.set_xlabel('x (m)')
ax4.set_ylabel('y (m)')
ax4.set_title('Incident Radiation (W/m²)')
plt.colorbar(im4, ax=ax4)

plt.tight_layout()
plt.savefig('case2_2d_cavity.png', dpi=150)
plt.close()
print("\n已保存: case2_2d_cavity.png")

# 生成温度分布动画
print("\n生成温度分布动画...")

# 简化动画，只显示温度分布的演变
fig, ax = plt.subplots(figsize=(8, 6))

im = ax.contourf(X, Y, T, levels=20, cmap='hot', vmin=Tc, vmax=Th)
ax.set_xlabel('x (m)')
ax.set_ylabel('y (m)')
ax.set_title('Temperature Distribution Evolution')
cbar = plt.colorbar(im, ax=ax)
cbar.set_label('Temperature (K)')

# 保存静态图作为动画的第一帧
plt.tight_layout()
plt.savefig('case2_temperature_animation.gif', dpi=100)
plt.close()

print("已保存: case2_temperature_animation.gif")

print("\n" + "=" * 70)
print("案例二完成！")
print("=" * 70)

---

9. 案例三：非灰体气体辐射

9.1 问题描述

考虑一维平板间充满非灰体气体（CO₂），压力为1 atm，温度为1000 K。使用WSGGM模型计算气体的有效发射率和辐射热流。

9.2 WSGGM模型

对于CO₂在1000 K，使用4灰气体模型：

灰气体	吸收系数 κk (m⁻¹)	权重系数 ak
1	0.1	0.3
2	1.0	0.4
3	5.0	0.2
4	20.0	0.1

9.3 Python实现

"""
案例三：非灰体气体辐射 - WSGGM模型
"""

import numpy as np
import matplotlib.pyplot as plt
import matplotlib
matplotlib.use('Agg')

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

print("=" * 70)
print("案例三：非灰体气体辐射 - WSGGM模型")
print("=" * 70)

# 物理参数
T_gas = 1000  # K, 气体温度
P_gas = 1.0   # atm, 气体压力
L = 1.0       # m, 几何尺度
sigma = 5.670374e-8  # Stefan-Boltzmann常数

# WSGGM参数 (CO2 at 1000K, 1 atm)
# 4灰气体模型
kappa_k = np.array([0.1, 1.0, 5.0, 20.0])  # m^-1
a_k = np.array([0.3, 0.4, 0.2, 0.1])        # 权重系数

print(f"\n物理参数:")
print(f"  气体温度: {T_gas} K")
print(f"  气体压力: {P_gas} atm")
print(f"  几何尺度: {L} m")

print(f"\nWSGGM参数 (4灰气体):")
for k in range(len(kappa_k)):
    print(f"  灰气体 {k+1}: κ={kappa_k[k]:.1f} m⁻¹, a={a_k[k]:.2f}")

# 计算各灰气体的发射率
epsilon_k = np.zeros(len(kappa_k))
for k in range(len(kappa_k)):
    tau_k = kappa_k[k] * L  # 光学厚度
    epsilon_k[k] = 1 - np.exp(-tau_k)

# 计算总发射率
epsilon_total = np.sum(a_k * epsilon_k)

print(f"\n各灰气体发射率:")
for k in range(len(kappa_k)):
    print(f"  灰气体 {k+1}: ε={epsilon_k[k]:.4f}")

print(f"\n总发射率: {epsilon_total:.4f}")

# 计算辐射热流 (简化模型)
# 假设两壁温度分别为 T1 和 T2
T1 = 1200  # K
T2 = 800   # K

# 各灰气体的辐射热流
q_k = np.zeros(len(kappa_k))
for k in range(len(kappa_k)):
    # 使用Schwarzschild-Schwarzschild近似
    q_k[k] = a_k[k] * sigma * (T1**4 - T2**4) * epsilon_k[k]

q_total = np.sum(q_k)

print(f"\n辐射热流计算:")
print(f"  壁面1温度: {T1} K")
print(f"  壁面2温度: {T2} K")
for k in range(len(kappa_k)):
    print(f"  灰气体 {k+1}: q={q_k[k]:.2f} W/m²")
print(f"  总辐射热流: {q_total:.2f} W/m²")

# 与灰体近似比较
kappa_gray = 1.0  # 灰体吸收系数
epsilon_gray = 1 - np.exp(-kappa_gray * L)
q_gray = sigma * (T1**4 - T2**4) * epsilon_gray

print(f"\n灰体近似比较:")
print(f"  灰体吸收系数: {kappa_gray} m⁻¹")
print(f"  灰体发射率: {epsilon_gray:.4f}")
print(f"  灰体热流: {q_gray:.2f} W/m²")
print(f"  相对误差: {abs(q_gray - q_total)/q_total * 100:.1f}%")

# 可视化
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 灰气体权重分布
ax1 = axes[0, 0]
ax1.bar(range(1, len(a_k)+1), a_k, color='steelblue', alpha=0.7)
ax1.set_xlabel('Gray Gas Index')
ax1.set_ylabel('Weight Coefficient')
ax1.set_title('WSGGM Weight Coefficients')
ax1.set_xticks(range(1, len(a_k)+1))
ax1.grid(True, alpha=0.3)

# 吸收系数分布
ax2 = axes[0, 1]
ax2.bar(range(1, len(kappa_k)+1), kappa_k, color='coral', alpha=0.7)
ax2.set_xlabel('Gray Gas Index')
ax2.set_ylabel('Absorption Coefficient (m⁻¹)')
ax2.set_title('Absorption Coefficients')
ax2.set_xticks(range(1, len(kappa_k)+1))
ax2.set_yscale('log')
ax2.grid(True, alpha=0.3)

# 各灰气体发射率
ax3 = axes[1, 0]
ax3.bar(range(1, len(epsilon_k)+1), epsilon_k, color='seagreen', alpha=0.7)
ax3.set_xlabel('Gray Gas Index')
ax3.set_ylabel('Emissivity')
ax3.set_title('Emissivity of Each Gray Gas')
ax3.set_xticks(range(1, len(epsilon_k)+1))
ax3.grid(True, alpha=0.3)

# 温度-发射率关系
ax4 = axes[1, 1]
T_range = np.linspace(300, 1500, 100)
epsilon_T = np.zeros_like(T_range)

# 简化的温度依赖模型
for i, T in enumerate(T_range):
    # 假设权重随温度变化
    a_T = a_k * (T / 1000)**0.5
    a_T = a_T / np.sum(a_T)  # 归一化
    epsilon_T[i] = np.sum(a_T * epsilon_k)

ax4.plot(T_range, epsilon_T, 'b-', linewidth=2)
ax4.axvline(x=T_gas, color='r', linestyle='--', alpha=0.5, label=f'T={T_gas}K')
ax4.set_xlabel('Temperature (K)')
ax4.set_ylabel('Total Emissivity')
ax4.set_title('Temperature Dependence of Emissivity')
ax4.legend()
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('case3_nongray_gas.png', dpi=150)
plt.close()
print("\n已保存: case3_nongray_gas.png")

# 生成波谱分布动画
print("\n生成波谱分布动画...")

from matplotlib.animation import FuncAnimation, PillowWriter

fig, ax = plt.subplots(figsize=(10, 6))

# 简化的波谱分布
wavelength = np.linspace(1, 20, 100)  # μm

# 不同温度下的波谱发射率
lines = []
colors = plt.cm.plasma(np.linspace(0, 1, 5))

for i, T in enumerate([500, 750, 1000, 1250, 1500]):
    # 简化的波谱模型
    epsilon_spectral = np.zeros_like(wavelength)
    for k in range(len(kappa_k)):
        a_T = a_k[k] * (T / 1000)**0.5
        # 简化的波谱分布
        peak_wl = 10 / kappa_k[k] if kappa_k[k] > 0 else 10
        epsilon_spectral += a_T * np.exp(-((wavelength - peak_wl)/5)**2)
    
    line, = ax.plot(wavelength, epsilon_spectral, color=colors[i], 
                    linewidth=2, label=f'T={T}K')
    lines.append(line)

ax.set_xlabel('Wavelength (μm)')
ax.set_ylabel('Spectral Emissivity')
ax.set_title('Spectral Emissivity Distribution')
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_ylim([0, 1])

plt.tight_layout()
plt.savefig('case3_spectral_distribution.gif', dpi=100)
plt.close()

print("已保存: case3_spectral_distribution.gif")

print("\n" + "=" * 70)
print("案例三完成！")
print("=" * 70)

---

10. 案例四：瞬态辐射传递

10.1 问题描述

考虑一维瞬态辐射传递问题。初始时刻介质温度为 $T_0=300$ K，左壁温度突然升高到 $T_h=1000$ K。研究辐射波的传播过程。

10.2 数学模型

瞬态辐射传递方程：

$$\frac{1}{c}\frac{\partial I}{\partial t}+\mu \frac{\partial I}{\partial x}=-\beta I+S$$

其中 $c$ 是光速。

10.3 Python实现

"""
案例四：瞬态辐射传递
"""

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation, PillowWriter
import matplotlib
matplotlib.use('Agg')

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

print("=" * 70)
print("案例四：瞬态辐射传递")
print("=" * 70)

# 物理参数
L = 1.0       # m, 板间距
T0 = 300      # K, 初始温度
Th = 1000     # K, 热壁温度
kappa = 1.0   # m^-1, 吸收系数
c = 3e8       # m/s, 光速
sigma = 5.670374e-8  # Stefan-Boltzmann常数

# 数值参数
Nx = 50
dx = L / (Nx - 1)
x = np.linspace(0, L, Nx)

# 时间参数
t_final = 1e-7  # s (辐射传播时间尺度)
dt = dx / c * 0.5  # CFL条件
Nt = int(t_final / dt)

print(f"\n物理参数:")
print(f"  板间距: {L} m")
print(f"  初始温度: {T0} K")
print(f"  热壁温度: {Th} K")
print(f"  吸收系数: {kappa} m⁻¹")

print(f"\n数值参数:")
print(f"  网格数: {Nx}")
print(f"  时间步长: {dt:.2e} s")
print(f"  总步数: {Nt}")

# S2离散坐标
mu = np.array([-0.5773503, 0.5773503])
w = np.array([np.pi, np.pi])

# 初始化
I = np.zeros((Nx, len(mu)))
T = np.ones(Nx) * T0

# 存储历史数据
I_history = []
T_history = []
t_history = []

print("\n开始瞬态计算...")

for n in range(Nt):
    t = n * dt
    
    # 保存数据
    if n % max(1, Nt//100) == 0:
        I_history.append(I.copy())
        T_history.append(T.copy())
        t_history.append(t)
    
    # 更新辐射强度 (显式格式)
    I_new = I.copy()
    
    for m in range(len(mu)):
        Ib = sigma * T**4 / np.pi
        S = kappa * Ib
        
        if mu[m] > 0:
            # 左边界
            I_new[0, m] = sigma * Th**4 / np.pi
            # 内部
            for i in range(1, Nx):
                I_new[i, m] = I[i, m] - c * dt * (
                    mu[m] * (I[i, m] - I[i-1, m]) / dx +
                    kappa * I[i, m] - S[i]
                )
        else:
            # 右边界
            I_new[-1, m] = sigma * T0**4 / np.pi
            # 内部
            for i in range(Nx-2, -1, -1):
                I_new[i, m] = I[i, m] - c * dt * (
                    mu[m] * (I[i+1, m] - I[i, m]) / dx +
                    kappa * I[i, m] - S[i]
                )
    
    I = I_new
    
    # 更新温度 (简化能量方程)
    G = np.zeros(Nx)
    for i in range(Nx):
        G[i] = np.sum(w * I[i, :])
    
    # 假设热容很小，温度快速达到平衡
    T = (G / (4 * sigma))**0.25

print(f"计算完成！保存了 {len(I_history)} 个时间步")

# 可视化
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 最终辐射强度分布
ax1 = axes[0, 0]
for m in range(len(mu)):
    ax1.plot(x, I_history[-1][:, m], linewidth=2, label=f'μ={mu[m]:.3f}')
ax1.set_xlabel('Position (m)')
ax1.set_ylabel('Radiation Intensity (W/(m²·sr))')
ax1.set_title('Final Radiation Intensity Distribution')
ax1.legend()
ax1.grid(True, alpha=0.3)

# 最终温度分布
ax2 = axes[0, 1]
ax2.plot(x, T_history[-1], 'r-', linewidth=2)
ax2.axhline(y=Th, color='r', linestyle='--', alpha=0.5, label=f'Th={Th}K')
ax2.axhline(y=T0, color='b', linestyle='--', alpha=0.5, label=f'T0={T0}K')
ax2.set_xlabel('Position (m)')
ax2.set_ylabel('Temperature (K)')
ax2.set_title('Final Temperature Distribution')
ax2.legend()
ax2.grid(True, alpha=0.3)

# 温度演化
ax3 = axes[1, 0]
for i in range(0, len(t_history), max(1, len(t_history)//10)):
    ax3.plot(x, T_history[i], alpha=0.7, label=f't={t_history[i]:.2e}s')
ax3.set_xlabel('Position (m)')
ax3.set_ylabel('Temperature (K)')
ax3.set_title('Temperature Evolution')
ax3.legend(fontsize=8)
ax3.grid(True, alpha=0.3)

# 热壁热流随时间变化
ax4 = axes[1, 1]
q_wall = []
for I_t in I_history:
    q = 2 * np.pi * np.sum(w * mu * I_t[0, :])
    q_wall.append(q)

ax4.plot(np.array(t_history) * 1e9, q_wall, 'b-', linewidth=2)
ax4.set_xlabel('Time (ns)')
ax4.set_ylabel('Heat Flux at Hot Wall (W/m²)')
ax4.set_title('Heat Flux Evolution')
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('case4_transient_radiation.png', dpi=150)
plt.close()
print("\n已保存: case4_transient_radiation.png")

# 生成瞬态演化动画
print("\n生成瞬态演化动画...")

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

# 辐射强度
lines1 = []
for m in range(len(mu)):
    line, = ax1.plot([], [], linewidth=2, label=f'μ={mu[m]:.3f}')
    lines1.append(line)
ax1.set_xlim([0, L])
ax1.set_ylim([0, np.max(I_history[-1]) * 1.1])
ax1.set_xlabel('Position (m)')
ax1.set_ylabel('Radiation Intensity (W/(m²·sr))')
ax1.set_title('Radiation Intensity Evolution')
ax1.legend()
ax1.grid(True, alpha=0.3)

# 温度
line2, = ax2.plot([], [], 'r-', linewidth=2)
ax2.set_xlim([0, L])
ax2.set_ylim([T0 * 0.9, Th * 1.1])
ax2.set_xlabel('Position (m)')
ax2.set_ylabel('Temperature (K)')
ax2.set_title('Temperature Evolution')
ax2.axhline(y=Th, color='r', linestyle='--', alpha=0.3)
ax2.axhline(y=T0, color='b', linestyle='--', alpha=0.3)
ax2.grid(True, alpha=0.3)

text = fig.text(0.5, 0.95, '', ha='center', fontsize=12)

def init():
    for line in lines1:
        line.set_data([], [])
    line2.set_data([], [])
    text.set_text('')
    return lines1 + [line2, text]

def update(frame):
    I_frame = I_history[frame]
    T_frame = T_history[frame]
    
    for m, line in enumerate(lines1):
        line.set_data(x, I_frame[:, m])
    
    line2.set_data(x, T_frame)
    
    text.set_text(f'Time: {t_history[frame]*1e9:.2f} ns')
    return lines1 + [line2, text]

ani = FuncAnimation(fig, update, frames=len(I_history), 
                    init_func=init, interval=100, blit=False)
writer = PillowWriter(fps=10)
ani.save('case4_transient_animation.gif', writer=writer)
plt.close()
print("已保存: case4_transient_animation.gif")

print("\n" + "=" * 70)