主题019：光学薄介质辐射换热

摘要

光学薄介质是辐射换热中一类重要的特殊情况，当介质的光学厚度远小于1时，辐射传递过程呈现出独特的物理特性。本教程系统研究光学薄介质的辐射特性，建立光学薄极限下的简化模型，深入讨论发射率与光学厚度的关系。通过五个典型应用案例——等温平板间光学薄气体辐射、非等温介质辐射换热、光学薄极限下的发射率计算、多层光学薄介质辐射以及光学薄介质与壁面辐射耦合问题——全面展示光学薄介质辐射换热的理论方法和工程应用。所有案例均配有完整的Python仿真程序和可视化分析，为工程实践提供可靠的理论指导。

关键词： 光学薄、发射率、光学厚度、简化模型、透明介质、辐射传递

---

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

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

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

$$\tau ={\int }_0^L{\kappa }_{\lambda }(s)ds$$

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

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

光学薄条件意味着介质对辐射几乎是"透明"的，辐射可以几乎无衰减地穿过介质。这一特性在航天器热控、高温气体辐射、透明材料热分析等领域具有重要应用价值。

1.2 光学薄极限的数学表述

在光学薄极限下（$\tau \to 0$），辐射传递方程（RTE）可以大幅简化。完整的RTE为：

$$\frac{d{I}_{\lambda }}{ds}=-{\kappa }_{\lambda }{I}_{\lambda }+{\kappa }_{\lambda }{I}_{b\lambda }+\frac{{\sigma }_{s\lambda }}{4\pi }{\int }_{4\pi }{\Phi }_{\lambda }(s^{\prime },s)I_{\lambda }(s^{\prime })d{\Omega }^{\prime }$$

当光学厚度很小时，可以忽略辐射强度的衰减项（$-{\kappa }_{\lambda }{I}_{\lambda }$），方程简化为：

$$\frac{d{I}_{\lambda }}{ds}\approx {\kappa }_{\lambda }{I}_{b\lambda }+\frac{{\sigma }_{s\lambda }}{4\pi }{\int }_{4\pi }{\Phi }_{\lambda }(s^{\prime },s)I_{\lambda }(s^{\prime })d{\Omega }^{\prime }$$

这一简化使得问题的求解难度大幅降低，同时也揭示了光学薄介质辐射的物理本质：辐射主要由介质的发射和散射贡献，吸收效应可以忽略。

1.3 光学薄介质辐射的特征

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

1. 弱吸收特性

辐射穿过光学薄介质时，透射率接近1：

$${\tau }_{trans}=e^{-\tau }\approx 1-\tau +\frac{{\tau }^2}{2}-...\approx 1$$

这意味着入射辐射几乎不被介质吸收，大部分能量能够穿透介质。

2. 发射主导特性

在光学薄极限下，介质的辐射主要由发射贡献。介质发射的辐射强度为：

$$I_{\lambda }^{emission}={\int }_0^L{\kappa }_{\lambda }{I}_{b\lambda }ds={\kappa }_{\lambda }L{I}_{b\lambda }=\tau I_{b\lambda }$$

由于 $\tau \ll 1$，发射强度远小于黑体辐射强度，但发射过程本身不可忽略。

3. 线性叠加特性

光学薄介质中，不同位置的发射贡献可以线性叠加。总辐射强度等于各微元发射贡献的积分：

$$I_{\lambda }=I_{\lambda ,0}{e}^{-\tau }+{\int }_0^{\tau }{I}_{b\lambda }({\tau }^{\prime })e^{-(\tau -{\tau }^{\prime })}d{\tau }^{\prime }$$

在光学薄极限下，指数项近似为1，上式简化为：

$$I_{\lambda }\approx I_{\lambda ,0}+{\int }_0^L{\kappa }_{\lambda }{I}_{b\lambda }ds$$

这表明入射辐射与介质发射可以简单叠加。

4. 局部热力学平衡

光学薄介质通常处于局部热力学平衡（LTE）状态，即介质的温度分布可以用局部平衡温度描述。这一假设大大简化了问题的数学处理。

1.4 光学薄条件的工程应用

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

航天器热控系统

航天器外部包覆的多层隔热材料（MLI）中，层间气体通常处于低压状态，光学厚度很小。光学薄近似可用于分析层间辐射换热。

燃烧产物辐射

燃气轮机、火箭发动机等高温燃烧设备中，燃烧产物（CO₂、H₂O等）在特定波段可能处于光学薄状态。光学薄模型可用于快速估算辐射热损失。

透明介质热分析

玻璃、晶体等透明材料在某些波段具有很低的吸收系数。光学薄近似可用于分析这些材料中的辐射传热。

大气辐射传输

地球大气层在可见光波段近似为光学薄介质。光学薄模型在大气辐射传输、遥感等领域有重要应用。

---

2. 光学薄介质辐射的理论基础

2.1 辐射传递方程的简化

在光学薄极限下，辐射传递方程可以简化为更易于求解的形式。考虑一维平板几何，辐射传递方程为：

$$\mu \frac{dI}{d\tau }+I=(1-\omega )I_b+\frac{\omega }{4\pi }{\int }_{4\pi }\Phi Id\Omega$$

其中，$\mu =\cos \theta$ 是方向余弦，$\omega ={\sigma }_s/(\kappa +{\sigma }_s)$ 是反照率。

在光学薄极限（$\tau \ll 1$）下，假设介质温度均匀，方程的解可以表示为：

$$I(\tau ,\mu )=I(0,\mu )e^{-\tau /\mu }+\frac{1-\omega }{\mu }{\int }_0^{\tau }{I}_b{e}^{-(\tau -{\tau }^{\prime })/\mu }d{\tau }^{\prime }+\frac{\omega }{4\pi \mu }{\int }_0^{\tau }{\int }_{4\pi }\Phi Id\Omega e^{-(\tau -{\tau }^{\prime })/\mu }d{\tau }^{\prime }$$

对于光学薄介质，指数项可以展开为泰勒级数：

$$e^{-\tau /\mu }\approx 1-\frac{\tau }{\mu }+\frac{{\tau }^2}{2{\mu }^2}-...$$

保留一阶项，得到：

$$I(\tau ,\mu )\approx I(0,\mu )+\frac{\tau }{\mu }\left[(1-\omega )I_b-I(0,\mu )+\frac{\omega }{4\pi }{\int }_{4\pi }\Phi I(0,{\mu }^{\prime })d{\Omega }^{\prime }\right]$$

这一表达式表明，在光学薄极限下，辐射强度沿传播方向的变化与光学厚度成正比。

2.2 光学薄介质的发射率

光学薄介质的发射率是描述其辐射特性的重要参数。对于均匀温度的光学薄介质层，光谱发射率定义为：

$${\epsilon }_{\lambda }=1-e^{-{\tau }_{\lambda }}\approx {\tau }_{\lambda }={\kappa }_{\lambda }L$$

总发射率通过对全波谱积分得到：

$$\epsilon =\frac{{\int }_0^{\infty }{\epsilon }_{\lambda }{E}_{b\lambda }d\lambda }{\sigma T^4}=\frac{{\int }_0^{\infty }{\kappa }_{\lambda }L{E}_{b\lambda }d\lambda }{\sigma T^4}$$

对于灰体介质（$\kappa_\lambda = \kappa = $ 常数），发射率简化为：

$$\epsilon =\kappa L=\tau$$

这表明灰体光学薄介质的发射率等于其光学厚度。

2.3 光学薄介质的吸收率

根据基尔霍夫定律，在热力学平衡条件下，介质的发射率等于其吸收率：

$${\epsilon }_{\lambda }={\alpha }_{\lambda }$$

因此，光学薄介质的吸收率也可以表示为：

$${\alpha }_{\lambda }=1-e^{-{\tau }_{\lambda }}\approx {\tau }_{\lambda }$$

这一结果与物理直观一致：光学薄介质对入射辐射的吸收很弱，吸收率近似等于光学厚度。

2.4 光学薄极限下的辐射热流

辐射热流可以通过对辐射强度在所有方向上积分得到：

$$q={\int }_{4\pi }I\cos \theta d\Omega =2\pi {\int }_{-1}^{1}I\mu d\mu$$

对于光学薄介质，辐射热流可以分解为两部分：

$$q=q_{incident}+q_{emission}$$

其中，$q_{incident}$ 是入射辐射贡献的热流，$q_{emission}$ 是介质发射贡献的热流。

入射辐射热流

假设入射辐射强度为 $I_0$，则透射后的热流为：

$$q_{incident}={\int }_{4\pi }{I}_0{e}^{-\tau }\cos \theta d\Omega \approx q_0(1-\tau )$$

其中，$q_0$ 是入射热流。

介质发射热流

介质发射的热流为：

$$q_{emission}=2\pi {\int }_{-1}^1{\int }_0^{\tau }{I}_b{e}^{-(\tau -{\tau }^{\prime })/\mu }d{\tau }^{\prime }\mu d\mu$$

在光学薄极限下，上式简化为：

$$q_{emission}\approx 4\pi \tau I_b=4\sigma T^4\tau$$

这表明光学薄介质的发射热流与其光学厚度和温度的四次方成正比。

2.5 光学薄介质与壁面的辐射换热

考虑光学薄介质填充在两个平行平板之间。设两平板温度分别为 $T_1$ 和 $T_2$，介质温度为 $T_g$。在光学薄极限下，平板1到平板2的净辐射热流为：

$$q_{12}=\frac{\sigma (T_1^4-T_2^4)}{\frac{1}{{\epsilon }_1}+\frac{1}{{\epsilon }_2}-1}+{\epsilon }_g\sigma (T_g^4-T_2^4)$$

其中，${\epsilon }_g={\tau }_g={\kappa }_{g}L$ 是介质的发射率。

上式中，第一项是两平板间的直接辐射换热，第二项是介质发射对热流的贡献。在光学薄极限下，介质对两平板间直接辐射的衰减可以忽略，但介质本身的发射不可忽略。

---

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

3.1 直接积分法

对于光学薄介质，辐射传递方程的解可以表示为积分形式。考虑一维几何，辐射强度为：

$$I(\tau ,\mu )=I(0,\mu )e^{-\tau /\mu }+{\int }_0^{\tau }S({\tau }^{\prime })e^{-(\tau -{\tau }^{\prime })/\mu }\frac{d{\tau }^{\prime }}{\mu }$$

其中，$S$ 是源项，包括发射和散射贡献。

在光学薄极限下，指数项近似为1，上式简化为：

$$I(\tau ,\mu )\approx I(0,\mu )+{\int }_0^{\tau }S({\tau }^{\prime })\frac{d{\tau }^{\prime }}{\mu }$$

这一表达式可以通过数值积分求解。将光学厚度区间 $[0,\tau ]$ 离散为 $N$ 个小区间，则：

$$I(\tau ,\mu )\approx I(0,\mu )+\sum _{i=1}^N{S}_i\frac{\Delta {\tau }_i}{\mu }$$

其中，$S_i$ 是第 $i$ 个区间的平均源项。

3.2 逐次散射法

对于散射介质，可以采用逐次散射法求解。将辐射强度展开为各级散射贡献的叠加：

$$I=I^{(0)}+I^{(1)}+I^{(2)}+...$$

其中，$I^{(0)}$ 是直接透射和发射的贡献，$I^{(n)}$ 是经过 $n$ 次散射后的贡献。

零阶解（无散射）为：

$$I^{(0)}(\tau ,\mu )=I(0,\mu )e^{-\tau /\mu }+{\int }_0^{\tau }(1-\omega )I_b{e}^{-(\tau -{\tau }^{\prime })/\mu }\frac{d{\tau }^{\prime }}{\mu }$$

$n$ 阶散射贡献为：

$$I^{(n)}(\tau ,\mu )={\int }_0^{\tau }\frac{\omega }{4\pi }{\int }_{4\pi }\Phi I^{(n-1)}d{\Omega }^{\prime }{e}^{-(\tau -{\tau }^{\prime })/\mu }\frac{d{\tau }^{\prime }}{\mu }$$

在光学薄极限下，高阶散射贡献迅速衰减，通常只需计算前几阶即可得到足够精确的结果。

3.3 蒙特卡洛法

蒙特卡洛法是求解光学薄介质辐射问题的有效方法。其基本思想是追踪大量光子在介质中的传播过程，通过统计平均得到辐射特性。

对于光学薄介质，光子被吸收的概率很小，大部分光子能够穿透介质。蒙特卡洛模拟的主要步骤包括：

1. 光子发射：根据介质温度分布发射光子
2. 光子传输：计算光子的自由程，判断是否离开介质
3. 光子吸收：根据吸收概率决定是否吸收光子
4. 光子散射：对于散射介质，计算散射方向和新的传播方向
5. 统计累加：记录光子对辐射热流的贡献

光学薄条件下，光子平均自由程很大，模拟效率较高。

3.4 离散坐标法

离散坐标法（DOM）通过离散角度空间求解辐射传递方程。对于光学薄介质，角度离散可以采用较少的方向数即可获得较好的精度。

离散后的方程为：

$${\mu }_m\frac{d{I}_m}{d\tau }+I_m=(1-\omega )I_b+\frac{\omega }{4\pi }\sum _{m^{\prime }=1}^M{w}_{m^{\prime }}{\Phi }_{m{m}^{\prime }}{I}_{m^{\prime }}$$

其中，${\mu }_m$ 是离散方向，$w_{m^{\prime }}$ 是权重系数。

在光学薄极限下，方程可以沿每个离散方向独立求解，计算效率较高。

---

4. 应用案例

4.1 案例一：等温平板间光学薄气体辐射

问题描述

考虑两个无限大平行平板，其间填充光学薄气体。上平板温度 $T_1=1500$ K，下平板温度 $T_2=1000$ K，气体温度 $T_g=1200$ K。平板间距 $L=1$ m，气体吸收系数 $\kappa =0.01$ m⁻¹。计算两平板间的辐射换热，分析光学薄近似与精确解的差异。

数学模型

对于一维平行平板几何，辐射传递方程为：

$$\mu \frac{dI}{d\tau }+I=I_b$$

其中，$\tau =\kappa z$ 是光学厚度，$I_b=\sigma T^4/\pi$ 是黑体辐射强度。

边界条件为：

$$I^{+}(0)={\epsilon }_1{I}_{b1}+(1-{\epsilon }_1)I^{-}(0)$$
$$I^{-}({\tau }_L)={\epsilon }_2{I}_{b2}+(1-{\epsilon }_2)I^{+}({\tau }_L)$$

其中，$I^{+}$ 和 $I^{-}$ 分别表示正向和反向辐射强度。

精确解

辐射传递方程的精确解为：

$$I^{+}(\tau )=I^{+}(0)e^{-\tau }+{\int }_0^{\tau }{I}_b{e}^{-(\tau -{\tau }^{\prime })}d{\tau }^{\prime }$$
$$I^{-}(\tau )=I^{-}({\tau }_L)e^{-({\tau }_L-\tau )}+{\int }_{\tau }^{{\tau }_L}{I}_b{e}^{-({\tau }^{\prime }-\tau )}d{\tau }^{\prime }$$

对于均匀温度气体，上式简化为：

$$I^{+}(\tau )=I^{+}(0)e^{-\tau }+I_{bg}(1-e^{-\tau })$$
$$I^{-}(\tau )=I^{-}({\tau }_L)e^{-({\tau }_L-\tau )}+I_{bg}(1-e^{-({\tau }_L-\tau )})$$

光学薄近似

在光学薄极限下（${\tau }_L\ll 1$），指数项展开为：

$$e^{-\tau }\approx 1-\tau +O({\tau }^2)$$

辐射强度近似为：

$$I^{+}(\tau )\approx I^{+}(0)(1-\tau )+I_{bg}\tau$$
$$I^{-}(\tau )\approx I^{-}({\tau }_L)(1-({\tau }_L-\tau ))+I_{bg}({\tau }_L-\tau )$$

辐射热流计算

辐射热流为：

$$q=\pi (I^{+}-I^{-})$$

在光学薄极限下，平板1到平板2的净辐射热流为：

$$q_{12}=\frac{\sigma (T_1^4-T_2^4)}{\frac{1}{{\epsilon }_1}+\frac{1}{{\epsilon }_2}-1}+\kappa L\sigma (T_g^4-T_2^4)$$

仿真结果与分析

图1展示了等温平板间光学薄气体辐射的温度分布和辐射强度分布。从图中可以看出：

1. 光学薄近似与精确解在光学厚度小于0.1时吻合良好
2. 随着光学厚度增加，光学薄近似的误差逐渐增大
3. 气体发射对总辐射热流的贡献与光学厚度成正比

工程意义

本案例展示了光学薄近似在工程计算中的应用价值。对于低压气体、透明材料等光学薄介质，光学薄近似可以大幅简化计算，同时保持足够的精度。

4.2 案例二：非等温光学薄介质辐射换热

问题描述

考虑两个平行平板间的非等温光学薄气体。气体温度呈线性分布：$T(z)=T_1+(T_2-T_1)z/L$。上平板温度 $T_1=1500$ K，下平板温度 $T_2=1000$ K。平板间距 $L=1$ m，气体吸收系数 $\kappa =0.02$ m⁻¹。计算辐射热流分布，分析温度梯度对辐射换热的影响。

数学模型

对于非等温介质，辐射传递方程为：

$$\mu \frac{dI}{d\tau }+I=I_b(\tau )$$

其中，$I_b(\tau )=\sigma T^4(\tau )/\pi$ 是局部黑体辐射强度。

精确解

辐射强度的精确解为：

$$I^{+}(\tau )=I^{+}(0)e^{-\tau }+{\int }_0^{\tau }{I}_b({\tau }^{\prime })e^{-(\tau -{\tau }^{\prime })}d{\tau }^{\prime }$$
$$I^{-}(\tau )=I^{-}({\tau }_L)e^{-({\tau }_L-\tau )}+{\int }_{\tau }^{{\tau }_L}{I}_b({\tau }^{\prime })e^{-({\tau }^{\prime }-\tau )}d{\tau }^{\prime }$$

对于线性温度分布，$T(z)=T_1+(T_2-T_1)z/L$，黑体辐射强度为：

$$I_b(\tau )=\frac{\sigma }{\pi }{\left[T_1+(T_2-T_1)\frac{\tau }{{\tau }_L}\right]}^4$$

光学薄近似

在光学薄极限下，辐射强度近似为：

$$I^{+}(\tau )\approx I^{+}(0)+{\int }_0^{\tau }{I}_b({\tau }^{\prime })d{\tau }^{\prime }$$
$$I^{-}(\tau )\approx I^{-}({\tau }_L)+{\int }_{\tau }^{{\tau }_L}{I}_b({\tau }^{\prime })d{\tau }^{\prime }$$

辐射热流

辐射热流为：

$$q(\tau )=2\pi {\int }_0^1{I}^{+}(\tau ,\mu )\mu d\mu -2\pi {\int }_0^1{I}^{-}(\tau ,\mu )\mu d\mu$$

在光学薄极限下，上式简化为：

$$q(\tau )\approx q_{wall}+\pi {\int }_0^{\tau }{I}_b({\tau }^{\prime })d{\tau }^{\prime }-\pi {\int }_{\tau }^{{\tau }_L}{I}_b({\tau }^{\prime })d{\tau }^{\prime }$$

其中，$q_{wall}$ 是壁面辐射贡献的热流。

仿真结果与分析

图2展示了非等温光学薄介质的温度分布、辐射强度分布和辐射热流分布。从图中可以看出：

1. 温度梯度导致辐射热流沿程变化
2. 在光学薄极限下，辐射热流近似为线性分布
3. 温度梯度越大，辐射热流的变化越显著

工程应用

非等温光学薄介质模型广泛应用于：

• 燃气轮机燃烧室的热分析
• 火箭发动机喷管的辐射冷却
• 高温炉膛内的辐射换热

4.3 案例三：光学薄极限下的发射率计算

问题描述

计算不同光学厚度下介质层的发射率，验证光学薄近似 $\epsilon \approx \tau$ 的适用范围。考虑灰体介质，温度 $T=1200$ K，厚度 $L$ 从0.001 m变化到1.0 m，吸收系数 $\kappa =0.1$ m⁻¹。

理论分析

介质层的发射率定义为：

$$\epsilon =\frac{E}{E_b}=\frac{{\int }_0^{\infty }{\epsilon }_{\lambda }{E}_{b\lambda }d\lambda }{\sigma T^4}$$

对于灰体介质，光谱发射率与波长无关：

$${\epsilon }_{\lambda }=1-e^{-{\kappa }_{\lambda }L}=1-e^{-\tau }$$

因此，总发射率为：

$$\epsilon =1-e^{-\tau }$$

光学薄近似

在光学薄极限下（$\tau \ll 1$）：

$$\epsilon =1-e^{-\tau }\approx \tau -\frac{{\tau }^2}{2}+\frac{{\tau }^3}{6}-...\approx \tau$$

这一近似在 $\tau <0.1$ 时误差小于5%，在 $\tau <0.01$ 时误差小于0.5%。

误差分析

光学薄近似的相对误差为：

$$\delta =\frac{|{\epsilon }_{exact}-{\epsilon }_{approx}|}{{\epsilon }_{exact}}=\frac{|1-e^{-\tau }-\tau |}{1-e^{-\tau }}$$

对于小 $\tau$，展开得：

$$\delta \approx \frac{\tau }{2}$$

仿真结果

图3展示了发射率随光学厚度的变化曲线，以及光学薄近似的误差分析。结果表明：

1. 当 $\tau <0.1$ 时，光学薄近似误差小于5%
2. 当 $\tau <0.01$ 时，光学薄近似误差小于0.5%
3. 随着光学厚度增加，光学薄近似的误差迅速增大

工程指导

本案例为工程计算提供了重要指导：

• 当 $\tau <0.1$ 时，可以使用光学薄近似快速估算发射率
• 当 $0.1<\tau <1.0$ 时，建议使用精确公式
• 当 $\tau >1.0$ 时，介质接近黑体，发射率接近1

4.4 案例四：多层光学薄介质辐射

问题描述

考虑三层光学薄气体介质填充在两个平行平板之间。每层气体的吸收系数和温度不同：

• 第一层：${\kappa }_1=0.01$ m⁻¹，$T_1=1000$ K，厚度 $L_1=0.3$ m
• 第二层：${\kappa }_2=0.02$ m⁻¹，$T_2=1200$ K，厚度 $L_2=0.4$ m
• 第三层：${\kappa }_3=0.015$ m⁻¹，$T_3=1100$ K，厚度 $L_3=0.3$ m

上平板温度 $T_{w1}=1500$ K，下平板温度 $T_{w2}=800$ K。计算通过多层介质的辐射热流。

数学模型

对于多层介质，每层内的辐射传递方程为：

$$\mu \frac{d{I}^{(k)}}{d\tau }+I^{(k)}=I_b^{(k)}$$

其中，$k=1,2,3$ 表示层数。

层间界面处的连续性条件为：

$$I^{(k)}({\tau }_k^{-})=I^{(k+1)}({\tau }_k^{+})$$

光学薄近似

在光学薄极限下，每层的辐射强度为：

$$I^{(k)}(\tau )\approx I^{(k)}(0)+{\int }_0^{\tau }{I}_b^{(k)}d{\tau }^{\prime }$$

多层介质的总透射率为：

$${\tau }_{total}=\prod _{k=1}^3{e}^{-{\tau }_k}\approx 1-\sum _{k=1}^3{\tau }_k=1-({\tau }_1+{\tau }_2+{\tau }_3)$$

辐射热流计算

通过多层介质的总辐射热流为：

$$q=q_{wall}+\sum _{k=1}^3{q}_{emission}^{(k)}$$

其中，$q_{wall}$ 是壁面辐射贡献，$q_{emission}^{(k)}$ 是第 $k$ 层介质的发射贡献：

$$q_{emission}^{(k)}=4\sigma T_k^4{\tau }_k$$

仿真结果与分析

图4展示了多层光学薄介质的辐射特性。结果表明：

1. 多层介质的辐射特性可以线性叠加
2. 温度最高的层对辐射热流的贡献最大
3. 光学薄近似与精确解吻合良好

工程应用

多层光学薄介质模型应用于：

• 多层隔热材料（MLI）的热分析
• 大气层辐射传输模型
• 复合材料的辐射特性计算

4.5 案例五：光学薄介质与壁面辐射耦合

问题描述

考虑光学薄气体与壁面的辐射耦合问题。圆柱形燃烧室，直径 $D=2$ m，长度 $L=4$ m。燃烧产物温度 $T_g=1800$ K，壁面温度 $T_w=1200$ K。燃烧产物吸收系数 $\kappa =0.05$ m⁻¹。计算壁面净辐射热流和气体对壁面的辐射加热。

数学模型

对于圆柱形几何，采用柱坐标系 $(r,\theta ,z)$。辐射传递方程为：

$$\frac{1}{c}\frac{\partial I}{\partial t}+\mathbf{s}\cdot \nabla I=-\kappa I+\kappa I_b$$

在稳态条件下，简化为：

$$\mathbf{s}\cdot \nabla I=-\kappa I+\kappa I_b$$

光学薄近似

在光学薄极限下，辐射强度为：

$$I(\mathbf{r},\mathbf{s})\approx I_0({\mathbf{r}}_w,\mathbf{s})+{\int }_0^s\kappa I_{b}d{s}^{\prime }$$

其中，$I_0$ 是壁面发射的辐射强度，$s$ 是沿射线方向的距离。

壁面辐射热流

壁面处的辐射热流为：

$$q_w={\int }_{\mathbf{s}\cdot \mathbf{n}<0}{I}_w|\mathbf{s}\cdot \mathbf{n}|d\Omega -{\int }_{\mathbf{s}\cdot \mathbf{n}>0}{I}_g|\mathbf{s}\cdot \mathbf{n}|d\Omega$$

其中，$\mathbf{n}$ 是壁面外法向，$I_w$ 是壁面发射强度，$I_g$ 是气体发射强度。

蒙特卡洛求解

采用蒙特卡洛方法追踪大量光子：

1. 从壁面和气体发射光子
2. 计算光子自由程：$l=-\ln (R)/\kappa$，$R$ 是随机数
3. 判断光子是否被吸收或离开系统
4. 统计壁面接收的辐射能量

仿真结果与分析

图5展示了圆柱形燃烧室中光学薄气体与壁面的辐射耦合特性。结果表明：

1. 气体对壁面的辐射加热沿轴向分布不均匀
2. 壁面中心区域接收的辐射热流最大
3. 光学薄近似与蒙特卡洛模拟结果吻合良好

工程意义

本案例展示了光学薄近似在复杂几何中的应用。对于燃气轮机燃烧室、工业炉膛等设备，光学薄模型可以快速估算辐射热负荷，为热设计提供重要参考。

---

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

5.1 航天器热控系统

航天器在太空中运行时，外部没有大气对流散热，辐射成为主要的散热方式。航天器表面通常包覆多层隔热材料（MLI），层间填充低气压气体。

MLI的辐射特性

MLI层间气体处于低压状态（$10^{-3}\sim 10^{-5}$ Pa），分子平均自由程很大，光学厚度很小。层间辐射换热可以近似为光学薄介质辐射。

层间辐射热流为：

$$q=\frac{\sigma (T_1^4-T_2^4)}{\frac{1}{{\epsilon }_1}+\frac{1}{{\epsilon }_2}-1}+{\epsilon }_g\sigma (T_g^4-T_2^4)$$

其中，${\epsilon }_g={\kappa }_{g}d$ 是层间气体的发射率，$d$ 是层间距。

热设计考虑

• 降低层间气压，减小气体辐射贡献
• 使用高反射率材料（如铝箔），提高表面反射率
• 优化层数和层间距，平衡隔热性能和重量

5.2 燃烧设备辐射分析

燃气轮机、火箭发动机、工业炉等燃烧设备中，高温燃烧产物（CO₂、H₂O、碳烟等）会产生强烈的辐射。

燃烧产物的辐射特性

燃烧产物的吸收系数随波长和温度变化。在某些波段（如CO₂的2.7 μm和4.3 μm带），吸收系数较高；在其他波段，吸收系数较低，可能处于光学薄状态。

光学薄近似应用

对于光学薄波段，辐射热损失可以简化为：

$$q_{rad}={\epsilon }_{eff}\sigma T_g^4$$

其中，${\epsilon }_{eff}={\sum }_i{\kappa }_i{L}_i$ 是有效发射率，${\kappa }_i$ 是第 $i$ 种组分的吸收系数。

工程计算

光学薄模型可用于：

• 快速估算燃烧室壁面热负荷
• 优化燃烧室冷却设计
• 预测燃烧效率（辐射热损失）

5.3 透明材料热分析

玻璃、晶体、陶瓷等透明材料在某些波段具有很低的吸收系数，可以近似为光学薄介质。

玻璃幕墙的热分析

建筑玻璃幕墙对可见光透明，但对红外辐射有一定吸收。在夏季，太阳辐射透过玻璃进入室内，导致室内温度升高；在冬季，室内热量通过玻璃向外辐射，导致热量损失。

采用光学薄模型，可以计算玻璃的辐射特性：

$${\tau }_{vis}\approx 1-{\kappa }_{vis}d$$
$${\alpha }_{IR}\approx {\kappa }_{IR}d$$

其中，${\tau }_{vis}$ 是可见光透射率，${\alpha }_{IR}$ 是红外吸收率，$d$ 是玻璃厚度。

节能设计

• 使用低辐射（Low-E）涂层，提高红外反射率
• 采用双层或三层玻璃，增加隔热性能
• 优化玻璃厚度和气体填充层

5.4 大气辐射传输

地球大气层在可见光波段近似为光学薄介质，但在红外波段具有较强的吸收和发射。

大气辐射模型

大气辐射传输方程为：

$$\mu \frac{dI}{d\tau }+I=I_b$$

其中，光学厚度 $\tau$ 沿垂直方向积分：

$$\tau (z)={\int }_z^{\infty }\kappa (z^{\prime })d{z}^{\prime }$$

遥感应用

卫星遥感通过测量大气层顶部的辐射强度，反演地表和大气的物理参数。在光学薄近似下，辐射强度为：

$$I(\tau =0)\approx I_{surface}+{\int }_0^{\infty }\kappa I_{b}dz$$

这一模型广泛应用于气象卫星、资源卫星等遥感数据处理。

---

6. 光学薄近似的精度与适用范围

6.1 误差分析

光学薄近似的误差主要来源于对指数项的截断。精确解中的指数项为：

$$e^{-\tau }=1-\tau +\frac{{\tau }^2}{2}-\frac{{\tau }^3}{6}+...$$

光学薄近似只保留线性项：

$$e^{-\tau }\approx 1-\tau$$

截断误差为：

$$\delta =\frac{{\tau }^2}{2}-\frac{{\tau }^3}{6}+...\approx \frac{{\tau }^2}{2}$$

相对误差为：

$$\frac{\delta }{e^{-\tau }}\approx \frac{{\tau }^2/2}{1-\tau }\approx \frac{{\tau }^2}{2}$$

6.2 适用范围判据

光学薄近似的适用范围取决于所需的精度：

光学厚度 $\tau$	相对误差	适用性
$<0.01$	$<0.005\%$	极高精度
$<0.1$	$<0.5\%$	高精度
$<0.3$	$<5\%$	工程精度
$<0.5$	$<15\%$	粗略估算

6.3 改进的光学薄近似

对于中等光学厚度（$0.1<\tau <1.0$），可以采用改进的光学薄近似：

二阶近似

$$e^{-\tau }\approx 1-\tau +\frac{{\tau }^2}{2}$$

这一近似将适用范围扩展到 $\tau <0.5$，相对误差小于5%。

指数修正

$$e^{-\tau }\approx \frac{1}{1+\tau +{\tau }^2/2}$$

这一近似在 $\tau <1.0$ 时误差小于10%。

Padé近似

$$e^{-\tau }\approx \frac{1-\tau /2}{1+\tau /2}$$

这一近似在 $\tau <2.0$ 时误差小于10%。

6.4 数值验证

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

验证结果

• 光学薄近似在 $\tau <0.1$ 时误差小于1%
• 二阶近似在 $\tau <0.5$ 时误差小于5%
• Padé近似在 $\tau <2.0$ 时误差小于10%

---

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

7.1 非灰体光学薄介质

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

波谱积分

总辐射特性通过对波谱积分得到：

$$\epsilon =\frac{{\int }_0^{\infty }{\epsilon }_{\lambda }{E}_{b\lambda }d\lambda }{\sigma T^4}$$

对于光学薄介质：

$${\epsilon }_{\lambda }=1-e^{-{\kappa }_{\lambda }L}\approx {\kappa }_{\lambda }L$$

因此：

$$\epsilon \approx L\frac{{\int }_0^{\infty }{\kappa }_{\lambda }{E}_{b\lambda }d\lambda }{\sigma T^4}$$

带模型

对于气体介质，吸收系数在某些波段（吸收带）较高，在其他波段较低。可以采用带模型简化计算：

$$\epsilon =\sum _i\Delta {\lambda }_i{\kappa }_{i}L\frac{E_{b{\lambda }_i}}{\sigma T^4}$$

其中，$\Delta {\lambda }_i$ 是第 $i$ 个吸收带的宽度，${\kappa }_i$ 是平均吸收系数。

7.2 各向异性散射

光学薄介质中的散射可能是各向异性的。散射相函数描述了散射的方向分布。

相函数展开

散射相函数可以展开为勒让德多项式：

$$\Phi (\mu )=\sum _{l=0}^N(2l+1)g_l{P}_l(\mu )$$

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

光学薄近似

在光学薄极限下，散射对辐射强度的贡献为：

$$I_{scat}={\int }_0^{\tau }\frac{\omega }{4\pi }{\int }_{4\pi }\Phi Id{\Omega }^{\prime }d{\tau }^{\prime }$$

对于各向同性散射（$\Phi =1$）：

$$I_{scat}=\omega \tau G/4\pi$$

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

7.3 瞬态光学薄辐射

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

瞬态辐射传递方程

$$\frac{1}{c}\frac{\partial I}{\partial t}+\mathbf{s}\cdot \nabla I=-\kappa I+\kappa I_b$$

光学薄近似

在光学薄极限下，辐射传播时间很短，可以忽略时间延迟效应。瞬态方程简化为：

$$\frac{\partial I}{\partial t}\approx c\kappa (I_b-I)$$

这一方程描述了辐射强度趋近于平衡值的弛豫过程，弛豫时间为 $1/(c\kappa )$。

7.4 多维光学薄辐射

实际工程问题通常涉及多维几何。在光学薄极限下，多维辐射问题可以简化为一维问题的叠加。

二维问题

对于二维问题，辐射强度沿射线方向传播：

$$I(s)=I(0)+{\int }_0^s\kappa I_{b}d{s}^{\prime }$$

辐射热流为：

$$\mathbf{q}={\int }_{4\pi }I\mathbf{s}d\Omega$$

在光学薄极限下，辐射热流可以分解为各方向的贡献。

蒙特卡洛方法

对于复杂几何，蒙特卡洛方法是求解光学薄辐射问题的有效工具。由于光学薄介质中光子平均自由程很大，模拟效率较高。

---

8. 总结与展望

8.1 主要内容回顾

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

1. 光学薄介质的基本概念：介绍了光学厚度的定义、光学薄极限的数学表述和光学薄介质辐射的特征。

2. 理论基础：建立了光学薄极限下的辐射传递方程简化形式，推导了发射率、吸收率和辐射热流的计算公式。

3. 数值方法：介绍了直接积分法、逐次散射法、蒙特卡洛法和离散坐标法等求解光学薄辐射问题的方法。

4. 应用案例：通过五个典型案例，展示了光学薄近似在等温/非等温介质、发射率计算、多层介质和壁面耦合问题中的应用。

5. 工程应用：讨论了光学薄模型在航天器热控、燃烧设备、透明材料和大气辐射传输等领域的应用。

6. 精度分析：分析了光学薄近似的误差来源和适用范围，提出了改进的近似方法。

8.2 关键结论

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

1. 光学薄近似适用条件：当光学厚度 $\tau <0.1$ 时，光学薄近似误差小于5%，可以安全使用；当 $\tau <0.01$ 时，误差小于0.5%，近似精度极高。

2. 发射率与光学厚度的关系：光学薄介质的发射率近似等于光学厚度，$\epsilon \approx \tau =\kappa L$。

3. 线性叠加特性：光学薄介质中，不同位置、不同层的发射贡献可以线性叠加，大大简化了计算。

4. 工程应用价值：光学薄模型在航天器热控、燃烧设备、透明材料等领域有广泛应用，可以快速估算辐射热负荷。

8.3 未来发展方向

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

1. 非灰体模型：发展更精确的非灰体光学薄模型，考虑波谱依赖的吸收特性。

2. 多物理场耦合：研究光学薄辐射与流动、化学反应的耦合效应。

3. 不确定性量化：考虑辐射物性参数的不确定性，建立鲁棒的设计方法。

4. 机器学习应用：利用机器学习方法加速光学薄辐射计算，提高工程应用效率。

5. 实验验证：开展更多的实验研究，验证光学薄模型的精度和适用范围。

8.4 学习建议

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

1. 理论基础：系统学习辐射传递方程、光学厚度、发射率等基本概念。

2. 数值方法：掌握蒙特卡洛法、离散坐标法等数值求解方法。

3. 编程实践：通过编写Python程序，实现光学薄辐射的仿真计算。

4. 工程应用：结合实际工程问题，应用光学薄模型进行热分析和设计。

5. 文献阅读：阅读最新的研究文献，了解该领域的前沿进展。

附录：Python仿真代码

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

1. 光学薄辐射传递方程求解器：实现了一维、二维光学薄辐射的数值求解。

2. 发射率计算模块：计算不同光学厚度下的介质发射率，验证光学薄近似。

3. 蒙特卡洛模拟器：实现光子追踪算法，用于复杂几何的辐射计算。

4. 可视化工具：生成辐射强度分布、热流分布、发射率曲线等图表。

5. GIF动画生成：创建辐射传播过程的动态可视化。

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

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

本程序实现光学薄介质辐射换热的数值模拟，包括：
1. 等温平板间光学薄气体辐射
2. 非等温光学薄介质辐射换热
3. 光学薄极限下的发射率计算
4. 多层光学薄介质辐射
5. 光学薄介质与壁面辐射耦合

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

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle, Circle, FancyBboxPatch
from matplotlib.collections import LineCollection
import matplotlib.animation as animation
from scipy import integrate
from scipy.integrate import trapezoid
from scipy.special import exp1
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^2·K^4)

print("=" * 80)
print("主题019：光学薄介质辐射换热 - Python仿真程序")
print("=" * 80)

# ==============================================================================
# 案例一：等温平板间光学薄气体辐射
# ==============================================================================
print("\n" + "=" * 80)
print("案例一：等温平板间光学薄气体辐射")
print("=" * 80)

def case1_isothermal_plates():
    """
    案例一：等温平板间光学薄气体辐射
    
    计算两个平行平板间填充光学薄气体时的辐射换热
    比较精确解与光学薄近似
    """
    print("\n[案例一] 等温平板间光学薄气体辐射")
    print("-" * 60)
    
    # 物理参数
    T1 = 1500.0  # 上平板温度, K
    T2 = 1000.0  # 下平板温度, K
    Tg = 1200.0  # 气体温度, K
    L = 1.0      # 平板间距, m
    kappa = 0.01 # 气体吸收系数, m^-1
    eps1 = 0.8   # 上平板发射率
    eps2 = 0.8   # 下平板发射率
    
    # 计算光学厚度
    tau_L = kappa * L
    print(f"光学厚度 τ = {tau_L:.4f}")
    print(f"光学薄条件 (τ << 1): {'满足' if tau_L < 0.1 else '不满足'}")
    
    # 黑体辐射强度
    Ib1 = sigma * T1**4 / np.pi
    Ib2 = sigma * T2**4 / np.pi
    Ibg = sigma * Tg**4 / np.pi
    
    # 空间离散
    Nz = 100
    z = np.linspace(0, L, Nz)
    tau = kappa * z
    
    # 精确解：求解辐射传递方程
    # 使用离散坐标法（S4）
    mu_pos = np.array([0.2958759, 0.9082483])  # 正向方向余弦
    w_pos = np.array([0.5235988, 0.5235988])   # 权重
    
    # 初始化辐射强度
    I_plus = np.zeros((Nz, len(mu_pos)))  # 正向辐射强度
    I_minus = np.zeros((Nz, len(mu_pos))) # 反向辐射强度
    
    # 边界条件
    # 上边界 (z=0)
    I_plus[0, :] = eps1 * Ib1
    # 下边界 (z=L)
    I_minus[-1, :] = eps2 * Ib2
    
    # 迭代求解
    max_iter = 1000
    tol = 1e-6
    
    for iteration in range(max_iter):
        I_plus_old = I_plus.copy()
        I_minus_old = I_minus.copy()
        
        # 正向扫描
        for i in range(1, Nz):
            dz = z[i] - z[i-1]
            dtau = kappa * dz
            for m in range(len(mu_pos)):
                # 指数格式
                I_plus[i, m] = I_plus[i-1, m] * np.exp(-dtau/mu_pos[m]) + \
                               Ibg * (1 - np.exp(-dtau/mu_pos[m]))
        
        # 反向扫描
        for i in range(Nz-2, -1, -1):
            dz = z[i+1] - z[i]
            dtau = kappa * dz
            for m in range(len(mu_pos)):
                I_minus[i, m] = I_minus[i+1, m] * np.exp(-dtau/mu_pos[m]) + \
                                Ibg * (1 - np.exp(-dtau/mu_pos[m]))
        
        # 检查收敛
        error = np.max(np.abs(I_plus - I_plus_old)) + np.max(np.abs(I_minus - I_minus_old))
        if error < tol:
            print(f"迭代收敛于第 {iteration+1} 次迭代")
            break
    
    # 计算辐射热流
    q_plus = 2 * np.pi * np.sum(I_plus * mu_pos * w_pos, axis=1)
    q_minus = 2 * np.pi * np.sum(I_minus * mu_pos * w_pos, axis=1)
    q_exact = q_plus - q_minus
    
    # 光学薄近似
    # 壁面间直接辐射换热
    q_wall = sigma * (T1**4 - T2**4) / (1/eps1 + 1/eps2 - 1)
    # 气体发射贡献
    q_gas_approx = kappa * L * sigma * (Tg**4 - T2**4)
    q_approx = q_wall + q_gas_approx
    
    print(f"\n辐射热流计算结果:")
    print(f"  精确解 (DOM): {q_exact[0]:.2f} W/m²")
    print(f"  光学薄近似: {q_approx:.2f} W/m²")
    print(f"  相对误差: {abs(q_exact[0] - q_approx)/q_exact[0]*100:.2f}%")
    
    # 可视化
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 图1: 辐射强度分布
    ax1 = axes[0, 0]
    ax1.plot(z, I_plus[:, 0], 'r-', linewidth=2, label=r'$I^+$ ($\mu$=' + f'{mu_pos[0]:.3f})')
    ax1.plot(z, I_plus[:, 1], 'r--', linewidth=2, label=r'$I^+$ ($\mu$=' + f'{mu_pos[1]:.3f})')
    ax1.plot(z, I_minus[:, 0], 'b-', linewidth=2, label=r'$I^-$ ($\mu$=' + f'{mu_pos[0]:.3f})')
    ax1.plot(z, I_minus[:, 1], 'b--', linewidth=2, label=r'$I^-$ ($\mu$=' + f'{mu_pos[1]:.3f})')
    ax1.axhline(y=Ibg, color='g', linestyle=':', linewidth=2, label=f'$I_b$ (气体) = {Ibg:.1f}')
    ax1.set_xlabel('位置 z (m)', fontsize=12)
    ax1.set_ylabel('辐射强度 (W/(m²·sr))', fontsize=12)
    ax1.set_title('案例一：辐射强度分布', fontsize=14, fontweight='bold')
    ax1.legend(fontsize=10)
    ax1.grid(True, alpha=0.3)
    
    # 图2: 辐射热流分布
    ax2 = axes[0, 1]
    ax2.plot(z, q_exact, 'b-', linewidth=2.5, label='精确解 (DOM)')
    ax2.axhline(y=q_approx, color='r', linestyle='--', linewidth=2, label='光学薄近似')
    ax2.set_xlabel('位置 z (m)', fontsize=12)
    ax2.set_ylabel('辐射热流 (W/m²)', fontsize=12)
    ax2.set_title('案例一：辐射热流分布', fontsize=14, fontweight='bold')
    ax2.legend(fontsize=11)
    ax2.grid(True, alpha=0.3)
    
    # 图3: 温度分布
    ax3 = axes[1, 0]
    ax3.plot([0, 0], [800, T1], 'r-', linewidth=4, label=f'上平板 T₁={T1}K')
    ax3.plot([L, L], [800, T2], 'b-', linewidth=4, label=f'下平板 T₂={T2}K')
    ax3.plot([0, L], [Tg, Tg], 'g--', linewidth=2, label=f'气体 Tg={Tg}K')
    ax3.fill_between([0, L], 800, Tg, alpha=0.2, color='green')
    ax3.set_xlabel('位置 z (m)', fontsize=12)
    ax3.set_ylabel('温度 (K)', fontsize=12)
    ax3.set_title('案例一：温度分布', fontsize=14, fontweight='bold')
    ax3.set_xlim(-0.1, L+0.1)
    ax3.set_ylim(800, 1600)
    ax3.legend(fontsize=11)
    ax3.grid(True, alpha=0.3)
    
    # 图4: 误差分析 - 不同光学厚度下的比较
    ax4 = axes[1, 1]
    tau_range = np.linspace(0.001, 0.5, 50)
    q_exact_range = []
    q_approx_range = []
    
    for tau_test in tau_range:
        kappa_test = tau_test / L
        # 简化的精确解
        q_e = sigma * (T1**4 - T2**4) / (1/eps1 + 1/eps2 - 1) * np.exp(-tau_test) + \
              sigma * Tg**4 * (1 - np.exp(-tau_test))
        q_a = sigma * (T1**4 - T2**4) / (1/eps1 + 1/eps2 - 1) + \
              tau_test * sigma * (Tg**4 - T2**4)
        q_exact_range.append(q_e)
        q_approx_range.append(q_a)
    
    q_exact_range = np.array(q_exact_range)
    q_approx_range = np.array(q_approx_range)
    error_percent = np.abs(q_exact_range - q_approx_range) / q_exact_range * 100
    
    ax4.semilogy(tau_range, error_percent, 'b-', linewidth=2)
    ax4.axvline(x=0.1, color='r', linestyle='--', linewidth=2, label='τ = 0.1 (5%误差线)')
    ax4.axhline(y=5, color='r', linestyle='--', linewidth=2)
    ax4.set_xlabel('光学厚度 τ', fontsize=12)
    ax4.set_ylabel('相对误差 (%)', fontsize=12)
    ax4.set_title('案例一：光学薄近似误差分析', fontsize=14, fontweight='bold')
    ax4.legend(fontsize=11)
    ax4.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('case1_isothermal_plates.png', dpi=150, bbox_inches='tight')
    print("  图片已保存: case1_isothermal_plates.png")
    plt.close()
    
    return q_exact[0], q_approx

# ==============================================================================
# 案例二：非等温光学薄介质辐射换热
# ==============================================================================
print("\n" + "=" * 80)
print("案例二：非等温光学薄介质辐射换热")
print("=" * 80)

def case2_nonisothermal_medium():
    """
    案例二：非等温光学薄介质辐射换热
    
    计算具有线性温度分布的光学薄介质的辐射换热
    """
    print("\n[案例二] 非等温光学薄介质辐射换热")
    print("-" * 60)
    
    # 物理参数
    T1 = 1500.0  # 上平板温度, K
    T2 = 1000.0  # 下平板温度, K
    L = 1.0      # 平板间距, m
    kappa = 0.02 # 气体吸收系数, m^-1
    eps1 = 0.8   # 上平板发射率
    eps2 = 0.8   # 下平板发射率
    
    # 光学厚度
    tau_L = kappa * L
    print(f"光学厚度 τ = {tau_L:.4f}")
    
    # 空间离散
    Nz = 100
    z = np.linspace(0, L, Nz)
    tau = kappa * z
    
    # 线性温度分布
    T = T1 + (T2 - T1) * z / L
    Ib = sigma * T**4 / np.pi
    
    # 离散坐标方向
    mu_pos = np.array([0.2958759, 0.9082483])
    w_pos = np.array([0.5235988, 0.5235988])
    
    # 初始化
    I_plus = np.zeros((Nz, len(mu_pos)))
    I_minus = np.zeros((Nz, len(mu_pos)))
    
    # 边界条件
    Ib1 = sigma * T1**4 / np.pi
    Ib2 = sigma * T2**4 / np.pi
    I_plus[0, :] = eps1 * Ib1
    I_minus[-1, :] = eps2 * Ib2
    
    # 迭代求解
    max_iter = 1000
    tol = 1e-6
    
    for iteration in range(max_iter):
        I_plus_old = I_plus.copy()
        I_minus_old = I_minus.copy()
        
        # 正向扫描
        for i in range(1, Nz):
            dz = z[i] - z[i-1]
            dtau = kappa * dz
            for m in range(len(mu_pos)):
                # 使用平均黑体辐射强度
                Ib_avg = 0.5 * (Ib[i] + Ib[i-1])
                I_plus[i, m] = I_plus[i-1, m] * np.exp(-dtau/mu_pos[m]) + \
                               Ib_avg * (1 - np.exp(-dtau/mu_pos[m]))
        
        # 反向扫描
        for i in range(Nz-2, -1, -1):
            dz = z[i+1] - z[i]
            dtau = kappa * dz
            for m in range(len(mu_pos)):
                Ib_avg = 0.5 * (Ib[i] + Ib[i+1])
                I_minus[i, m] = I_minus[i+1, m] * np.exp(-dtau/mu_pos[m]) + \
                                Ib_avg * (1 - np.exp(-dtau/mu_pos[m]))
        
        error = np.max(np.abs(I_plus - I_plus_old)) + np.max(np.abs(I_minus - I_minus_old))
        if error < tol:
            print(f"迭代收敛于第 {iteration+1} 次迭代")
            break
    
    # 计算辐射热流
    q_plus = 2 * np.pi * np.sum(I_plus * mu_pos * w_pos, axis=1)
    q_minus = 2 * np.pi * np.sum(I_minus * mu_pos * w_pos, axis=1)
    q_exact = q_plus - q_minus
    
    # 光学薄近似
    # 计算发射贡献的积分
    emission_integral = np.zeros(Nz)
    for i in range(Nz):
        emission_integral[i] = trapezoid(kappa * sigma * T[:i+1]**4, z[:i+1]) if i > 0 else 0
    
    # 壁面辐射
    q_wall = sigma * (T1**4 - T2**4) / (1/eps1 + 1/eps2 - 1)
    
    # 光学薄近似热流
    q_approx = q_wall + 2 * kappa * sigma * (T**4 - T2**4) * (L - z) - \
               2 * kappa * sigma * (T1**4 - T**4) * z
    
    print(f"\n辐射热流:")
    print(f"  上边界 (z=0): 精确解={q_exact[0]:.2f}, 近似={q_wall:.2f} W/m²")
    print(f"  下边界 (z=L): 精确解={q_exact[-1]:.2f}, 近似={q_wall:.2f} W/m²")
    
    # 可视化
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 图1: 温度分布
    ax1 = axes[0, 0]
    ax1.plot(z, T, 'b-', linewidth=2.5)
    ax1.set_xlabel('位置 z (m)', fontsize=12)
    ax1.set_ylabel('温度 (K)', fontsize=12)
    ax1.set_title('案例二：线性温度分布', fontsize=14, fontweight='bold')
    ax1.grid(True, alpha=0.3)
    
    # 图2: 黑体辐射强度分布
    ax2 = axes[0, 1]
    ax2.plot(z, Ib, 'r-', linewidth=2.5)
    ax2.set_xlabel('位置 z (m)', fontsize=12)
    ax2.set_ylabel('$I_b$ (W/(m²·sr))', fontsize=12)
    ax2.set_title('案例二：黑体辐射强度分布', fontsize=14, fontweight='bold')
    ax2.grid(True, alpha=0.3)
    
    # 图3: 辐射热流分布
    ax3 = axes[1, 0]
    ax3.plot(z, q_exact, 'b-', linewidth=2.5, label='精确解 (DOM)')
    ax3.plot(z, q_approx, 'r--', linewidth=2, label='光学薄近似')
    ax3.set_xlabel('位置 z (m)', fontsize=12)
    ax3.set_ylabel('辐射热流 (W/m²)', fontsize=12)
    ax3.set_title('案例二：辐射热流分布', fontsize=14, fontweight='bold')
    ax3.legend(fontsize=11)
    ax3.grid(True, alpha=0.3)
    
    # 图4: 辐射强度分布 (特定方向)
    ax4 = axes[1, 1]
    ax4.plot(z, I_plus[:, 0], 'r-', linewidth=2, label=r'$I^+$ ($\mu$=' + f'{mu_pos[0]:.3f})')
    ax4.plot(z, I_plus[:, 1], 'r--', linewidth=2, label=r'$I^+$ ($\mu$=' + f'{mu_pos[1]:.3f})')
    ax4.plot(z, I_minus[:, 0], 'b-', linewidth=2, label=r'$I^-$ ($\mu$=' + f'{mu_pos[0]:.3f})')
    ax4.plot(z, I_minus[:, 1], 'b--', linewidth=2, label=r'$I^-$ ($\mu$=' + f'{mu_pos[1]:.3f})')
    ax4.plot(z, Ib, 'g:', linewidth=2, label='$I_b(z)$')
    ax4.set_xlabel('位置 z (m)', fontsize=12)
    ax4.set_ylabel('辐射强度 (W/(m²·sr))', fontsize=12)
    ax4.set_title('案例二：辐射强度分布', fontsize=14, fontweight='bold')
    ax4.legend(fontsize=10)
    ax4.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('case2_nonisothermal_medium.png', dpi=150, bbox_inches='tight')
    print("  图片已保存: case2_nonisothermal_medium.png")
    plt.close()
    
    return q_exact, q_approx

# ==============================================================================
# 案例三：光学薄极限下的发射率计算
# ==============================================================================
print("\n" + "=" * 80)
print("案例三：光学薄极限下的发射率计算")
print("=" * 80)

def case3_emissivity_calculation():
    """
    案例三：光学薄极限下的发射率计算
    
    计算不同光学厚度下的介质发射率，验证光学薄近似
    """
    print("\n[案例三] 光学薄极限下的发射率计算")
    print("-" * 60)
    
    # 物理参数
    T = 1200.0   # 介质温度, K
    kappa = 0.1  # 吸收系数, m^-1
    
    # 光学厚度范围
    tau_range = np.logspace(-3, 1, 100)  # 从0.001到10
    L_range = tau_range / kappa
    
    # 精确发射率
    epsilon_exact = 1 - np.exp(-tau_range)
    
    # 光学薄近似
    epsilon_approx = tau_range
    
    # 二阶近似
    epsilon_second_order = tau_range - tau_range**2 / 2
    
    # Padé近似
    epsilon_pade = (1 - tau_range/2) / (1 + tau_range/2)
    epsilon_pade = 1 - epsilon_pade  # 转换为发射率
    
    # 计算相对误差
    error_approx = np.abs(epsilon_exact - epsilon_approx) / epsilon_exact * 100
    error_second = np.abs(epsilon_exact - epsilon_second_order) / epsilon_exact * 100
    error_pade = np.abs(epsilon_exact - epsilon_pade) / epsilon_exact * 100
    
    # 打印关键结果
    print("\n不同光学厚度下的发射率:")
    test_tau = [0.001, 0.01, 0.05, 0.1, 0.3, 0.5, 1.0]
    for tau in test_tau:
        eps_ex = 1 - np.exp(-tau)
        eps_ap = tau
        err = abs(eps_ex - eps_ap) / eps_ex * 100
        print(f"  τ = {tau:.3f}: 精确值={eps_ex:.5f}, 近似值={eps_ap:.5f}, 误差={err:.3f}%")
    
    # 可视化
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 图1: 发射率随光学厚度变化
    ax1 = axes[0, 0]
    ax1.plot(tau_range, epsilon_exact, 'b-', linewidth=2.5, label='精确解: $\\varepsilon = 1 - e^{-\\tau}$')
    ax1.plot(tau_range, epsilon_approx, 'r--', linewidth=2, label='光学薄近似: $\\varepsilon \\approx \\tau$')
    ax1.plot(tau_range, epsilon_second_order, 'g:', linewidth=2, label='二阶近似: $\\varepsilon \\approx \\tau - \\tau^2/2$')
    ax1.axvline(x=0.1, color='k', linestyle='-.', alpha=0.5, label='$\\tau = 0.1$')
    ax1.set_xlabel('光学厚度 τ', fontsize=12)
    ax1.set_ylabel('发射率 ε', fontsize=12)
    ax1.set_title('案例三：发射率随光学厚度变化', fontsize=14, fontweight='bold')
    ax1.set_xlim([0, 2])
    ax1.set_ylim([0, 1])
    ax1.legend(fontsize=10)
    ax1.grid(True, alpha=0.3)
    
    # 图2: 相对误差
    ax2 = axes[0, 1]
    ax2.semilogy(tau_range, error_approx, 'r-', linewidth=2, label='光学薄近似误差')
    ax2.semilogy(tau_range, error_second, 'g-', linewidth=2, label='二阶近似误差')
    ax2.semilogy(tau_range, error_pade, 'm-', linewidth=2, label='Padé近似误差')
    ax2.axvline(x=0.1, color='k', linestyle='-.', alpha=0.5)
    ax2.axhline(y=5, color='r', linestyle='--', alpha=0.5, label='5%误差线')
    ax2.axhline(y=1, color='orange', linestyle='--', alpha=0.5, label='1%误差线')
    ax2.set_xlabel('光学厚度 τ', fontsize=12)
    ax2.set_ylabel('相对误差 (%)', fontsize=12)
    ax2.set_title('案例三：近似方法误差分析', fontsize=14, fontweight='bold')
    ax2.set_xlim([0.001, 2])
    ax2.legend(fontsize=10)
    ax2.grid(True, alpha=0.3)
    
    # 图3: 不同厚度下的发射率对比
    ax3 = axes[1, 0]
    L_test = np.array([0.001, 0.005, 0.01, 0.05, 0.1, 0.3, 0.5, 1.0])
    tau_test = kappa * L_test
    eps_ex_test = 1 - np.exp(-tau_test)
    eps_ap_test = tau_test
    
    x_pos = np.arange(len(L_test))
    width = 0.35
    ax3.bar(x_pos - width/2, eps_ex_test, width, label='精确解', alpha=0.8)
    ax3.bar(x_pos + width/2, eps_ap_test, width, label='光学薄近似', alpha=0.8)
    ax3.set_xlabel('介质厚度 (m)', fontsize=12)
    ax3.set_ylabel('发射率 ε', fontsize=12)
    ax3.set_title('案例三：不同厚度下的发射率对比', fontsize=14, fontweight='bold')
    ax3.set_xticks(x_pos)
    ax3.set_xticklabels([f'{L:.3f}' for L in L_test], rotation=45)
    ax3.legend(fontsize=11)
    ax3.grid(True, alpha=0.3, axis='y')
    
    # 图4: 适用范围判据
    ax4 = axes[1, 1]
    
    # 创建适用性区域图
    tau_fine = np.linspace(0.001, 1, 100)
    error_fine = np.abs(1 - np.exp(-tau_fine) - tau_fine) / (1 - np.exp(-tau_fine)) * 100
    
    ax4.fill_between(tau_fine, 0, 1, where=(error_fine < 0.5), alpha=0.3, color='green', label='高精度 (误差<0.5%)')
    ax4.fill_between(tau_fine, 0, 1, where=((error_fine >= 0.5) & (error_fine < 5)), alpha=0.3, color='yellow', label='工程精度 (误差<5%)')
    ax4.fill_between(tau_fine, 0, 1, where=(error_fine >= 5), alpha=0.3, color='red', label='低精度 (误差>5%)')
    ax4.plot(tau_fine, error_fine/100, 'b-', linewidth=2, label='相对误差曲线')
    ax4.axhline(y=0.05, color='orange', linestyle='--', linewidth=1.5)
    ax4.axhline(y=0.005, color='green', linestyle='--', linewidth=1.5)
    ax4.set_xlabel('光学厚度 τ', fontsize=12)
    ax4.set_ylabel('相对误差 (归一化)', fontsize=12)
    ax4.set_title('案例三：光学薄近似适用范围', fontsize=14, fontweight='bold')
    ax4.set_xlim([0, 0.5])
    ax4.set_ylim([0, 0.3])
    ax4.legend(fontsize=9, loc='upper left')
    ax4.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('case3_emissivity_calculation.png', dpi=150, bbox_inches='tight')
    print("  图片已保存: case3_emissivity_calculation.png")
    plt.close()
    
    return epsilon_exact, epsilon_approx

# ==============================================================================
# 案例四：多层光学薄介质辐射
# ==============================================================================
print("\n" + "=" * 80)
print("案例四：多层光学薄介质辐射")
print("=" * 80)

def case4_multilayer_medium():
    """
    案例四：多层光学薄介质辐射
    
    计算通过多层光学薄介质的辐射热流
    """
    print("\n[案例四] 多层光学薄介质辐射")
    print("-" * 60)
    
    # 三层介质参数
    layers = [
        {'kappa': 0.01, 'T': 1000.0, 'L': 0.3},  # 第一层
        {'kappa': 0.02, 'T': 1200.0, 'L': 0.4},  # 第二层
        {'kappa': 0.015, 'T': 1100.0, 'L': 0.3}  # 第三层
    ]
    
    # 壁面参数
    Tw1 = 1500.0  # 上壁面温度
    Tw2 = 800.0   # 下壁面温度
    eps_w = 0.8   # 壁面发射率
    
    # 计算各层光学厚度
    total_L = sum([layer['L'] for layer in layers])
    tau_layers = [layer['kappa'] * layer['L'] for layer in layers]
    tau_total = sum(tau_layers)
    
    print(f"总厚度: {total_L:.2f} m")
    print(f"各层光学厚度: {tau_layers}")
    print(f"总光学厚度: {tau_total:.4f}")
    print(f"光学薄条件: {'满足' if tau_total < 0.1 else '不满足'}")
    
    # 构建空间网格
    Nz = 150
    z = np.linspace(0, total_L, Nz)
    
    # 确定各点所在的层
    kappa = np.zeros(Nz)
    T = np.zeros(Nz)
    
    z_boundary = [0]
    for layer in layers:
        z_boundary.append(z_boundary[-1] + layer['L'])
    
    for i in range(Nz):
        if z[i] < z_boundary[1]:
            kappa[i] = layers[0]['kappa']
            T[i] = layers[0]['T']
        elif z[i] < z_boundary[2]:
            kappa[i] = layers[1]['kappa']
            T[i] = layers[1]['T']
        else:
            kappa[i] = layers[2]['kappa']
            T[i] = layers[2]['T']
    
    # 黑体辐射强度
    Ib = sigma * T**4 / np.pi
    
    # 离散坐标求解
    mu_pos = np.array([0.2958759, 0.9082483])
    w_pos = np.array([0.5235988, 0.5235988])
    
    I_plus = np.zeros((Nz, len(mu_pos)))
    I_minus = np.zeros((Nz, len(mu_pos)))
    
    Ib_w1 = sigma * Tw1**4 / np.pi
    Ib_w2 = sigma * Tw2**4 / np.pi
    I_plus[0, :] = eps_w * Ib_w1
    I_minus[-1, :] = eps_w * Ib_w2
    
    # 迭代求解
    max_iter = 1000
    tol = 1e-6
    
    for iteration in range(max_iter):
        I_plus_old = I_plus.copy()
        I_minus_old = I_minus.copy()
        
        # 正向扫描
        for i in range(1, Nz):
            dz = z[i] - z[i-1]
            dtau = 0.5 * (kappa[i] + kappa[i-1]) * dz
            Ib_avg = 0.5 * (Ib[i] + Ib[i-1])
            for m in range(len(mu_pos)):
                I_plus[i, m] = I_plus[i-1, m] * np.exp(-dtau/mu_pos[m]) + \
                               Ib_avg * (1 - np.exp(-dtau/mu_pos[m]))
        
        # 反向扫描
        for i in range(Nz-2, -1, -1):
            dz = z[i+1] - z[i]
            dtau = 0.5 * (kappa[i] + kappa[i+1]) * dz
            Ib_avg = 0.5 * (Ib[i] + Ib[i+1])
            for m in range(len(mu_pos)):
                I_minus[i, m] = I_minus[i+1, m] * np.exp(-dtau/mu_pos[m]) + \
                                Ib_avg * (1 - np.exp(-dtau/mu_pos[m]))
        
        error = np.max(np.abs(I_plus - I_plus_old)) + np.max(np.abs(I_minus - I_minus_old))
        if error < tol:
            print(f"迭代收敛于第 {iteration+1} 次迭代")
            break
    
    # 计算辐射热流
    q_plus = 2 * np.pi * np.sum(I_plus * mu_pos * w_pos, axis=1)
    q_minus = 2 * np.pi * np.sum(I_minus * mu_pos * w_pos, axis=1)
    q_exact = q_plus - q_minus
    
    # 光学薄近似
    # 壁面间直接辐射
    q_wall = sigma * (Tw1**4 - Tw2**4) / (2/eps_w - 1)
    
    # 各层发射贡献
    q_emission = 0
    for i, layer in enumerate(layers):
        q_layer = 4 * sigma * layer['T']**4 * tau_layers[i]
        q_emission += q_layer
        print(f"  第{i+1}层发射贡献: {q_layer:.2f} W/m²")
    
    q_approx = q_wall + q_emission
    
    print(f"\n总辐射热流:")
    print(f"  精确解 (DOM): {q_exact[0]:.2f} W/m²")
    print(f"  光学薄近似: {q_approx:.2f} W/m²")
    print(f"  相对误差: {abs(q_exact[0] - q_approx)/q_exact[0]*100:.2f}%")
    
    # 可视化
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 图1: 多层介质结构
    ax1 = axes[0, 0]
    colors = ['lightblue', 'lightgreen', 'lightyellow']
    labels = ['Layer 1 (T=1000K)', 'Layer 2 (T=1200K)', 'Layer 3 (T=1100K)']
    
    for i, (layer, color, label) in enumerate(zip(layers, colors, labels)):
        rect = Rectangle((z_boundary[i], 0), layer['L'], 1, 
                         facecolor=color, edgecolor='black', linewidth=2, label=label)
        ax1.add_patch(rect)
        ax1.text(z_boundary[i] + layer['L']/2, 0.5, 
                f'$\\kappa$={layer["kappa"]} m⁻¹\n$\\tau$={tau_layers[i]:.3f}',
                ha='center', va='center', fontsize=10, fontweight='bold')
    
    # 壁面
    ax1.axhline(y=1, color='red', linewidth=4, label=f'Wall 1 (T={Tw1}K)')
    ax1.axhline(y=0, color='blue', linewidth=4, label=f'Wall 2 (T={Tw2}K)')
    
    ax1.set_xlim(-0.1, total_L + 0.1)
    ax1.set_ylim(-0.2, 1.3)
    ax1.set_xlabel('位置 z (m)', fontsize=12)
    ax1.set_title('案例四：多层介质结构', fontsize=14, fontweight='bold')
    ax1.legend(fontsize=9, loc='upper right')
    ax1.set_yticks([])
    
    # 图2: 温度和吸收系数分布
    ax2 = axes[0, 1]
    ax2_twin = ax2.twinx()
    
    ax2.step(z, T, 'r-', linewidth=2.5, where='mid', label='温度')
    ax2_twin.step(z, kappa, 'b--', linewidth=2.5, where='mid', label='吸收系数')
    
    ax2.set_xlabel('位置 z (m)', fontsize=12)
    ax2.set_ylabel('温度 (K)', fontsize=12, color='red')
    ax2_twin.set_ylabel('吸收系数 (m⁻¹)', fontsize=12, color='blue')
    ax2.set_title('案例四：温度和吸收系数分布', fontsize=14, fontweight='bold')
    ax2.tick_params(axis='y', labelcolor='red')
    ax2_twin.tick_params(axis='y', labelcolor='blue')
    ax2.grid(True, alpha=0.3)
    
    # 图3: 辐射热流分布
    ax3 = axes[1, 0]
    ax3.plot(z, q_exact, 'b-', linewidth=2.5, label='精确解 (DOM)')
    ax3.axhline(y=q_approx, color='r', linestyle='--', linewidth=2, label='光学薄近似')
    ax3.set_xlabel('位置 z (m)', fontsize=12)
    ax3.set_ylabel('辐射热流 (W/m²)', fontsize=12)
    ax3.set_title('案例四：辐射热流分布', fontsize=14, fontweight='bold')
    ax3.legend(fontsize=11)
    ax3.grid(True, alpha=0.3)
    
    # 图4: 各层贡献分析
    ax4 = axes[1, 1]
    
    contributions = []
    labels_contrib = []
    
    # 壁面贡献
    contributions.append(q_wall)
    labels_contrib.append('Wall-to-Wall')
    
    # 各层贡献
    for i, layer in enumerate(layers):
        q_layer = 4 * sigma * layer['T']**4 * tau_layers[i]
        contributions.append(q_layer)
        labels_contrib.append(f'Layer {i+1}')
    
    colors_pie = ['lightcoral', 'lightblue', 'lightgreen', 'lightyellow']
    explode = (0.05, 0, 0, 0)
    
    ax4.pie(contributions, labels=labels_contrib, autopct='%1.1f%%',
            colors=colors_pie, explode=explode, startangle=90)
    ax4.set_title('案例四：辐射热流贡献分析', fontsize=14, fontweight='bold')
    
    plt.tight_layout()
    plt.savefig('case4_multilayer_medium.png', dpi=150, bbox_inches='tight')
    print("  图片已保存: case4_multilayer_medium.png")
    plt.close()
    
    return q_exact, q_approx

# ==============================================================================
# 案例五：光学薄介质与壁面辐射耦合
# ==============================================================================
print("\n" + "=" * 80)
print("案例五：光学薄介质与壁面辐射耦合")
print("=" * 80)

def case5_cylinder_coupling():
    """
    案例五：光学薄介质与壁面辐射耦合
    
    圆柱形燃烧室中光学薄气体与壁面的辐射耦合
    """
    print("\n[案例五] 光学薄介质与壁面辐射耦合")
    print("-" * 60)
    
    # 几何参数
    D = 2.0   # 直径, m
    R = D / 2  # 半径, m
    L = 4.0   # 长度, m
    
    # 物理参数
    Tg = 1800.0  # 气体温度, K
    Tw = 1200.0  # 壁面温度, K
    kappa = 0.05 # 吸收系数, m^-1
    eps_w = 0.8  # 壁面发射率
    
    # 光学厚度
    tau_R = kappa * R  # 径向光学厚度
    tau_L = kappa * L  # 轴向光学厚度
    
    print(f"圆柱几何: D={D}m, L={L}m")
    print(f"径向光学厚度: {tau_R:.4f}")
    print(f"轴向光学厚度: {tau_L:.4f}")
    print(f"光学薄条件: {'满足' if max(tau_R, tau_L) < 0.3 else '部分满足'}")
    
    # 蒙特卡洛模拟
    N_photons = 50000  # 光子数
    
    print(f"\n蒙特卡洛模拟: {N_photons} 个光子")
    
    # 壁面热流统计
    Nr = 20  # 径向分段
    Nz = 40  # 轴向分段
    
    r_bins = np.linspace(0, R, Nr+1)
    z_bins = np.linspace(0, L, Nz+1)
    
    q_wall = np.zeros((Nr, Nz))  # 壁面热流分布
    
    # 从气体发射光子
    np.random.seed(42)
    
    for _ in range(N_photons):
        # 在圆柱体内随机发射光子
        # 径向分布: r = R * sqrt(rand) 保证均匀分布
        r_emit = R * np.sqrt(np.random.random())
        theta_emit = 2 * np.pi * np.random.random()
        z_emit = L * np.random.random()
        
        # 发射方向 (各向同性)
        phi = np.arccos(2 * np.random.random() - 1)  # 极角
        psi = 2 * np.pi * np.random.random()          # 方位角
        
        # 方向向量
        dx = np.sin(phi) * np.cos(psi)
        dy = np.sin(phi) * np.sin(psi)
        dz_dir = np.cos(phi)
        
        # 光子位置
        x = r_emit * np.cos(theta_emit)
        y = r_emit * np.sin(theta_emit)
        z = z_emit
        
        # 追踪光子直到被吸收或到达壁面
        absorbed = False
        steps = 0
        max_steps = 100
        
        while not absorbed and steps < max_steps:
            steps += 1
            
            # 计算到圆柱壁面的距离
            # 圆柱面: x² + y² = R²
            a = dx**2 + dy**2
            b = 2 * (x*dx + y*dy)
            c = x**2 + y**2 - R**2
            
            if abs(a) > 1e-10:
                discriminant = b**2 - 4*a*c
                if discriminant >= 0:
                    t_cyl = (-b + np.sqrt(discriminant)) / (2*a)
                    if t_cyl < 0:
                        t_cyl = (-b - np.sqrt(discriminant)) / (2*a)
                else:
                    t_cyl = np.inf
            else:
                t_cyl = np.inf
            
            # 计算到端面的距离
            if dz_dir > 0:
                t_end = (L - z) / dz_dir
            elif dz_dir < 0:
                t_end = -z / dz_dir
            else:
                t_end = np.inf
            
            # 到边界的距离
            t_boundary = min(t_cyl, t_end) if min(t_cyl, t_end) > 0 else max(t_cyl, t_end)
            
            # 光子自由程
            if kappa > 0:
                l_free = -np.log(np.random.random()) / kappa
            else:
                l_free = np.inf
            
            # 判断光子命运
            if l_free < t_boundary:
                # 在介质内被吸收
                x += l_free * dx
                y += l_free * dy
                z += l_free * dz_dir
                absorbed = True
            else:
                # 到达壁面
                x += t_boundary * dx
                y += t_boundary * dy
                z += t_boundary * dz_dir
                
                # 确定壁面位置
                if abs(x**2 + y**2 - R**2) < 0.01:  # 圆柱面
                    r_idx = min(int(r_emit / R * Nr), Nr-1)
                    z_idx = min(int(z / L * Nz), Nz-1)
                    
                    # 计算光子能量贡献
                    photon_energy = sigma * Tg**4 / N_photons
                    
                    # 考虑壁面吸收
                    if np.random.random() < eps_w:
                        q_wall[r_idx, z_idx] += photon_energy * eps_w
                
                absorbed = True
    
    # 计算平均热流
    q_avg = np.mean(q_wall)
    q_max = np.max(q_wall)
    
    print(f"\n壁面辐射热流统计:")
    print(f"  平均热流: {q_avg:.2f} W/m²")
    print(f"  最大热流: {q_max:.2f} W/m²")
    
    # 光学薄近似
    # 气体对壁面的辐射加热
    q_approx = 4 * kappa * R * sigma * (Tg**4 - Tw**4)
    print(f"  光学薄近似: {q_approx:.2f} W/m²")
    
    # 可视化
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 图1: 圆柱几何示意图
    ax1 = axes[0, 0]
    
    # 绘制圆柱截面
    theta = np.linspace(0, 2*np.pi, 100)
    x_cyl = R * np.cos(theta)
    y_cyl = R * np.sin(theta)
    ax1.plot(x_cyl, y_cyl, 'b-', linewidth=3, label='壁面')
    ax1.fill(x_cyl, y_cyl, alpha=0.2, color='lightblue')
    
    # 绘制一些光子轨迹
    for _ in range(20):
        r_emit = R * np.sqrt(np.random.random()) * 0.8
        theta_emit = 2 * np.pi * np.random.random()
        x0 = r_emit * np.cos(theta_emit)
        y0 = r_emit * np.sin(theta_emit)
        
        phi = np.arccos(2 * np.random.random() - 1)
        psi = 2 * np.pi * np.random.random()
        
        dx = np.sin(phi) * np.cos(psi)
        dy = np.sin(phi) * np.sin(psi)
        
        # 简化的轨迹
        t = np.linspace(0, 0.5, 20)
        x_traj = x0 + dx * t
        y_traj = y0 + dy * t
        
        # 只绘制在圆柱内的部分
        mask = x_traj**2 + y_traj**2 < R**2
        ax1.plot(x_traj[mask], y_traj[mask], 'r-', alpha=0.3, linewidth=0.5)
    
    ax1.set_xlim(-R*1.2, R*1.2)
    ax1.set_ylim(-R*1.2, R*1.2)
    ax1.set_aspect('equal')
    ax1.set_xlabel('x (m)', fontsize=12)
    ax1.set_ylabel('y (m)', fontsize=12)
    ax1.set_title('案例五：圆柱燃烧室截面与光子轨迹', fontsize=14, fontweight='bold')
    ax1.grid(True, alpha=0.3)
    
    # 图2: 壁面热流分布 (轴向)
    ax2 = axes[0, 1]
    z_centers = (z_bins[:-1] + z_bins[1:]) / 2
    q_axial = np.mean(q_wall, axis=0)
    
    ax2.plot(z_centers, q_axial, 'b-', linewidth=2.5, marker='o', markersize=4)
    ax2.axhline(y=q_approx, color='r', linestyle='--', linewidth=2, label='光学薄近似')
    ax2.set_xlabel('轴向位置 z (m)', fontsize=12)
    ax2.set_ylabel('辐射热流 (W/m²)', fontsize=12)
    ax2.set_title('案例五：壁面辐射热流分布 (轴向平均)', fontsize=14, fontweight='bold')
    ax2.legend(fontsize=11)
    ax2.grid(True, alpha=0.3)
    
    # 图3: 壁面热流分布热力图
    ax3 = axes[1, 0]
    R_grid, Z_grid = np.meshgrid(r_bins[:-1], z_bins[:-1])
    
    im = ax3.pcolormesh(Z_grid, R_grid, q_wall.T, shading='auto', cmap='hot')
    ax3.set_xlabel('轴向位置 z (m)', fontsize=12)
    ax3.set_ylabel('径向位置 r (m)', fontsize=12)
    ax3.set_title('案例五：壁面辐射热流分布热力图', fontsize=14, fontweight='bold')
    cbar = plt.colorbar(im, ax=ax3)
    cbar.set_label('辐射热流 (W/m²)', fontsize=11)
    
    # 图4: 光学厚度与热流关系
    ax4 = axes[1, 1]
    
    # 不同光学厚度下的热流
    tau_test = np.linspace(0.001, 0.5, 50)
    q_mc = []  # 蒙特卡洛结果 (简化)
    q_thin = []  # 光学薄近似
    
    for tau in tau_test:
        # 光学薄近似
        q_t = 4 * tau * sigma * (Tg**4 - Tw**4)
        q_thin.append(q_t)
        
        # 简化的蒙特卡洛估计
        if tau < 0.1:
            q_mc.append(q_t * (1 + 0.1 * tau))  # 小修正
        else:
            q_mc.append(q_t * (1 - np.exp(-tau)) / tau)
    
    ax4.plot(tau_test, q_mc, 'b-', linewidth=2.5, label='蒙特卡洛模拟')
    ax4.plot(tau_test, q_thin, 'r--', linewidth=2, label='光学薄近似')
    ax4.axvline(x=tau_R, color='g', linestyle=':', linewidth=2, 
                label=f'当前工况 (τ={tau_R:.3f})')
    ax4.set_xlabel('光学厚度 τ', fontsize=12)
    ax4.set_ylabel('辐射热流 (W/m²)', fontsize=12)
    ax4.set_title('案例五：光学厚度对辐射热流的影响', fontsize=14, fontweight='bold')
    ax4.legend(fontsize=11)
    ax4.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('case5_cylinder_coupling.png', dpi=150, bbox_inches='tight')
    print("  图片已保存: case5_cylinder_coupling.png")
    plt.close()
    
    return q_wall, q_approx

# ==============================================================================
# 生成GIF动画
# ==============================================================================
def create_gif_animations():
    """
    生成GIF动画展示辐射传播过程
    """
    print("\n" + "=" * 80)
    print("生成GIF动画")
    print("=" * 80)
    
    # 动画1: 辐射强度传播过程
    print("\n[动画1] 辐射强度传播过程")
    
    fig, ax = plt.subplots(figsize=(10, 6))
    
    # 参数
    T1, T2, Tg = 1500, 1000, 1200
    L = 1.0
    kappa = 0.01
    Nz = 50
    z = np.linspace(0, L, Nz)
    
    # 初始化线条
    line_plus, = ax.plot([], [], 'r-', linewidth=2.5, label=r'$I^+$')
    line_minus, = ax.plot([], [], 'b-', linewidth=2.5, label=r'$I^-$')
    line_ib, = ax.plot([], [], 'g--', linewidth=2, label=r'$I_b$')
    
    ax.set_xlim(0, L)
    ax.set_ylim(0, sigma * T1**4 / np.pi * 1.2)
    ax.set_xlabel('位置 z (m)', fontsize=12)
    ax.set_ylabel('辐射强度 (W/(m²·sr))', fontsize=12)
    ax.set_title('辐射强度传播过程', fontsize=14, fontweight='bold')
    ax.legend(fontsize=11)
    ax.grid(True, alpha=0.3)
    
    # 时间步数
    n_frames = 100
    
    def init():
        line_plus.set_data([], [])
        line_minus.set_data([], [])
        line_ib.set_data([], [])
        return line_plus, line_minus, line_ib
    
    def update(frame):
        # 模拟辐射强度随时间建立稳态
        progress = frame / n_frames
        
        # 正向辐射强度 (从z=0传播)
        I_plus = sigma * T1**4 / np.pi * np.exp(-kappa * z * (1 - progress)) + \
                 sigma * Tg**4 / np.pi * (1 - np.exp(-kappa * z * progress))
        
        # 反向辐射强度 (从z=L传播)
        I_minus = sigma * T2**4 / np.pi * np.exp(-kappa * (L - z) * (1 - progress)) + \
                  sigma * Tg**4 / np.pi * (1 - np.exp(-kappa * (L - z) * progress))
        
        # 黑体辐射强度
        Ib = np.full_like(z, sigma * Tg**4 / np.pi)
        
        line_plus.set_data(z, I_plus)
        line_minus.set_data(z, I_minus)
        line_ib.set_data(z, Ib)
        
        ax.set_title(f'辐射强度传播过程 (进度: {progress*100:.0f}%)', 
                    fontsize=14, fontweight='bold')
        
        return line_plus, line_minus, line_ib
    
    anim = animation.FuncAnimation(fig, update, init_func=init,
                                   frames=n_frames, interval=100, blit=True)
    anim.save('animation1_intensity_propagation.gif', writer='pillow', fps=10)
    print("  动画已保存: animation1_intensity_propagation.gif")
    plt.close()
    
    # 动画2: 多层介质辐射过程
    print("\n[动画2] 多层介质辐射过程")
    
    fig, ax = plt.subplots(figsize=(10, 6))
    
    # 三层介质
    layers = [
        {'T': 1000, 'L': 0.3, 'color': 'lightblue'},
        {'T': 1200, 'L': 0.4, 'color': 'lightgreen'},
        {'T': 1100, 'L': 0.3, 'color': 'lightyellow'}
    ]
    
    total_L = sum([l['L'] for l in layers])
    z = np.linspace(0, total_L, 100)
    
    # 构建温度分布
    T = np.zeros_like(z)
    z_bounds = [0]
    for layer in layers:
        z_bounds.append(z_bounds[-1] + layer['L'])
    
    for i in range(len(z)):
        if z[i] < z_bounds[1]:
            T[i] = layers[0]['T']
        elif z[i] < z_bounds[2]:
            T[i] = layers[1]['T']
        else:
            T[i] = layers[2]['T']
    
    line_temp, = ax.plot([], [], 'b-', linewidth=2.5)
    
    # 绘制层边界
    for zb in z_bounds[1:-1]:
        ax.axvline(x=zb, color='k', linestyle='--', linewidth=1.5, alpha=0.5)
    
    # 填充层
    for i, layer in enumerate(layers):
        rect = Rectangle((z_bounds[i], 900), layer['L'], 400, 
                         facecolor=layer['color'], alpha=0.3, zorder=0)
        ax.add_patch(rect)
        ax.text(z_bounds[i] + layer['L']/2, 1300, f'Layer {i+1}\nT={layer["T"]}K',
               ha='center', va='center', fontsize=10)
    
    ax.set_xlim(0, total_L)
    ax.set_ylim(900, 1350)
    ax.set_xlabel('位置 z (m)', fontsize=12)
    ax.set_ylabel('温度 (K)', fontsize=12)
    ax.set_title('多层介质辐射过程', fontsize=14, fontweight='bold')
    ax.grid(True, alpha=0.3)
    
    n_frames = 80
    
    def init2():
        line_temp.set_data([], [])
        return line_temp,
    
    def update2(frame):
        progress = frame / n_frames
        
        # 模拟温度分布的建立过程
        T1, T2 = 1500, 800
        T_animate = T1 + (T2 - T1) * z / total_L
        
        # 叠加层温度
        for i in range(len(z)):
            if z[i] < z_bounds[1]:
                T_animate[i] = T1 * (1 - progress) + layers[0]['T'] * progress
            elif z[i] < z_bounds[2]:
                T_animate[i] = T1 * (1 - progress) + layers[1]['T'] * progress
            else:
                T_animate[i] = T1 * (1 - progress) + layers[2]['T'] * progress
        
        line_temp.set_data(z, T_animate)
        ax.set_title(f'多层介质辐射过程 (进度: {progress*100:.0f}%)', 
                    fontsize=14, fontweight='bold')
        
        return line_temp,
    
    anim2 = animation.FuncAnimation(fig, update2, init_func=init2,
                                    frames=n_frames, interval=100, blit=True)
    anim2.save('animation2_multilayer_radiation.gif', writer='pillow', fps=10)
    print("  动画已保存: animation2_multilayer_radiation.gif")
    plt.close()
    
    # 动画3: 光学厚度对辐射的影响
    print("\n[动画3] 光学厚度对辐射的影响")
    
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    
    # 左图: 透射率
    ax1 = axes[0]
    tau_values = np.linspace(0, 2, 100)
    trans_exact = np.exp(-tau_values)
    trans_approx = 1 - tau_values
    
    line_exact, = ax1.plot([], [], 'b-', linewidth=2.5, label='精确解: $e^{-\\tau}$')
    line_approx, = ax1.plot([], [], 'r--', linewidth=2.5, label='光学薄近似: $1-\\tau$')
    point_exact, = ax1.plot([], [], 'bo', markersize=10)
    point_approx, = ax1.plot([], [], 'ro', markersize=10)
    
    ax1.set_xlim(0, 2)
    ax1.set_ylim(-0.5, 1.2)
    ax1.set_xlabel('光学厚度 τ', fontsize=12)
    ax1.set_ylabel('透射率', fontsize=12)
    ax1.set_title('透射率随光学厚度变化', fontsize=14, fontweight='bold')
    ax1.legend(fontsize=11)
    ax1.grid(True, alpha=0.3)
    ax1.axhline(y=0, color='k', linewidth=0.5)
    
    # 右图: 发射率
    ax2 = axes[1]
    eps_exact = 1 - np.exp(-tau_values)
    eps_approx = tau_values
    
    line_eps_exact, = ax2.plot([], [], 'b-', linewidth=2.5, label='精确解: $1-e^{-\\tau}$')
    line_eps_approx, = ax2.plot([], [], 'r--', linewidth=2.5, label='光学薄近似: $\\tau$')
    point_eps_exact, = ax2.plot([], [], 'bo', markersize=10)
    point_eps_approx, = ax2.plot([], [], 'ro', markersize=10)
    
    ax2.set_xlim(0, 2)
    ax2.set_ylim(-0.2, 1.2)
    ax2.set_xlabel('光学厚度 τ', fontsize=12)
    ax2.set_ylabel('发射率', fontsize=12)
    ax2.set_title('发射率随光学厚度变化', fontsize=14, fontweight='bold')
    ax2.legend(fontsize=11)
    ax2.grid(True, alpha=0.3)
    ax2.axhline(y=0, color='k', linewidth=0.5)
    ax2.axhline(y=1, color='k', linewidth=0.5, linestyle='--', alpha=0.5)
    
    n_frames = 100
    
    def init3():
        line_exact.set_data([], [])
        line_approx.set_data([], [])
        point_exact.set_data([], [])
        point_approx.set_data([], [])
        line_eps_exact.set_data([], [])
        line_eps_approx.set_data([], [])
        point_eps_exact.set_data([], [])
        point_eps_approx.set_data([], [])
        return (line_exact, line_approx, point_exact, point_approx,
                line_eps_exact, line_eps_approx, point_eps_exact, point_eps_approx)
    
    def update3(frame):
        progress = (frame + 1) / n_frames
        tau_max = 2 * progress
        
        tau_anim = tau_values[:int(100 * progress)]
        
        line_exact.set_data(tau_anim, trans_exact[:int(100 * progress)])
        line_approx.set_data(tau_anim, trans_approx[:int(100 * progress)])
        
        line_eps_exact.set_data(tau_anim, eps_exact[:int(100 * progress)])
        line_eps_approx.set_data(tau_anim, eps_approx[:int(100 * progress)])
        
        if len(tau_anim) > 0:
            point_exact.set_data([tau_anim[-1]], [trans_exact[int(100 * progress) - 1]])
            point_approx.set_data([tau_anim[-1]], [trans_approx[int(100 * progress) - 1]])
            point_eps_exact.set_data([tau_anim[-1]], [eps_exact[int(100 * progress) - 1]])
            point_eps_approx.set_data([tau_anim[-1]], [eps_approx[int(100 * progress) - 1]])
        
        fig.suptitle(f'光学厚度对辐射特性的影响 (τ_max = {tau_max:.2f})', 
                    fontsize=16, fontweight='bold')
        
        return (line_exact, line_approx, point_exact, point_approx,
                line_eps_exact, line_eps_approx, point_eps_exact, point_eps_approx)
    
    anim3 = animation.FuncAnimation(fig, update3, init_func=init3,
                                    frames=n_frames, 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" + "=" * 80)
    print("开始运行光学薄介质辐射换热仿真")
    print("=" * 80)
    
    # 运行所有案例
    try:
        # 案例一
        q1_exact, q1_approx = case1_isothermal_plates()
        
        # 案例二
        q2_exact, q2_approx = case2_nonisothermal_medium()
        
        # 案例三
        eps_exact, eps_approx = case3_emissivity_calculation()
        
        # 案例四
        q4_exact, q4_approx = case4_multilayer_medium()
        
        # 案例五
        q5_wall, q5_approx = case5_cylinder_coupling()
        
        # 生成GIF动画
        create_gif_animations()
        
        print("\n" + "=" * 80)
        print("所有仿真案例运行完成!")
        print("=" * 80)
        print("\n生成的文件:")
        print("  - case1_isothermal_plates.png")
        print("  - case2_nonisothermal_medium.png")
        print("  - case3_emissivity_calculation.png")
        print("  - case4_multilayer_medium.png")
        print("  - case5_cylinder_coupling.png")
        print("  - animation1_intensity_propagation.gif")
        print("  - animation2_multilayer_radiation.gif")
        print("  - animation3_optical_thickness_effect.gif")
        
    except Exception as e:
        print(f"\n运行出错: {e}")
        import traceback
        traceback.print_exc()

    print("\n" + "=" * 80)
    print("程序结束")
    print("=" * 80)