辐射换热仿真-主题020_光学厚介质辐射换热-主题020_光学厚介质辐射换热

kkchenjj

61人浏览 · 2026-04-24 21:46:15

kkchenjj · 2026-04-24 21:46:15 发布

主题020：光学厚介质辐射换热

摘要

光学厚介质是辐射换热中的另一重要极端情况，当介质的光学厚度远大于1时，辐射传递呈现出扩散特性。本教程深入分析光学厚介质的辐射特性，详细推导罗斯兰扩散近似（Rosseland Approximation），建立光学厚极限下的辐射导热模型。通过五个典型应用案例——光学厚平板辐射换热、球形容器辐射传热、高温炉膛辐射导热、辐射与导热耦合传热以及非均匀光学厚介质辐射——全面展示光学厚介质辐射换热的理论方法和工程应用。所有案例均配有完整的Python仿真程序和可视化分析，为工程实践提供可靠的理论指导。

关键词： 光学厚、罗斯兰近似、扩散模型、辐射导热、有效导热系数、辐射传递

在这里插入图片描述

1. 光学厚介质的基本概念

1.1 光学厚度的定义与物理意义

光学厚度（Optical Thickness）是描述介质对辐射吸收能力的无量纲参数，定义为：

$τ=∫0Lκλ(s)ds\tau = \int_0^L \kappa_\lambda(s) ds$

其中， $κλ\kappa_\lambda$ 是光谱吸收系数， $L$ 是辐射传播路径长度。光学厚度表征了辐射在穿过介质时被吸收的程度：

当 $τ≪1\tau \ll 1$ 时，介质为光学薄，辐射几乎不被吸收
当 $τ≈1\tau \approx 1$ 时，介质为光学中等厚度
当 $τ≫1\tau \gg 1$ 时，介质为光学厚，辐射被强烈吸收

光学厚条件意味着介质对辐射几乎是"不透明"的，辐射在穿过介质很短距离后就会被完全吸收。这一特性在高温炉膛、恒星内部、核反应堆等场景中具有重要应用价值。

1.2 光学厚极限的物理图像

在光学厚极限下（ $τ→∞\tau \to \infty$ ），辐射传递过程呈现出独特的物理特性：

1. 局部热力学平衡

在光学厚介质中，辐射与物质达到局部热力学平衡（Local Thermodynamic Equilibrium, LTE）。此时，辐射强度接近局部黑体辐射强度：

$Iλ≈Ibλ(T)I_\lambda \approx I_{b\lambda}(T)$

2. 扩散特性

辐射在光学厚介质中的传播类似于热传导过程，呈现出扩散特性。辐射能量通过多次吸收和再发射逐步传递，而不是直接穿透介质。

3. 辐射热流与温度梯度成正比

类似于傅里叶导热定律，光学厚介质中的辐射热流与温度梯度成正比：

$qr=−kr∇T\mathbf{q}_r = -k_r \nabla T$

其中， $k_r$ 是辐射导热系数。

4. 边界层效应

在光学厚介质与边界壁面交界处，存在一个薄层（边界层），在该层内辐射强度从壁面值逐渐过渡到介质内部的黑体辐射强度。边界层厚度约为几个平均自由程。

1.3 光学厚介质辐射的特征

光学厚介质辐射具有以下显著特征：

1. 强吸收特性

辐射穿过光学厚介质时，透射率接近0：

$τtrans=e−τ≈0\tau_{trans} = e^{-\tau} \approx 0$

这意味着入射辐射几乎被完全吸收，无法穿透介质。

2. 发射主导特性

在光学厚极限下，介质的辐射主要由发射贡献。介质发射的辐射强度接近黑体辐射强度：

$Iλemission≈Ibλ(T)I_\lambda^{emission} \approx I_{b\lambda}(T)$

3. 各向同性特性

在光学厚介质内部，辐射强度接近各向同性，即辐射强度在所有方向上近似相等。

4. 温度梯度驱动

辐射热流由温度梯度驱动，类似于导热过程。温度梯度越大，辐射热流越大。

1.4 光学厚条件的工程应用

光学厚近似在以下工程领域有广泛应用：

高温炉膛

工业炉膛、玻璃熔炉等设备中，高温气体和颗粒物形成光学厚介质。罗斯兰近似可用于快速估算炉膛内的辐射传热。

恒星内部

恒星内部等离子体处于极端高温高压状态，光学厚度极大。辐射扩散是恒星能量传输的主要机制。

核反应堆

核反应堆堆芯中，高温冷却剂和结构材料形成光学厚环境。辐射传热分析对于反应堆安全设计至关重要。

大气科学

地球大气层在某些波段（如红外波段）呈现光学厚特性。辐射扩散模型用于气候模拟和天气预报。

2. 罗斯兰扩散近似

2.1 辐射传递方程的扩散近似

辐射传递方程（RTE）为：

$dIλds=−κλIλ+κλIbλ\frac{dI_\lambda}{ds} = -\kappa_\lambda I_\lambda + \kappa_\lambda I_{b\lambda}$

在光学厚极限下，可以假设辐射强度接近黑体辐射强度，但存在小的偏离：

$Iλ=Ibλ+Iλ(1)I_\lambda = I_{b\lambda} + I_\lambda^{(1)}$

其中， $Iλ(1)I_\lambda^{(1)}$ 是一阶修正项， $∣Iλ(1)∣≪Ibλ|I_\lambda^{(1)}| \ll I_{b\lambda}$ 。

将上式代入RTE，得到：

$dIbλds+dIλ(1)ds=−κλIλ(1)\frac{dI_{b\lambda}}{ds} + \frac{dI_\lambda^{(1)}}{ds} = -\kappa_\lambda I_\lambda^{(1)}$

由于 $Iλ(1)I_\lambda^{(1)}$ 是小量，可以忽略其二阶导数项：

$dIbλds≈−κλIλ(1)\frac{dI_{b\lambda}}{ds} \approx -\kappa_\lambda I_\lambda^{(1)}$

因此：

$Iλ(1)=−1κλdIbλdsI_\lambda^{(1)} = -\frac{1}{\kappa_\lambda}\frac{dI_{b\lambda}}{ds}$

2.2 罗斯兰平均吸收系数

辐射热流为：

$qr=∫0∞∫4πIλsdΩdλ\mathbf{q}_r = \int_0^\infty \int_{4\pi} I_\lambda \mathbf{s} d\Omega d\lambda$

在光学厚极限下，主要贡献来自一阶修正项：

$qr=∫0∞∫4πIλ(1)sdΩdλ\mathbf{q}_r = \int_0^\infty \int_{4\pi} I_\lambda^{(1)} \mathbf{s} d\Omega d\lambda$

将 $Iλ(1)I_\lambda^{(1)}$ 的表达式代入：

$qr=−∫0∞1κλ∫4πdIbλdssdΩdλ\mathbf{q}_r = -\int_0^\infty \frac{1}{\kappa_\lambda} \int_{4\pi} \frac{dI_{b\lambda}}{ds} \mathbf{s} d\Omega d\lambda$

利用链式法则：

$dIbλds=dIbλdTdTds=dIbλdT(∇T⋅s)\frac{dI_{b\lambda}}{ds} = \frac{dI_{b\lambda}}{dT} \frac{dT}{ds} = \frac{dI_{b\lambda}}{dT} (\nabla T \cdot \mathbf{s})$

因此：

$qr=−∫0∞1κλdIbλdT∫4π(∇T⋅s)sdΩdλ\mathbf{q}_r = -\int_0^\infty \frac{1}{\kappa_\lambda} \frac{dI_{b\lambda}}{dT} \int_{4\pi} (\nabla T \cdot \mathbf{s}) \mathbf{s} d\Omega d\lambda$

计算角度积分：

$∫4π(∇T⋅s)sdΩ=4π3∇T\int_{4\pi} (\nabla T \cdot \mathbf{s}) \mathbf{s} d\Omega = \frac{4\pi}{3} \nabla T$

因此：

$qr=−4π3∇T∫0∞1κλdIbλdTdλ\mathbf{q}_r = -\frac{4\pi}{3} \nabla T \int_0^\infty \frac{1}{\kappa_\lambda} \frac{dI_{b\lambda}}{dT} d\lambda$

定义罗斯兰平均吸收系数：

$1κR=∫0∞1κλdIbλdTdλ∫0∞dIbλdTdλ\frac{1}{\kappa_R} = \frac{\int_0^\infty \frac{1}{\kappa_\lambda} \frac{dI_{b\lambda}}{dT} d\lambda}{\int_0^\infty \frac{dI_{b\lambda}}{dT} d\lambda}$

利用普朗克定律：

$∫0∞dIbλdTdλ=ddT(σT4π)=4σT3π\int_0^\infty \frac{dI_{b\lambda}}{dT} d\lambda = \frac{d}{dT}\left(\frac{\sigma T^4}{\pi}\right) = \frac{4\sigma T^3}{\pi}$

因此，辐射热流可以表示为：

$qr=−16σT33κR∇T\mathbf{q}_r = -\frac{16\sigma T^3}{3\kappa_R} \nabla T$

2.3 辐射导热系数

定义辐射导热系数：

$kr=16σT33κRk_r = \frac{16\sigma T^3}{3\kappa_R}$

辐射热流可以写成类似于傅里叶导热定律的形式：

$qr=−kr∇T\mathbf{q}_r = -k_r \nabla T$

辐射导热系数与温度的三次方成正比，与罗斯兰平均吸收系数成反比。在高温条件下，辐射导热系数可以远大于材料的本征导热系数。

2.4 有效导热系数

在同时存在导热和辐射的情况下，可以定义有效导热系数：

$keff=kc+kr=kc+16σT33κRk_{eff} = k_c + k_r = k_c + \frac{16\sigma T^3}{3\kappa_R}$

其中， $k_c$ 是材料的本征导热系数（导热贡献）。

总热流为：

$q=−(kc+kr)∇T=−keff∇T\mathbf{q} = -(k_c + k_r) \nabla T = -k_{eff} \nabla T$

2.5 边界条件处理

在光学厚介质与壁面交界处，需要特殊的边界条件处理。由于边界层效应，壁面处的辐射热流不等于介质内部的扩散热流。

跳跃边界条件

在光学厚极限下，可以采用跳跃边界条件（Jump Boundary Condition）：

$Twall−Tgas=23κR(∂T∂n)wallT_{wall} - T_{gas} = \frac{2}{3\kappa_R} \left(\frac{\partial T}{\partial n}\right)_{wall}$

其中， $T_{wall}$ 是壁面温度， $T_{gas}$ 是紧邻壁面的气体温度， $∂T∂n\frac{\partial T}{\partial n}$ 是温度在壁面法向的梯度。

这一边界条件反映了辐射在边界层内的非扩散行为。

** slip 边界条件**

对于灰体介质，可以采用更精确的 slip 边界条件：

$Twall−Tgas=23κR(∂T∂n)wall⋅2−εεT_{wall} - T_{gas} = \frac{2}{3\kappa_R} \left(\frac{\partial T}{\partial n}\right)_{wall} \cdot \frac{2 - \varepsilon}{\varepsilon}$

其中， $ε\varepsilon$ 是壁面发射率。

3. 光学厚介质辐射的数值方法

3.1 有限差分法

对于一维光学厚介质，能量方程为：

$ddz(keffdTdz)=0\frac{d}{dz}\left(k_{eff} \frac{dT}{dz}\right) = 0$

采用中心差分格式：

$ki+1/2(Ti+1−Ti)−ki−1/2(Ti−Ti−1)(Δz)2=0\frac{k_{i+1/2}(T_{i+1} - T_i) - k_{i-1/2}(T_i - T_{i-1})}{(\Delta z)^2} = 0$

其中， $ki+1/2=ki+ki+12k_{i+1/2} = \frac{k_i + k_{i+1}}{2}$ 是界面导热系数。

3.2 有限元法

对于复杂几何，可以采用有限元法求解辐射导热问题。弱形式为：

$∫Ωkeff∇T⋅∇vdΩ=∫ΓqbvdΓ\int_\Omega k_{eff} \nabla T \cdot \nabla v d\Omega = \int_{\Gamma} q_b v d\Gamma$

其中， $v$ 是测试函数， $q_b$ 是边界热流。

3.3 谱方法

对于高精度需求，可以采用谱方法求解。将温度展开为基函数的线性组合：

$\sum_{n=0}^N a_n \phi_n(z)$

其中， $ϕn(z)\phi_n(z)$ 可以是切比雪夫多项式、勒让德多项式等正交基函数。

3.4 蒙特卡洛法

对于非均匀介质或复杂边界条件，可以采用蒙特卡洛方法。但由于光学厚介质中光子平均自由程很短，蒙特卡洛方法的计算效率较低。

4. 应用案例

4.1 案例一：光学厚平板辐射换热

问题描述

考虑两个无限大平行平板，其间填充光学厚气体。上平板温度 $T_1 = 2000$ K，下平板温度 $T_2 = 1500$ K。平板间距 $L = 0.5$ m，气体吸收系数 $κ=10\kappa = 10$ m⁻¹。计算两平板间的辐射换热，比较罗斯兰近似与精确解的差异。

数学模型

对于一维光学厚介质，能量方程为：

$ddz(krdTdz)=0\frac{d}{dz}\left(k_r \frac{dT}{dz}\right) = 0$

其中，辐射导热系数为：

$kr=16σT33κk_r = \frac{16\sigma T^3}{3\kappa}$

罗斯兰近似解

假设辐射导热系数为常数（取平均温度下的值），则温度分布为线性：

$T_1 + (T_2 - T_1)\frac{z}{L}$

辐射热流为：

$qr=−krT2−T1L=16σTavg33κT1−T2Lq_r = -k_r \frac{T_2 - T_1}{L} = \frac{16\sigma T_{avg}^3}{3\kappa} \frac{T_1 - T_2}{L}$

其中， $T_{avg} = (T_1 + T_2)/2$ 是平均温度。

考虑温度变化的精确解

考虑辐射导热系数随温度变化，能量方程变为非线性：

$ddz(16σT33κdTdz)=0\frac{d}{dz}\left(\frac{16\sigma T^3}{3\kappa} \frac{dT}{dz}\right) = 0$

积分得到：

$T4=T14+(T24−T14)zLT^4 = T_1^4 + (T_2^4 - T_1^4)\frac{z}{L}$

辐射热流为：

$qr=4σ3κL(T14−T24)q_r = \frac{4\sigma}{3\kappa L}(T_1^4 - T_2^4)$

仿真结果与分析

图1展示了光学厚平板辐射换热的温度分布和热流分布。从图中可以看出：

温度分布近似为线性，但在高温区域略有偏差
罗斯兰近似与精确解吻合良好
辐射热流与温度差的四次方成正比

工程意义

本案例展示了罗斯兰近似在工程计算中的应用价值。对于高温光学厚介质，罗斯兰近似可以大幅简化计算，同时保持足够的精度。

4.2 案例二：球形容器辐射传热

问题描述

考虑一个球形容器，内半径 $R_1 = 0.2$ m，外半径 $R_2 = 0.5$ m。容器内填充光学厚气体，内壁温度 $T_1 = 2200$ K，外壁温度 $T_2 = 1800$ K。气体吸收系数 $κ=15\kappa = 15$ m⁻¹。计算球形容器内的温度分布和辐射热流。

数学模型

对于球坐标系，能量方程为：

$1r2ddr(r2krdTdr)=0\frac{1}{r^2}\frac{d}{dr}\left(r^2 k_r \frac{dT}{dr}\right) = 0$

罗斯兰近似解

假设辐射导热系数为常数，温度分布为：

$T_1 + (T_2 - T_1)\frac{1/r - 1/R_1}{1/R_2 - 1/R_1}$

辐射热流为：

$qr=−krdTdr=kr(T1−T2)r2(1/R1−1/R2)q_r = -k_r \frac{dT}{dr} = \frac{k_r(T_1 - T_2)}{r^2(1/R_1 - 1/R_2)}$

总热流量为：

$4\pi r^2 q_r = \frac{4\pi k_r(T_1 - T_2)}{1/R_1 - 1/R_2}$

考虑温度变化的精确解

考虑辐射导热系数随温度变化：

$1r2ddr(r216σT33κdTdr)=0\frac{1}{r^2}\frac{d}{dr}\left(r^2 \frac{16\sigma T^3}{3\kappa} \frac{dT}{dr}\right) = 0$

积分得到：

$T4=T14+(T24−T14)1/r−1/R11/R2−1/R1T^4 = T_1^4 + (T_2^4 - T_1^4)\frac{1/r - 1/R_1}{1/R_2 - 1/R_1}$

仿真结果与分析

图2展示了球形容器辐射传热的温度分布和热流分布。结果表明：

温度分布呈非线性，在靠近内壁处温度梯度较大
辐射热流与半径的平方成反比
总热流量与几何形状有关

工程应用

球形容器模型广泛应用于：

核反应堆压力容器的热分析
高温储热罐的设计
行星内部热传输研究

4.3 案例三：高温炉膛辐射导热

问题描述

考虑一个工业炉膛，炉膛尺寸为 $2m \times 2m \times 4m$ 。炉膛内充满高温燃烧气体，温度分布不均匀。炉膛壁面温度为 $T_w = 1500$ K，炉膛中心温度为 $T_c = 2000$ K。气体吸收系数 $κ=8\kappa = 8$ m⁻¹。采用罗斯兰近似计算炉膛内的温度分布和辐射热流。

数学模型

对于三维光学厚介质，能量方程为：

$∇⋅(kr∇T)=0\nabla \cdot (k_r \nabla T) = 0$

数值求解

采用有限差分法求解。将炉膛离散为 $Nx×Ny×NzN_x \times N_y \times N_z$ 个网格单元。

离散方程为：

$ki+1/2(Ti+1−Ti)−ki−1/2(Ti−Ti−1)(Δx)2+...=0\frac{k_{i+1/2}(T_{i+1} - T_i) - k_{i-1/2}(T_i - T_{i-1})}{(\Delta x)^2} + ... = 0$

边界条件

在炉膛壁面处，采用跳跃边界条件：

$Twall−Tgas=23κ∂T∂nT_{wall} - T_{gas} = \frac{2}{3\kappa} \frac{\partial T}{\partial n}$

仿真结果与分析

图3展示了高温炉膛辐射导热的温度分布和热流分布。结果表明：

温度从炉膛中心向壁面逐渐降低
温度梯度在壁面附近最大
辐射热流主要沿径向传播

工程指导

本案例为工业炉膛设计提供了重要指导：

炉膛壁面需要耐高温材料
温度梯度大的区域需要加强冷却
辐射传热是炉膛内主要的传热方式

4.4 案例四：辐射与导热耦合传热

问题描述

考虑一个复合墙体，由两层材料组成：外层为耐火砖（厚度 $L_1 = 0.3$ m，导热系数 $k_1 = 1.5$ W/(m·K)），内层为隔热材料（厚度 $L_2 = 0.2$ m，导热系数 $k_2 = 0.1$ W/(m·K)）。内表面温度 $T_i = 1800$ K，外表面温度 $T_o = 300$ K。在高温侧，辐射传热不可忽略，气体吸收系数 $κ=5\kappa = 5$ m⁻¹。计算通过墙体的总热流。

数学模型

对于耦合传热问题，总热流为导热和辐射之和：

$q_c + q_r = -k_{eff} \frac{dT}{dz}$

其中，有效导热系数为：

$keff=kc+16σT33κk_{eff} = k_c + \frac{16\sigma T^3}{3\kappa}$

分层求解

对于多层结构，总热阻为各层热阻之和：

$Rtotal=∑iLikeff,iR_{total} = \sum_i \frac{L_i}{k_{eff,i}}$

总热流为：

$\frac{T_i - T_o}{R_{total}}$

仿真结果与分析

图4展示了辐射与导热耦合传热的温度分布和热流分布。结果表明：

在高温区域，辐射传热占主导地位
有效导热系数随温度变化显著
隔热层显著降低了热流

工程应用

辐射与导热耦合模型应用于：

高温炉墙设计
航天器热防护系统
核反应堆安全壳

4.5 案例五：非均匀光学厚介质辐射

问题描述

考虑一个光学厚介质层，吸收系数随温度变化： $κ=κ0(T/T0)n\kappa = \kappa_0 (T/T_0)^n$ ，其中 $κ0=10\kappa_0 = 10$ m⁻¹， $T_0 = 1500$ K， $n = - 1.5$ （表示吸收系数随温度升高而降低）。介质层厚度 $L = 0.4$ m，上表面温度 $T_1 = 2000$ K，下表面温度 $T_2 = 1400$ K。计算温度分布和辐射热流。

数学模型

对于非均匀吸收系数，能量方程为：

$ddz(16σT33κ(T)dTdz)=0\frac{d}{dz}\left(\frac{16\sigma T^3}{3\kappa(T)} \frac{dT}{dz}\right) = 0$

数值求解

采用迭代法求解。首先假设温度分布，计算吸收系数，然后求解能量方程，更新温度分布，直到收敛。

仿真结果与分析

图5展示了非均匀光学厚介质的温度分布和热流分布。结果表明：

温度分布呈非线性，与常吸收系数情况不同
吸收系数随温度变化导致辐射导热系数变化
高温区域辐射传热更强

工程意义

非均匀介质模型应用于：

燃烧产物辐射（吸收系数随温度和组分变化）
等离子体辐射
多孔介质辐射

5. 光学厚介质辐射的工程应用

5.1 工业炉膛设计

工业炉膛是光学厚介质辐射的典型应用场景。炉膛内的高温气体和颗粒物形成强吸收环境。

炉膛热分析

采用罗斯兰近似，炉膛内的辐射热流为：

$qr=−16σT33κ∇Tq_r = -\frac{16\sigma T^3}{3\kappa} \nabla T$

炉膛壁面热负荷为：

$qw=4σ3κTg4−Tw4δq_w = \frac{4\sigma}{3\kappa} \frac{T_g^4 - T_w^4}{\delta}$

其中， $δ\delta$ 是边界层厚度。

设计考虑

选择合适的炉膛尺寸和形状
优化燃烧器布置，控制温度分布
采用耐高温材料和冷却系统

5.2 玻璃熔炉

玻璃熔炉中，熔融玻璃是强吸收介质，光学厚度很大。

玻璃辐射特性

熔融玻璃的吸收系数在红外波段很高，呈现光学厚特性。辐射是玻璃熔炉内主要的传热方式。

温度控制

玻璃熔炉需要精确控制温度分布，以保证玻璃质量。罗斯兰近似可用于快速估算温度场。

5.3 核反应堆

核反应堆堆芯中，冷却剂和结构材料形成光学厚环境。

辐射传热分析

在事故工况下，辐射传热可能成为主要的传热机制。光学厚近似可用于安全分析。

设计准则

保证足够的冷却能力
控制燃料温度，防止熔化
考虑辐射传热对结构材料的影响

5.4 恒星内部

恒星内部是极端的光学厚环境，光学厚度可达 $10^8$ 以上。

能量传输

在恒星内部，辐射扩散是主要的能量传输机制。能量从核心向外传递，最终从表面辐射出去。

辐射压

在极端光学厚条件下，辐射压对恒星结构有重要影响：

$Pr=4σ3cT4P_r = \frac{4\sigma}{3c} T^4$

其中， $c$ 是光速。

6. 罗斯兰近似的精度与适用范围

6.1 误差分析

罗斯兰近似的误差主要来源于：

扩散假设误差：假设辐射强度接近各向同性，在边界层内不成立
温度梯度假设：假设温度梯度较小，在温度剧变区域不成立
灰体假设误差：假设吸收系数与波长无关，对于非灰体介质不精确

6.2 适用范围判据

罗斯兰近似的适用范围：

光学厚度 $τ\tau$	适用性
$> 10$	极好
$\sim 10$	良好
$\sim 5$	中等（需要修正）
$< 1$	不适用

6.3 改进方法

对于中等光学厚度，可以采用改进的扩散近似：

P-1近似

P-1近似是球谐函数法的一阶近似，适用于中等光学厚度：

$∇⋅(13κ∇G)=κ(G−4πIb)\nabla \cdot \left(\frac{1}{3\kappa} \nabla G\right) = \kappa(G - 4\pi I_b)$

其中， $\int_{4\pi} I d\Omega$ 是入射辐射。

通量限制器

在温度梯度大的区域，可以采用通量限制器（Flux Limiter）修正扩散近似：

$qr=−kr1+∣∇T∣Tκ∇T\mathbf{q}_r = -\frac{k_r}{1 + \frac{|\nabla T|}{T \kappa}} \nabla T$

6.4 数值验证

通过数值计算验证罗斯兰近似的精度。考虑一维平板间辐射换热，比较精确解、罗斯兰近似和改进近似的计算结果。

验证结果

罗斯兰近似在 $τ>10\tau > 10$ 时误差小于1%
P-1近似在 $τ>1\tau > 1$ 时误差小于5%
通量限制器可以改善大温度梯度区域的精度

7. 光学厚介质辐射的高级主题

7.1 非灰体光学厚介质

实际介质通常具有波长依赖的吸收特性（非灰体）。对于非灰体光学厚介质，需要对波谱进行积分。

波谱平均

罗斯兰平均吸收系数已经对波谱进行了加权平均：

$1κR=∫0∞1κλdIbλdTdλ∫0∞dIbλdTdλ\frac{1}{\kappa_R} = \frac{\int_0^\infty \frac{1}{\kappa_\lambda} \frac{dI_{b\lambda}}{dT} d\lambda}{\int_0^\infty \frac{dI_{b\lambda}}{dT} d\lambda}$

普朗克平均与罗斯兰平均

普朗克平均吸收系数（用于发射计算）：

$κP=∫0∞κλIbλdλ∫0∞Ibλdλ\kappa_P = \frac{\int_0^\infty \kappa_\lambda I_{b\lambda} d\lambda}{\int_0^\infty I_{b\lambda} d\lambda}$

罗斯兰平均吸收系数（用于传输计算）：

$1κR=∫0∞1κλdIbλdTdλ∫0∞dIbλdTdλ\frac{1}{\kappa_R} = \frac{\int_0^\infty \frac{1}{\kappa_\lambda} \frac{dI_{b\lambda}}{dT} d\lambda}{\int_0^\infty \frac{dI_{b\lambda}}{dT} d\lambda}$

两种平均方法适用于不同的物理过程。

7.2 各向异性散射

光学厚介质中的散射可能是各向异性的。在光学厚极限下，散射对辐射传输的影响可以忽略，因为辐射在传播很短距离后就会被吸收。

散射修正

对于中等光学厚度，可以考虑散射修正：

$κeff=κa+κs(1−g)\kappa_{eff} = \kappa_a + \kappa_s(1 - g)$

其中， $κa\kappa_a$ 是吸收系数， $κs\kappa_s$ 是散射系数， $g$ 是不对称因子。

7.3 瞬态光学厚辐射

对于瞬态问题，需要考虑辐射传播的时间效应。

瞬态扩散方程

$ρcp∂T∂t=∇⋅(keff∇T)\rho c_p \frac{\partial T}{\partial t} = \nabla \cdot (k_{eff} \nabla T)$

其中， $ρ\rho$ 是密度， $c_p$ 是比热容。

特征时间

辐射扩散的特征时间为：

$τrad=ρcpL2kr\tau_{rad} = \frac{\rho c_p L^2}{k_r}$

7.4 多维光学厚辐射

实际工程问题通常涉及多维几何。在光学厚极限下，多维辐射问题可以简化为扩散问题。

二维/三维扩散方程

$ρcp∂T∂t=∂∂x(keff∂T∂x)+∂∂y(keff∂T∂y)+∂∂z(keff∂T∂z)\rho c_p \frac{\partial T}{\partial t} = \frac{\partial}{\partial x}\left(k_{eff} \frac{\partial T}{\partial x}\right) + \frac{\partial}{\partial y}\left(k_{eff} \frac{\partial T}{\partial y}\right) + \frac{\partial}{\partial z}\left(k_{eff} \frac{\partial T}{\partial z}\right)$

边界条件

在复杂几何边界上，需要采用适当的边界条件：

第一类边界条件：给定温度
第二类边界条件：给定热流
第三类边界条件：对流换热
跳跃边界条件：辐射扩散边界

8. 总结与展望

8.1 主要内容回顾

本教程系统研究了光学厚介质的辐射换热特性，主要内容包括：

光学厚介质的基本概念：介绍了光学厚度的定义、光学厚极限的物理图像和辐射特征。
罗斯兰扩散近似：详细推导了罗斯兰近似，定义了辐射导热系数和有效导热系数。
数值方法：介绍了有限差分法、有限元法、谱方法和蒙特卡洛法等求解光学厚辐射问题的方法。
应用案例：通过五个典型案例，展示了罗斯兰近似在平板、球体、炉膛、耦合传热和非均匀介质中的应用。
工程应用：讨论了光学厚模型在工业炉膛、玻璃熔炉、核反应堆和恒星内部等领域的应用。
精度分析：分析了罗斯兰近似的误差来源和适用范围，提出了改进方法。

8.2 关键结论

通过本教程的学习，可以得出以下关键结论：

罗斯兰近似适用条件：当光学厚度 $τ>10\tau > 10$ 时，罗斯兰近似精度很高；当 $\tau < 10$ 时，可以采用P-1近似等改进方法。
辐射导热系数：光学厚介质的辐射导热系数与温度的三次方成正比，与吸收系数成反比。
有效导热系数：在同时存在导热和辐射的情况下，有效导热系数等于两者之和。
工程应用价值：罗斯兰模型在高温光学厚介质（如炉膛、恒星内部）中有广泛应用，可以大幅简化计算。

8.3 未来发展方向

光学厚介质辐射换热研究的发展方向包括：

非灰体模型：发展更精确的非灰体光学厚模型，考虑波谱依赖的吸收特性。
多物理场耦合：研究光学厚辐射与流动、化学反应的耦合效应。
边界层模型：发展更精确的边界层模型，改善边界附近的计算精度。
机器学习应用：利用机器学习方法加速光学厚辐射计算，提高工程应用效率。
实验验证：开展更多的实验研究，验证罗斯兰模型的精度和适用范围。

8.4 学习建议

对于希望深入学习光学厚介质辐射换热的读者，建议：

理论基础：系统学习辐射传递方程、罗斯兰近似、扩散理论等基本概念。
数值方法：掌握有限差分法、有限元法等数值求解方法。
编程实践：通过编写Python程序，实现光学厚辐射的仿真计算。
工程应用：结合实际工程问题，应用罗斯兰模型进行热分析和设计。
文献阅读：阅读最新的研究文献，了解该领域的前沿进展。

附录：Python仿真代码

本文所有案例的Python仿真代码详见同目录下的 run_simulation.py 文件。代码包括以下主要功能：

罗斯兰扩散方程求解器：实现了一维、二维光学厚辐射的数值求解。
有效导热系数计算模块：计算不同温度下的辐射导热系数和有效导热系数。
边界条件处理模块：实现跳跃边界条件和 slip 边界条件。
可视化工具：生成温度分布、热流分布、导热系数曲线等图表。
GIF动画生成：创建辐射扩散过程的动态可视化。

代码采用模块化设计，便于扩展和修改。读者可以根据具体需求，修改介质参数、几何形状和边界条件，进行定制化的仿真分析。

"""
主题020：光学厚介质辐射换热 - Python仿真程序

本程序实现了五个光学厚介质辐射换热的典型案例：
1. 光学厚平板辐射换热
2. 球形容器辐射传热
3. 高温炉膛辐射导热
4. 辐射与导热耦合传热
5. 非均匀光学厚介质辐射

作者：仿真教学团队
日期：2026-03-01
"""

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.colors import LinearSegmentedColormap
import matplotlib.animation as animation
from scipy import integrate
from scipy.integrate import trapezoid
import warnings
warnings.filterwarnings('ignore')

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

# 物理常数
sigma = 5.670374419e-8  # 斯特藩-玻尔兹曼常数, W/(m²·K⁴)

print("=" * 80)
print("主题020：光学厚介质辐射换热仿真")
print("=" * 80)

# ==============================================================================
# 案例一：光学厚平板辐射换热
# ==============================================================================
def case1_optically_thick_plates():
    """
    案例一：光学厚平板辐射换热
    
    计算两个平行平板间光学厚气体的辐射换热，比较罗斯兰近似与精确解
    """
    print("\n[案例一] 光学厚平板辐射换热")
    print("-" * 60)
    
    # 物理参数
    T1 = 2000.0  # 上平板温度, K
    T2 = 1500.0  # 下平板温度, K
    L = 0.5      # 平板间距, m
    kappa = 10.0 # 气体吸收系数, m^-1
    
    # 光学厚度
    tau = kappa * L
    print(f"光学厚度 τ = {tau:.2f}")
    print(f"光学厚条件: {'满足' if tau > 10 else '不满足'}")
    
    # 空间离散
    Nz = 100
    z = np.linspace(0, L, Nz)
    
    # 罗斯兰近似解（常导热系数）
    T_avg = (T1 + T2) / 2
    k_r_const = 16 * sigma * T_avg**3 / (3 * kappa)  # 常辐射导热系数
    T_rosseland_const = T1 + (T2 - T1) * z / L
    q_rosseland_const = k_r_const * (T1 - T2) / L
    
    # 考虑温度变化的精确解
    # T^4 = T1^4 + (T2^4 - T1^4) * z/L
    T_exact = (T1**4 + (T2**4 - T1**4) * z / L)**0.25
    q_exact = 4 * sigma * (T1**4 - T2**4) / (3 * kappa * L)
    
    # 辐射导热系数分布
    k_r_exact = 16 * sigma * T_exact**3 / (3 * kappa)
    
    print(f"\n辐射热流:")
    print(f"  罗斯兰近似(常k): {q_rosseland_const:.2f} W/m²")
    print(f"  精确解: {q_exact:.2f} W/m²")
    print(f"  相对误差: {abs(q_rosseland_const - q_exact)/q_exact * 100:.2f}%")
    
    # 可视化
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 温度分布
    ax1 = axes[0, 0]
    ax1.plot(z, T_rosseland_const, 'b-', linewidth=2, label='Rosseland (constant $k_r$)')
    ax1.plot(z, T_exact, 'r--', linewidth=2, label='Exact solution')
    ax1.set_xlabel('Position z (m)', fontsize=11)
    ax1.set_ylabel('Temperature (K)', fontsize=11)
    ax1.set_title('Temperature Distribution', fontsize=12)
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # 辐射导热系数分布
    ax2 = axes[0, 1]
    ax2.plot(z, k_r_exact, 'g-', linewidth=2)
    ax2.axhline(y=k_r_const, color='b', linestyle='--', label=f'Average $k_r$ = {k_r_const:.2f}')
    ax2.set_xlabel('Position z (m)', fontsize=11)
    ax2.set_ylabel('Radiative Conductivity (W/(m·K))', fontsize=11)
    ax2.set_title('Radiative Conductivity Distribution', fontsize=12)
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    # T^4分布（显示线性关系）
    ax3 = axes[1, 0]
    ax3.plot(z, T_exact**4 / 1e12, 'r-', linewidth=2)
    ax3.set_xlabel('Position z (m)', fontsize=11)
    ax3.set_ylabel('T⁴ (×10¹² K⁴)', fontsize=11)
    ax3.set_title('T⁴ Distribution (Linear in Exact Solution)', fontsize=12)
    ax3.grid(True, alpha=0.3)
    
    # 热流对比
    ax4 = axes[1, 1]
    methods = ['Rosseland\n(constant k)', 'Exact']
    heat_fluxes = [q_rosseland_const, q_exact]
    colors = ['blue', 'red']
    bars = ax4.bar(methods, heat_fluxes, color=colors, alpha=0.7)
    ax4.set_ylabel('Heat Flux (W/m²)', fontsize=11)
    ax4.set_title('Heat Flux Comparison', fontsize=12)
    ax4.grid(True, alpha=0.3, axis='y')
    
    # 添加数值标签
    for bar, q in zip(bars, heat_fluxes):
        height = bar.get_height()
        ax4.text(bar.get_x() + bar.get_width()/2., height,
                f'{q:.1f}',
                ha='center', va='bottom', fontsize=10)
    
    plt.tight_layout()
    plt.savefig('case1_optically_thick_plates.png', dpi=150, bbox_inches='tight')
    print("  图片已保存: case1_optically_thick_plates.png")
    plt.close()
    
    return z, T_exact, q_exact, k_r_exact

# ==============================================================================
# 案例二：球形容器辐射传热
# ==============================================================================
def case2_spherical_container():
    """
    案例二：球形容器辐射传热
    
    计算球形容器内光学厚气体的温度分布和辐射热流
    """
    print("\n[案例二] 球形容器辐射传热")
    print("-" * 60)
    
    # 物理参数
    R1 = 0.2  # 内半径, m
    R2 = 0.5  # 外半径, m
    T1 = 2200.0  # 内壁温度, K
    T2 = 1800.0  # 外壁温度, K
    kappa = 15.0  # 气体吸收系数, m^-1
    
    # 光学厚度
    tau = kappa * (R2 - R1)
    print(f"径向光学厚度 τ = {tau:.2f}")
    print(f"光学厚条件: {'满足' if tau > 10 else '不满足'}")
    
    # 径向离散
    Nr = 100
    r = np.linspace(R1, R2, Nr)
    
    # 罗斯兰近似解（常导热系数）
    T_avg = (T1 + T2) / 2
    k_r = 16 * sigma * T_avg**3 / (3 * kappa)
    
    # 温度分布: T = T1 + (T2 - T1) * (1/r - 1/R1) / (1/R2 - 1/R1)
    T_rosseland = T1 + (T2 - T1) * (1/r - 1/R1) / (1/R2 - 1/R1)
    
    # 辐射热流密度
    q_r = k_r * (T1 - T2) / (r**2 * (1/R1 - 1/R2))
    
    # 总热流量
    Q_total = 4 * np.pi * k_r * (T1 - T2) / (1/R1 - 1/R2)
    
    # 考虑温度变化的精确解
    # T^4 = T1^4 + (T2^4 - T1^4) * (1/r - 1/R1) / (1/R2 - 1/R1)
    T_exact = (T1**4 + (T2**4 - T1**4) * (1/r - 1/R1) / (1/R2 - 1/R1))**0.25
    
    # 精确热流
    q_exact = 4 * sigma * (T1**4 - T2**4) / (3 * kappa * r**2 * (1/R1 - 1/R2))
    Q_exact = 4 * np.pi * sigma * (T1**4 - T2**4) / (3 * kappa * (1/R1 - 1/R2))
    
    print(f"\n温度分布:")
    print(f"  内壁(r={R1}m): T = {T_exact[0]:.1f} K")
    print(f"  外壁(r={R2}m): T = {T_exact[-1]:.1f} K")
    print(f"\n热流量:")
    print(f"  罗斯兰近似: {Q_total:.2f} W")
    print(f"  精确解: {Q_exact:.2f} W")
    print(f"  相对误差: {abs(Q_total - Q_exact)/Q_exact * 100:.2f}%")
    
    # 可视化
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 温度分布
    ax1 = axes[0, 0]
    ax1.plot(r, T_rosseland, 'b-', linewidth=2, label='Rosseland')
    ax1.plot(r, T_exact, 'r--', linewidth=2, label='Exact')
    ax1.set_xlabel('Radius r (m)', fontsize=11)
    ax1.set_ylabel('Temperature (K)', fontsize=11)
    ax1.set_title('Temperature Distribution in Spherical Container', fontsize=12)
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # 热流密度分布
    ax2 = axes[0, 1]
    ax2.plot(r, q_r, 'b-', linewidth=2, label='Rosseland')
    ax2.plot(r, q_exact, 'r--', linewidth=2, label='Exact')
    ax2.set_xlabel('Radius r (m)', fontsize=11)
    ax2.set_ylabel('Heat Flux (W/m²)', fontsize=11)
    ax2.set_title('Radial Heat Flux Distribution', fontsize=12)
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    # T^4分布
    ax3 = axes[1, 0]
    ax3.plot(r, T_exact**4 / 1e13, 'g-', linewidth=2)
    ax3.set_xlabel('Radius r (m)', fontsize=11)
    ax3.set_ylabel('T⁴ (×10¹³ K⁴)', fontsize=11)
    ax3.set_title('T⁴ Distribution', fontsize=12)
    ax3.grid(True, alpha=0.3)
    
    # 热流量对比
    ax4 = axes[1, 1]
    methods = ['Rosseland', 'Exact']
    heat_flows = [Q_total, Q_exact]
    colors = ['blue', 'red']
    bars = ax4.bar(methods, heat_flows, color=colors, alpha=0.7)
    ax4.set_ylabel('Total Heat Flow (W)', fontsize=11)
    ax4.set_title('Total Heat Flow Comparison', fontsize=12)
    ax4.grid(True, alpha=0.3, axis='y')
    
    for bar, Q in zip(bars, heat_flows):
        height = bar.get_height()
        ax4.text(bar.get_x() + bar.get_width()/2., height,
                f'{Q:.1f}',
                ha='center', va='bottom', fontsize=10)
    
    plt.tight_layout()
    plt.savefig('case2_spherical_container.png', dpi=150, bbox_inches='tight')
    print("  图片已保存: case2_spherical_container.png")
    plt.close()
    
    return r, T_exact, q_exact, Q_exact

# ==============================================================================
# 案例三：高温炉膛辐射导热
# ==============================================================================
def case3_furnace_radiation():
    """
    案例三：高温炉膛辐射导热
    
    计算三维炉膛内的温度分布和辐射热流
    """
    print("\n[案例三] 高温炉膛辐射导热")
    print("-" * 60)
    
    # 炉膛参数
    Lx, Ly, Lz = 2.0, 2.0, 4.0  # 炉膛尺寸, m
    T_wall = 1500.0  # 壁面温度, K
    T_center = 2000.0  # 中心温度, K
    kappa = 8.0  # 气体吸收系数, m^-1
    
    # 光学厚度
    tau_x = kappa * Lx
    tau_y = kappa * Ly
    tau_z = kappa * Lz
    print(f"光学厚度: τ_x = {tau_x:.1f}, τ_y = {tau_y:.1f}, τ_z = {tau_z:.1f}")
    print(f"光学厚条件: {'满足' if min(tau_x, tau_y, tau_z) > 10 else '部分满足'}")
    
    # 网格离散
    Nx, Ny, Nz = 30, 30, 40
    x = np.linspace(-Lx/2, Lx/2, Nx)
    y = np.linspace(-Ly/2, Ly/2, Ny)
    z = np.linspace(-Lz/2, Lz/2, Nz)
    
    dx, dy, dz = x[1]-x[0], y[1]-y[0], z[1]-z[0]
    
    # 初始化温度场
    T = np.zeros((Nx, Ny, Nz))
    
    # 使用近似解析解：T^4 = T_center^4 - (T_center^4 - T_wall^4) * (r/r_max)^2
    # 其中 r 是到中心的归一化距离
    for i in range(Nx):
        for j in range(Ny):
            for k in range(Nz):
                # 归一化距离
                r_norm = np.sqrt((x[i]/(Lx/2))**2 + (y[j]/(Ly/2))**2 + (z[k]/(Lz/2))**2)
                r_norm = min(r_norm, 1.0)
                # 温度分布
                T4 = T_center**4 - (T_center**4 - T_wall**4) * r_norm**2
                T[i, j, k] = T4**0.25
    
    # 计算辐射导热系数
    k_r = 16 * sigma * T**3 / (3 * kappa)
    
    # 计算热流（中心截面）
    k = Nz // 2
    T_slice = T[:, :, k]
    k_r_slice = k_r[:, :, k]
    
    # 热流计算（使用中心差分）
    q_x = np.zeros((Nx, Ny))
    q_y = np.zeros((Nx, Ny))
    
    for i in range(1, Nx-1):
        for j in range(1, Ny-1):
            k_avg_x = (k_r_slice[i+1, j] + k_r_slice[i-1, j]) / 2
            k_avg_y = (k_r_slice[i, j+1] + k_r_slice[i, j-1]) / 2
            q_x[i, j] = -k_avg_x * (T_slice[i+1, j] - T_slice[i-1, j]) / (2*dx)
            q_y[i, j] = -k_avg_y * (T_slice[i, j+1] - T_slice[i, j-1]) / (2*dy)
    
    q_magnitude = np.sqrt(q_x**2 + q_y**2)
    
    print(f"\n温度场统计:")
    print(f"  中心温度: {T[Nx//2, Ny//2, Nz//2]:.1f} K")
    print(f"  壁面温度: {T[0, Ny//2, Nz//2]:.1f} K")
    print(f"  最大热流: {np.max(q_magnitude):.2f} W/m²")
    
    # 可视化
    fig = plt.figure(figsize=(16, 12))
    
    # 温度分布3D
    ax1 = fig.add_subplot(2, 3, 1, projection='3d')
    X, Y = np.meshgrid(x, y)
    surf1 = ax1.plot_surface(X, Y, T_slice.T, cmap='hot', alpha=0.8)
    ax1.set_xlabel('x (m)')
    ax1.set_ylabel('y (m)')
    ax1.set_zlabel('T (K)')
    ax1.set_title('Temperature Distribution (z=0)')
    fig.colorbar(surf1, ax=ax1, shrink=0.5)
    
    # 温度等高线
    ax2 = fig.add_subplot(2, 3, 2)
    contour = ax2.contourf(X, Y, T_slice.T, levels=20, cmap='hot')
    ax2.set_xlabel('x (m)')
    ax2.set_ylabel('y (m)')
    ax2.set_title('Temperature Contours')
    ax2.set_aspect('equal')
    fig.colorbar(contour, ax=ax2)
    
    # 辐射导热系数
    ax3 = fig.add_subplot(2, 3, 3)
    k_r_center = k_r[Nx//2, Ny//2, :]
    ax3.plot(z, k_r_center, 'r-', linewidth=2)
    ax3.set_xlabel('z (m)')
    ax3.set_ylabel('$k_r$ (W/(m·K))')
    ax3.set_title('Radiative Conductivity (center line)')
    ax3.grid(True, alpha=0.3)
    
    # 热流矢量图
    ax4 = fig.add_subplot(2, 3, 4)
    skip = 3
    ax4.quiver(X[::skip, ::skip], Y[::skip, ::skip], 
               q_x[::skip, ::skip].T, q_y[::skip, ::skip].T,
               q_magnitude[::skip, ::skip].T, cmap='jet', scale=1000)
    ax4.set_xlabel('x (m)')
    ax4.set_ylabel('y (m)')
    ax4.set_title('Heat Flux Vectors')
    ax4.set_aspect('equal')
    
    # 热流大小分布
    ax5 = fig.add_subplot(2, 3, 5)
    im = ax5.imshow(q_magnitude.T, extent=[-Lx/2, Lx/2, -Ly/2, Ly/2], 
                    origin='lower', cmap='jet')
    ax5.set_xlabel('x (m)')
    ax5.set_ylabel('y (m)')
    ax5.set_title('Heat Flux Magnitude')
    fig.colorbar(im, ax=ax5)
    
    # 沿x轴的温度分布
    ax6 = fig.add_subplot(2, 3, 6)
    ax6.plot(x, T_slice[:, Ny//2], 'b-', linewidth=2, label='T(x, y=0)')
    ax6.plot(x, T_slice[0, :], 'r--', linewidth=2, label='T(x, y=±1)')
    ax6.set_xlabel('x (m)')
    ax6.set_ylabel('Temperature (K)')
    ax6.set_title('Temperature Profiles')
    ax6.legend()
    ax6.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('case3_furnace_radiation.png', dpi=150, bbox_inches='tight')
    print("  图片已保存: case3_furnace_radiation.png")
    plt.close()
    
    return T, k_r, q_magnitude

# ==============================================================================
# 案例四：辐射与导热耦合传热
# ==============================================================================
def case4_coupled_conduction_radiation():
    """
    案例四：辐射与导热耦合传热
    
    计算复合墙体中导热和辐射的耦合传热
    """
    print("\n[案例四] 辐射与导热耦合传热")
    print("-" * 60)
    
    # 墙体参数
    L1 = 0.3  # 耐火砖厚度, m
    L2 = 0.2  # 隔热材料厚度, m
    L_total = L1 + L2
    
    k1 = 1.5  # 耐火砖导热系数, W/(m·K)
    k2 = 0.1  # 隔热材料导热系数, W/(m·K)
    
    T_i = 1800.0  # 内表面温度, K
    T_o = 300.0   # 外表面温度, K
    
    kappa = 5.0  # 气体吸收系数, m^-1 (仅在内层高温区有效)
    
    # 空间离散
    Nz = 200
    z = np.linspace(0, L_total, Nz)
    dz = z[1] - z[0]
    
    # 材料属性数组
    k_c = np.zeros(Nz)  # 导热系数
    kappa_arr = np.zeros(Nz)  # 吸收系数
    
    for i in range(Nz):
        if z[i] <= L1:
            k_c[i] = k1
            kappa_arr[i] = kappa
        else:
            k_c[i] = k2
            kappa_arr[i] = 0  # 隔热层无辐射
    
    # 迭代求解温度分布
    T = np.linspace(T_i, T_o, Nz)  # 初始猜测
    max_iter = 1000
    tol = 1e-6
    
    for iteration in range(max_iter):
        T_old = T.copy()
        
        # 计算辐射导热系数
        k_r = np.zeros(Nz)
        for i in range(Nz):
            if kappa_arr[i] > 0:
                k_r[i] = 16 * sigma * T[i]**3 / (3 * kappa_arr[i])
        
        # 有效导热系数
        k_eff = k_c + k_r
        
        # 有限差分求解
        for i in range(1, Nz-1):
            k_plus = (k_eff[i+1] + k_eff[i]) / 2
            k_minus = (k_eff[i] + k_eff[i-1]) / 2
            
            a_p = k_plus + k_minus
            a_e = k_plus
            a_w = k_minus
            
            T[i] = (a_e * T[i+1] + a_w * T[i-1]) / a_p
        
        # 边界条件
        T[0] = T_i
        T[-1] = T_o
        
        # 检查收敛
        error = np.max(np.abs(T - T_old))
        if error < tol:
            print(f"迭代收敛于第 {iteration+1} 次迭代")
            break
    
    # 计算热流
    q_cond = k_c * (T[0] - T[-1]) / L_total  # 仅导热热流
    q_total = np.zeros(Nz)
    for i in range(Nz-1):
        k_eff_avg = (k_eff[i] + k_eff[i+1]) / 2
        q_total[i] = k_eff_avg * (T[i] - T[i+1]) / dz
    
    q_avg = np.mean(q_total[:-1])
    
    # 计算各层热阻
    R_cond_only = L1/k1 + L2/k2
    q_cond_only = (T_i - T_o) / R_cond_only
    
    print(f"\n温度分布:")
    print(f"  界面温度: {T[int(Nz*L1/L_total)]:.1f} K")
    print(f"\n热流:")
    print(f"  仅导热: {q_cond_only:.2f} W/m²")
    print(f"  耦合传热: {q_avg:.2f} W/m²")
    print(f"  辐射增强: {(q_avg/q_cond_only - 1)*100:.1f}%")
    
    # 可视化
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 温度分布
    ax1 = axes[0, 0]
    ax1.plot(z, T, 'b-', linewidth=2, label='Coupled')
    
    # 仅导热情况
    T_cond = np.zeros(Nz)
    for i in range(Nz):
        if z[i] <= L1:
            T_cond[i] = T_i - (T_i - T_o) * (z[i]/k1) / (L1/k1 + L2/k2)
        else:
            T_cond[i] = T_i - (T_i - T_o) * (L1/k1 + (z[i]-L1)/k2) / (L1/k1 + L2/k2)
    ax1.plot(z, T_cond, 'r--', linewidth=2, label='Conduction only')
    ax1.axvline(x=L1, color='k', linestyle=':', alpha=0.5, label='Interface')
    ax1.set_xlabel('Position z (m)', fontsize=11)
    ax1.set_ylabel('Temperature (K)', fontsize=11)
    ax1.set_title('Temperature Distribution', fontsize=12)
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # 有效导热系数
    ax2 = axes[0, 1]
    ax2.semilogy(z, k_eff, 'g-', linewidth=2, label='$k_{eff}$')
    ax2.plot(z, k_c, 'b--', linewidth=2, label='$k_c$')
    ax2.plot(z, k_r, 'r:', linewidth=2, label='$k_r$')
    ax2.axvline(x=L1, color='k', linestyle=':', alpha=0.5)
    ax2.set_xlabel('Position z (m)', fontsize=11)
    ax2.set_ylabel('Conductivity (W/(m·K))', fontsize=11)
    ax2.set_title('Conductivity Distribution', fontsize=12)
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    # 热流分布
    ax3 = axes[1, 0]
    ax3.plot(z[:-1], q_total[:-1], 'm-', linewidth=2)
    ax3.axhline(y=q_cond_only, color='r', linestyle='--', label=f'Conduction only: {q_cond_only:.1f}')
    ax3.axvline(x=L1, color='k', linestyle=':', alpha=0.5)
    ax3.set_xlabel('Position z (m)', fontsize=11)
    ax3.set_ylabel('Heat Flux (W/m²)', fontsize=11)
    ax3.set_title('Heat Flux Distribution', fontsize=12)
    ax3.legend()
    ax3.grid(True, alpha=0.3)
    
    # 辐射导热系数vs温度
    ax4 = axes[1, 1]
    T_range = np.linspace(300, 2000, 100)
    k_r_range = 16 * sigma * T_range**3 / (3 * kappa)
    ax4.semilogy(T_range, k_r_range, 'r-', linewidth=2)
    ax4.axhline(y=k1, color='b', linestyle='--', label=f'$k_1$ = {k1}')
    ax4.axhline(y=k2, color='g', linestyle='--', label=f'$k_2$ = {k2}')
    ax4.set_xlabel('Temperature (K)', fontsize=11)
    ax4.set_ylabel('Radiative Conductivity (W/(m·K))', fontsize=11)
    ax4.set_title('$k_r$ vs Temperature', fontsize=12)
    ax4.legend()
    ax4.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('case4_coupled_conduction_radiation.png', dpi=150, bbox_inches='tight')
    print("  图片已保存: case4_coupled_conduction_radiation.png")
    plt.close()
    
    return z, T, k_eff, q_avg

# ==============================================================================
# 案例五：非均匀光学厚介质辐射
# ==============================================================================
def case5_nonuniform_medium():
    """
    案例五：非均匀光学厚介质辐射
    
    计算吸收系数随温度变化的光学厚介质辐射
    """
    print("\n[案例五] 非均匀光学厚介质辐射")
    print("-" * 60)
    
    # 物理参数
    L = 0.4  # 介质层厚度, m
    T1 = 2000.0  # 上表面温度, K
    T2 = 1400.0  # 下表面温度, K
    
    kappa_0 = 10.0  # 参考吸收系数, m^-1
    T_0 = 1500.0    # 参考温度, K
    n = -1.5        # 温度指数（吸收系数随温度升高而降低）
    
    # 空间离散
    Nz = 100
    z = np.linspace(0, L, Nz)
    dz = z[1] - z[0]
    
    # 迭代求解
    T = np.linspace(T1, T2, Nz)  # 初始猜测
    max_iter = 1000
    tol = 1e-8
    
    for iteration in range(max_iter):
        T_old = T.copy()
        
        # 计算吸收系数（随温度变化）
        kappa = kappa_0 * (T / T_0)**n
        
        # 计算辐射导热系数
        k_r = 16 * sigma * T**3 / (3 * kappa)
        
        # 有限差分求解
        for i in range(1, Nz-1):
            k_plus = (k_r[i+1] + k_r[i]) / 2
            k_minus = (k_r[i] + k_r[i-1]) / 2
            
            a_p = k_plus + k_minus
            a_e = k_plus
            a_w = k_minus
            
            T[i] = (a_e * T[i+1] + a_w * T[i-1]) / a_p
        
        # 边界条件
        T[0] = T1
        T[-1] = T2
        
        # 检查收敛
        error = np.max(np.abs(T - T_old))
        if error < tol:
            print(f"迭代收敛于第 {iteration+1} 次迭代")
            break
    
    # 计算最终参数
    kappa_final = kappa_0 * (T / T_0)**n
    k_r_final = 16 * sigma * T**3 / (3 * kappa_final)
    
    # 热流
    q = np.zeros(Nz)
    for i in range(Nz-1):
        k_avg = (k_r_final[i] + k_r_final[i+1]) / 2
        q[i] = k_avg * (T[i] - T[i+1]) / dz
    q_avg = np.mean(q[:-1])
    
    # 与常吸收系数情况对比
    kappa_const = kappa_0
    k_r_const = 16 * sigma * T**3 / (3 * kappa_const)
    
    # 常吸收系数的解析解: T^4 = T1^4 + (T2^4 - T1^4) * z/L
    T_const = (T1**4 + (T2**4 - T1**4) * z / L)**0.25
    q_const = 4 * sigma * (T1**4 - T2**4) / (3 * kappa_const * L)
    
    print(f"\n温度分布:")
    print(f"  上表面: {T[0]:.1f} K")
    print(f"  下表面: {T[-1]:.1f} K")
    print(f"  中间位置: {T[Nz//2]:.1f} K")
    print(f"\n吸收系数:")
    print(f"  上表面: {kappa_final[0]:.2f} m⁻¹")
    print(f"  下表面: {kappa_final[-1]:.2f} m⁻¹")
    print(f"\n热流:")
    print(f"  变吸收系数: {q_avg:.2f} W/m²")
    print(f"  常吸收系数: {q_const:.2f} W/m²")
    print(f"  差异: {abs(q_avg - q_const)/q_const * 100:.2f}%")
    
    # 可视化
    fig, axes = plt.subplots(2, 3, figsize=(16, 10))
    
    # 温度分布对比
    ax1 = axes[0, 0]
    ax1.plot(z, T, 'b-', linewidth=2, label='Variable κ')
    ax1.plot(z, T_const, 'r--', linewidth=2, label='Constant κ')
    ax1.set_xlabel('Position z (m)', fontsize=11)
    ax1.set_ylabel('Temperature (K)', fontsize=11)
    ax1.set_title('Temperature Distribution', fontsize=12)
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # 吸收系数分布
    ax2 = axes[0, 1]
    ax2.plot(z, kappa_final, 'g-', linewidth=2)
    ax2.axhline(y=kappa_const, color='r', linestyle='--', label=f'$\kappa_0$ = {kappa_const}')
    ax2.set_xlabel('Position z (m)', fontsize=11)
    ax2.set_ylabel('Absorption Coefficient (m⁻¹)', fontsize=11)
    ax2.set_title('Absorption Coefficient Distribution', fontsize=12)
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    # 辐射导热系数
    ax3 = axes[0, 2]
    ax3.semilogy(z, k_r_final, 'm-', linewidth=2, label='Variable κ')
    ax3.semilogy(z, k_r_const, 'c--', linewidth=2, label='Constant κ')
    ax3.set_xlabel('Position z (m)', fontsize=11)
    ax3.set_ylabel('Radiative Conductivity (W/(m·K))', fontsize=11)
    ax3.set_title('Radiative Conductivity', fontsize=12)
    ax3.legend()
    ax3.grid(True, alpha=0.3)
    
    # T^4分布
    ax4 = axes[1, 0]
    ax4.plot(z, T**4 / 1e12, 'b-', linewidth=2, label='Variable κ')
    ax4.plot(z, T_const**4 / 1e12, 'r--', linewidth=2, label='Constant κ')
    ax4.set_xlabel('Position z (m)', fontsize=11)
    ax4.set_ylabel('T⁴ (×10¹² K⁴)', fontsize=11)
    ax4.set_title('T⁴ Distribution', fontsize=12)
    ax4.legend()
    ax4.grid(True, alpha=0.3)
    
    # 热流分布
    ax5 = axes[1, 1]
    ax5.plot(z[:-1], q[:-1], 'b-', linewidth=2, label='Variable κ')
    ax5.axhline(y=q_const, color='r', linestyle='--', label=f'Constant κ: {q_const:.1f}')
    ax5.set_xlabel('Position z (m)', fontsize=11)
    ax5.set_ylabel('Heat Flux (W/m²)', fontsize=11)
    ax5.set_title('Heat Flux Distribution', fontsize=12)
    ax5.legend()
    ax5.grid(True, alpha=0.3)
    
    # 吸收系数vs温度
    ax6 = axes[1, 2]
    T_range = np.linspace(1000, 2500, 100)
    kappa_range = kappa_0 * (T_range / T_0)**n
    ax6.plot(T_range, kappa_range, 'g-', linewidth=2)
    ax6.axvline(x=T_0, color='k', linestyle=':', alpha=0.5, label=f'$T_0$ = {T_0} K')
    ax6.axhline(y=kappa_0, color='r', linestyle=':', alpha=0.5, label=f'$\kappa_0$ = {kappa_0}')
    ax6.set_xlabel('Temperature (K)', fontsize=11)
    ax6.set_ylabel('Absorption Coefficient (m⁻¹)', fontsize=11)
    ax6.set_title('$\kappa$ vs T ($\kappa \propto T^{-1.5}$)', fontsize=12)
    ax6.legend()
    ax6.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('case5_nonuniform_medium.png', dpi=150, bbox_inches='tight')
    print("  图片已保存: case5_nonuniform_medium.png")
    plt.close()
    
    return z, T, kappa_final, k_r_final

# ==============================================================================
# 生成GIF动画
# ==============================================================================
def create_animations():
    """生成GIF动画展示辐射扩散过程"""
    print("\n" + "=" * 80)
    print("生成GIF动画")
    print("=" * 80)
    
    # 动画1：辐射扩散过程
    print("\n[动画1] 辐射扩散过程")
    fig, ax = plt.subplots(figsize=(10, 6))
    
    L = 0.5
    kappa = 10.0
    Nz = 50
    z = np.linspace(0, L, Nz)
    
    # 初始温度分布（阶跃）
    T1, T2 = 2000.0, 1500.0
    T_initial = np.where(z < L/2, T1, T2)
    
    # 时间演化
    dt = 0.001
    Nt = 100
    
    line, = ax.plot(z, T_initial, 'b-', linewidth=2)
    ax.set_xlabel('Position z (m)', fontsize=12)
    ax.set_ylabel('Temperature (K)', fontsize=12)
    ax.set_title('Radiation Diffusion Process', fontsize=14)
    ax.set_ylim(1400, 2100)
    ax.grid(True, alpha=0.3)
    
    T = T_initial.copy()
    
    def animate(frame):
        nonlocal T
        for _ in range(5):  # 每帧迭代5次
            k_r = 16 * sigma * T**3 / (3 * kappa)
            T_new = T.copy()
            for i in range(1, Nz-1):
                k_avg = (k_r[i+1] + 2*k_r[i] + k_r[i-1]) / 4
                T_new[i] = T[i] + dt * k_avg * (T[i+1] - 2*T[i] + T[i-1]) / ((z[1]-z[0])**2)
            T_new[0] = T1
            T_new[-1] = T2
            T = T_new
        
        line.set_ydata(T)
        ax.set_title(f'Radiation Diffusion (t = {frame*dt*5:.4f} s)', fontsize=14)
        return line,
    
    anim = animation.FuncAnimation(fig, animate, frames=Nt, interval=100, blit=True)
    anim.save('animation1_radiation_diffusion.gif', writer='pillow', fps=10)
    print("  动画已保存: animation1_radiation_diffusion.gif")
    plt.close()
    
    # 动画2：有效导热系数随温度变化
    print("\n[动画2] 有效导热系数随温度变化")
    fig, axes = plt.subplots(1, 2, figsize=(14, 6))
    
    T_range = np.linspace(300, 2500, 100)
    kappa_values = [5, 10, 20]
    colors = ['blue', 'green', 'red']
    
    lines = []
    for ax in axes:
        for kappa, color in zip(kappa_values, colors):
            line, = ax.plot([], [], color=color, linewidth=2, label=f'κ = {kappa}')
            lines.append(line)
        ax.set_xlabel('Temperature (K)', fontsize=11)
        ax.set_ylabel('Conductivity (W/(m·K))', fontsize=11)
        ax.legend()
        ax.grid(True, alpha=0.3)
    
    axes[0].set_title('Radiative Conductivity $k_r$', fontsize=12)
    axes[0].set_yscale('log')
    axes[1].set_title('Effective Conductivity $k_{eff}$ (with $k_c$ = 1)', fontsize=12)
    axes[1].set_yscale('log')
    
    k_c = 1.0
    
    def animate2(frame):
        n = frame / 20  # 温度指数变化
        for i, (kappa, color) in enumerate(zip(kappa_values, colors)):
            k_r = 16 * sigma * T_range**3 / (3 * kappa)
            k_eff = k_c + k_r
            lines[i].set_data(T_range, k_r)
            lines[i+3].set_data(T_range, k_eff)
        
        for ax in axes:
            ax.set_xlim(300, 2500)
            ax.set_ylim(1e-2, 1e4)
        
        return lines
    
    anim2 = animation.FuncAnimation(fig, animate2, frames=50, interval=100, blit=True)
    anim2.save('animation2_conductivity_variation.gif', writer='pillow', fps=10)
    print("  动画已保存: animation2_conductivity_variation.gif")
    plt.close()
    
    # 动画3：光学厚度对温度分布的影响
    print("\n[动画3] 光学厚度对温度分布的影响")
    fig, ax = plt.subplots(figsize=(10, 6))
    
    L = 0.5
    T1, T2 = 2000.0, 1500.0
    Nz = 50
    z = np.linspace(0, L, Nz)
    
    line, = ax.plot([], [], 'b-', linewidth=2)
    ax.set_xlabel('Position z (m)', fontsize=12)
    ax.set_ylabel('Temperature (K)', fontsize=12)
    ax.set_title('Effect of Optical Thickness on Temperature Profile', fontsize=14)
    ax.set_xlim(0, L)
    ax.set_ylim(1400, 2100)
    ax.grid(True, alpha=0.3)
    
    def animate3(frame):
        kappa = 2 + frame * 0.5  # 吸收系数从2增加到27
        tau = kappa * L
        
        # 温度分布: T^4 = T1^4 + (T2^4 - T1^4) * z/L
        T = (T1**4 + (T2**4 - T1**4) * z / L)**0.25
        
        line.set_data(z, T)
        ax.set_title(f'Optical Thickness τ = {tau:.1f}, κ = {kappa:.1f} m⁻¹', fontsize=14)
        return line,
    
    anim3 = animation.FuncAnimation(fig, animate3, frames=50, interval=100, blit=True)
    anim3.save('animation3_optical_thickness_effect.gif', writer='pillow', fps=10)
    print("  动画已保存: animation3_optical_thickness_effect.gif")
    plt.close()

# ==============================================================================
# 主程序
# ==============================================================================
if __name__ == "__main__":
    # 运行所有案例
    print("\n运行所有仿真案例...")
    print("=" * 80)
    
    # 案例一：光学厚平板辐射换热
    z1, T1, q1, kr1 = case1_optically_thick_plates()
    
    # 案例二：球形容器辐射传热
    r2, T2, q2, Q2 = case2_spherical_container()
    
    # 案例三：高温炉膛辐射导热
    T3, kr3, q3 = case3_furnace_radiation()
    
    # 案例四：辐射与导热耦合传热
    z4, T4, keff4, q4 = case4_coupled_conduction_radiation()
    
    # 案例五：非均匀光学厚介质辐射
    z5, T5, kappa5, kr5 = case5_nonuniform_medium()
    
    # 生成GIF动画
    create_animations()
    
    print("\n" + "=" * 80)
    print("所有仿真案例运行完成!")
    print("=" * 80)
    print("\n生成的文件:")
    print("  - case1_optically_thick_plates.png")
    print("  - case2_spherical_container.png")
    print("  - case3_furnace_radiation.png")
    print("  - case4_coupled_conduction_radiation.png")
    print("  - case5_nonuniform_medium.png")
    print("  - animation1_radiation_diffusion.gif")
    print("  - animation2_conductivity_variation.gif")
    print("  - animation3_optical_thickness_effect.gif")
    print("\n" + "=" * 80)