主题040：辐射传热模型

目录

1. 引言
2. 热辐射基础理论
3. 辐射传递方程(RTE)
4. 气体辐射特性
5. 颗粒辐射特性
6. 辐射传热数值方法
7. 燃烧中的辐射传热
8. Python仿真实现
9. 工程应用案例
10. 总结与展望

---

1. 引言

1.1 辐射传热在燃烧中的重要性

辐射传热是燃烧过程中最重要的传热方式之一，在高温燃烧系统中往往占据主导地位。与导热和对流不同，辐射传热不需要介质，可以在真空中传播，且与温度的四次方成正比，因此在高温条件下变得极其重要。

燃烧系统中辐射传热的重要性体现在：

1. 能量传递：在锅炉、炉窑等燃烧设备中，辐射传热占总传热量的70-90%
2. 温度分布：辐射影响燃烧室内的温度分布，进而影响燃烧效率和污染物排放
3. 火焰可视化：辐射是火焰发光的原因，也是火焰检测和诊断的基础
4. 热应力：不均匀的辐射热流会导致设备热应力和变形
5. 安全分析：辐射热流是火灾安全分析的关键参数

1.2 辐射传热的主要特点

与导热和对流的区别：

特性	导热	对流	辐射
传播介质	需要物质	需要流体	不需要介质
温度依赖	线性	线性	四次方
方向性	无	有	各向异性
波长特性	无	无	光谱分布

辐射传热的复杂性：

• 非线性：与温度的四次方成正比
• 光谱依赖性：不同波长有不同的辐射特性
• 方向性：辐射强度随方向变化
• 参与性介质：气体和颗粒吸收、发射和散射辐射

1.3 辐射传热模型的应用

工业应用：

• 锅炉和炉窑设计优化
• 燃烧室热负荷计算
• 火焰检测和监控
• 火灾安全评估
• 航天器热控设计

科学研究：

• 火焰结构和温度测量
• 燃烧诊断技术
• 新型燃烧技术开发
• 污染物生成预测

---

2. 热辐射基础理论

2.1 黑体辐射

2.1.1 黑体定义

黑体是理想化的辐射体，具有以下特性：

• 吸收所有入射辐射（吸收率α = 1）
• 在相同温度下发射最大的辐射能量
• 辐射特性仅取决于温度，与材料无关

2.1.2 Planck定律

Planck定律描述了黑体辐射的光谱分布：

$$E_{b\lambda }(T)=\frac{2\pi h{c}^2}{{\lambda }^5}\frac{1}{\exp \left(\frac{hc}{\lambda k_{B}T}\right)-1}$$

其中：

• $E_{b\lambda }$：光谱辐射力 (W/(m²·μm))
• $h$：Planck常数 (6.626×10⁻³⁴ J·s)
• $c$：光速 (2.998×10⁸ m/s)
• $k_B$：Boltzmann常数 (1.381×10⁻²³ J/K)
• $\lambda$：波长 (m)
• $T$：温度 (K)

Wien位移定律：

峰值波长与温度的关系：

KaTeX parse error: Undefined control sequence: \cdotp at position 1: \̲c̲d̲o̲t̲p̲

这表明温度越高，辐射峰值向短波方向移动。

Stefan-Boltzmann定律：

黑体总辐射力：

$$E_b(T)=\sigma T^4$$

其中Stefan-Boltzmann常数：

KaTeX parse error: Undefined control sequence: \cdotp at position 1: \̲c̲d̲o̲t̲p̲

2.1.3 实际表面的辐射特性

发射率(ε)：

实际表面的辐射力与黑体辐射力之比：

$$\epsilon =\frac{E}{E_b}$$

发射率取决于：

• 材料种类
• 表面状况（粗糙度、氧化层）
• 温度
• 波长（光谱发射率）
• 方向（方向发射率）

Kirchhoff定律：

在热平衡条件下：

$$\epsilon =\alpha$$

即发射率等于吸收率。对于灰体（光谱特性与波长无关），此关系在非平衡条件下也近似成立。

灰体近似：

当光谱发射率不随波长变化时，称为灰体：

$${\epsilon }_{\lambda }=\epsilon =\text{const}$$

这在工程计算中常用，可大大简化计算。

2.2 辐射强度与辐射力

2.2.1 辐射强度

辐射强度$I$定义为单位面积、单位立体角、单位时间的辐射能量：

KaTeX parse error: Undefined control sequence: \cdotp at position 1: \̲c̲d̲o̲t̲p̲

对于黑体，辐射强度与方向无关（Lambert表面）：

$$I_b=\frac{\sigma T^4}{\pi }$$

2.2.2 辐射力

半球辐射力：单位表面积向半球空间发射的总辐射功率：

$$E={\int }_{2\pi }I\cos \theta d\Omega$$

对于Lambert表面：

$$E=\pi I$$

光谱辐射力：单位波长间隔内的辐射力：

$$E_{\lambda }=\frac{dE}{d\lambda }$$

2.3 视角因子

2.3.1 定义

视角因子$F_{ij}$表示从表面$i$发射的辐射能直接到达表面$j$的比例：

$$F_{ij}=\frac{1}{A_i}{\int }_{A_i}{\int }_{A_j}\frac{\cos {\theta }_i\cos {\theta }_j}{\pi r^2}d{A}_{j}d{A}_i$$

其中：

• $A_i$, $A_j$：表面面积
• $r$：两微元面之间的距离
• ${\theta }_i$, ${\theta }_j$：连线与法线的夹角

2.3.2 视角因子性质

互易关系：

$$A_i{F}_{ij}=A_j{F}_{ji}$$

完整性：

$$\sum _{j=1}^N{F}_{ij}=1$$

自视角因子：对于凸表面或平面，$F_{ii}=0$

2.3.3 常见几何形状的视角因子

平行平板：

$$F_{12}=F_{21}=1$$

同心球或无限长同心圆柱：

$$F_{12}=1,F_{21}=\frac{A_1}{A_2}$$

两个无限大平行平面：

$$F_{12}=F_{21}=1$$

微元面与矩形面：有专门的图表和公式可查

---

3. 辐射传递方程(RTE)

3.1 方程推导

3.1.1 辐射在介质中的传播

考虑沿方向$s$传播的辐射，辐射强度的变化由以下因素引起：

吸收：介质吸收辐射能量

$$d{I}_{\lambda }=-{\kappa }_{\lambda }{I}_{\lambda }ds$$

其中${\kappa }_{\lambda }$是光谱吸收系数。

发射：介质发射辐射能量

$$d{I}_{\lambda }={\kappa }_{\lambda }{I}_{b\lambda }ds$$

根据Kirchhoff定律，发射系数等于吸收系数。

散射：辐射被散射到其他方向

$$d{I}_{\lambda }=-{\sigma }_{s\lambda }{I}_{\lambda }ds$$

其中${\sigma }_{s\lambda }$是光谱散射系数。

入散射：从其他方向散射来的辐射

$$d{I}_{\lambda }=\frac{{\sigma }_{s\lambda }}{4\pi }{\int }_{4\pi }{I}_{\lambda }({\mathbf{s}}^{\prime })\Phi ({\mathbf{s}}^{\prime },\mathbf{s})d{\Omega }^{\prime }ds$$

其中$\Phi$是散射相函数。

3.1.2 辐射传递方程

综合以上各项，得到辐射传递方程(RTE)：

$$\frac{d{I}_{\lambda }}{ds}=\underset{\text{吸收}}{\underbrace{-{\kappa }_{\lambda }{I}_{\lambda }}}+\underset{\text{发射}}{\underbrace{{\kappa }_{\lambda }{I}_{b\lambda }}}\underset{\text{散射损失}}{\underbrace{-{\sigma }_{s\lambda }{I}_{\lambda }}}+\underset{\text{散射增益}}{\underbrace{\frac{{\sigma }_{s\lambda }}{4\pi }{\int }_{4\pi }{I}_{\lambda }\Phi d{\Omega }^{\prime }}}$$

定义消光系数：

$${\beta }_{\lambda }={\kappa }_{\lambda }+{\sigma }_{s\lambda }$$

和单次散射反照率：

$${\omega }_{\lambda }=\frac{{\sigma }_{s\lambda }}{{\beta }_{\lambda }}$$

RTE可改写为：

$$\frac{d{I}_{\lambda }}{ds}=-{\beta }_{\lambda }{I}_{\lambda }+{\beta }_{\lambda }(1-{\omega }_{\lambda })I_{b\lambda }+\frac{{\beta }_{\lambda }{\omega }_{\lambda }}{4\pi }{\int }_{4\pi }{I}_{\lambda }\Phi d{\Omega }^{\prime }$$

3.1.3 方程的物理意义

RTE描述了辐射强度在参与性介质中的输运过程：

• 左边：辐射强度的空间变化率
• 右边第一项：消光（吸收+散射）导致的衰减
• 右边第二项：介质热发射
• 右边第三项：散射引起的强度变化

3.2 光学厚度与发射率

3.2.1 光学厚度

光学厚度定义为：

$${\tau }_{\lambda }={\int }_0^L{\beta }_{\lambda }ds$$

物理意义：无量纲路径长度，表征介质对辐射的阻挡能力。

分类：

• 光学薄：$\tau \ll 1$，辐射几乎不受阻碍地通过
• 光学厚：$\tau \gg 1$，辐射在短距离内被吸收

3.2.2 介质发射率

对于等温均匀介质层，光谱发射率：

$${\epsilon }_{\lambda }=1-\exp (-{\tau }_{\lambda })$$

总发射率需要对光谱和方向积分。

3.3 边界条件

3.3.1 表面边界条件

在壁面处，出射辐射强度包括：

1. 表面发射的辐射
2. 入射辐射的反射

$$I_{\lambda }^{+}({\mathbf{r}}_w,\mathbf{s})={\epsilon }_{\lambda }({\mathbf{r}}_w)I_{b\lambda }(T_w)+\frac{1-{\epsilon }_{\lambda }({\mathbf{r}}_w)}{\pi }{\int }_{\mathbf{n}\cdot {\mathbf{s}}^{\prime }<0}{I}_{\lambda }^{-}({\mathbf{r}}_w,{\mathbf{s}}^{\prime })|\mathbf{n}\cdot {\mathbf{s}}^{\prime }|d{\Omega }^{\prime }$$

其中：

• $I_{\lambda }^{+}$：出射辐射强度（$\mathbf{n}\cdot \mathbf{s}>0$）
• $I_{\lambda }^{-}$：入射辐射强度（$\mathbf{n}\cdot \mathbf{s}<0$）
• ${\epsilon }_{\lambda }$：表面光谱发射率

3.3.2 对称边界条件

在对称边界上：

$$I(\mathbf{r},\mathbf{s})=I(\mathbf{r},{\mathbf{s}}^{\prime })$$

其中${\mathbf{s}}^{\prime }$是$\mathbf{s}$关于对称面的镜像方向。

3.3.3 入口/出口边界条件

对于开放边界：

• 入口：给定入射辐射强度（通常来自环境）
• 出口：无入射辐射（$I^{-}=0$）

---

4. 气体辐射特性

4.1 气体辐射的特点

与固体和液体不同，气体辐射具有以下特点：

1. 选择性：只在特定波段（谱带）发射和吸收辐射
2. 体积性：辐射来自整个体积，不仅是表面
3. 温度依赖性：辐射特性强烈依赖于温度
4. 压力依赖性：高压下谱线展宽效应显著

4.2 主要辐射性气体

4.2.1 双原子分子

N₂和O₂：

• 对称分子，无永久偶极矩
• 辐射很弱，通常可忽略

CO：

• 在2.35 μm和4.67 μm有强吸收带
• 燃烧产物中的重要辐射气体

4.2.2 三原子分子

CO₂：
主要吸收/发射波段：

• 2.7 μm带（最强）
• 4.3 μm带（很重要）
• 15 μm带（重要）

H₂O：
主要吸收/发射波段：

• 1.38 μm和1.87 μm（较弱）
• 2.7 μm（强，与CO₂重叠）
• 6.3 μm（强）

4.2.3 碳氢化合物

CH₄：

• 3.3 μm和7.6 μm有强吸收带
• 天然气燃烧中的重要辐射气体

其他烃类：

• 在3.4 μm附近有C-H键振动带

4.3 谱带模型

4.3.1 谱线结构

分子辐射/吸收由大量离散的谱线组成，谱线位置由量子力学决定。

谱线展宽机制：

1. 自然展宽：由Heisenberg不确定性原理决定
2. Doppler展宽：由分子热运动引起，高温下显著
3. 碰撞展宽（压力展宽）：由分子间碰撞引起，高压下显著

Voigt线型：综合考虑Doppler和压力展宽：

$${\varphi }_V(\nu )={\int }_{-\infty }^{\infty }{\varphi }_D({\nu }^{\prime }){\varphi }_C(\nu -{\nu }^{\prime })d{\nu }^{\prime }$$

4.3.2 窄谱带模型

当光谱分辨率足够高时，可以分辨单个谱线。常用模型：

Elsasser模型：假设等强度、等间距的谱线

统计模型（Goody模型）：假设谱线位置和强度随机分布

窄谱带平均吸收率：

$$\bar{A}=1-\exp \left(-\frac{\bar{W}}{d}\right)$$

其中$\bar{W}$是平均等效线宽，$d$是平均谱线间距。

4.3.3 宽谱带模型

在宽光谱范围内平均，常用Edwards指数宽谱带模型：

$$\frac{\bar{W}}{d}=\frac{\alpha }{\omega }{\left[\ln \left(\frac{{\tau }_0}{\alpha /\omega }+1\right)\right]}^{1/2}$$

其中：

• $\alpha$：积分吸收率
• $\omega$：谱带宽
• ${\tau }_0$：光学厚度

4.3.4 总发射率模型

Hottel图表：

基于实验数据的总发射率图表，考虑了：

• 气体温度
• 分压力
• 射线长度（平均光束长度）

Leckner模型：

经验关联式，将总发射率表示为：

$${\epsilon }_g=\exp \left(a_0+a_1\ln T+a_2{\ln }^{2}T+\cdots \right)$$

系数$a_i$是压力和射线长度的函数。

4.4 加权求和灰气体模型(WSGGM)

4.4.1 模型原理

将非灰气体近似为多个灰气体的加权和：

$${\epsilon }_g(T_g,T_b)=\sum _{i=0}^N{a}_i(T_g)\left[1-\exp (-{\kappa }_{i}pL)\right]$$

其中：

• $a_i$：权重系数，与气体温度有关
• ${\kappa }_i$：灰气体吸收系数
• $p$：气体分压
• $L$：射线长度

4.4.2 模型参数

对于CO₂-H₂O混合气体，Smith等人给出的四灰气体模型参数：

i	$a_i$	${\kappa }_i$ (atm⁻¹m⁻¹)
0	0.350	0.0 (透明窗口)
1	0.350	0.5
2	0.225	5.0
3	0.075	50.0

权重系数$a_i$随温度变化。

4.4.3 模型优点

1. 计算效率高，适合CFD耦合
2. 精度较好，误差通常在10-20%
3. 可以考虑气体温度与壁温不同的情况

---

5. 颗粒辐射特性

5.1 颗粒辐射的重要性

在燃烧系统中，颗粒（煤粉、炭、灰、碳烟）对辐射有重要影响：

1. 增强辐射：颗粒通常是强辐射体，增强传热
2. 改变光谱特性：颗粒辐射通常是连续的（灰体特性）
3. 散射效应：颗粒散射辐射，改变辐射方向

5.2 单颗粒辐射特性

5.2.1 颗粒尺寸参数

定义无量纲尺寸参数：

$$x=\frac{\pi D}{\lambda }$$

其中$D$是颗粒直径，$\lambda$是辐射波长。

分类：

• $x\ll 1$：Rayleigh散射区（小颗粒）
• $x\approx 1$：Mie散射区（中等颗粒）
• $x\gg 1$：几何光学区（大颗粒）

5.2.2 Mie理论

Mie理论是描述球形颗粒对电磁波散射的精确解。

消光效率：

$$Q_{ext}=Q_{abs}+Q_{sca}$$

其中：

• $Q_{ext}$：消光效率
• $Q_{abs}$：吸收效率
• $Q_{sca}$：散射效率

效率因子与尺寸参数的关系：

对于大颗粒（$x\gg 1$）：

$$Q_{abs}\approx 1,Q_{sca}\approx 1$$

对于小颗粒（$x\ll 1$，Rayleigh区）：

$$Q_{abs}\propto x,Q_{sca}\propto x^4$$

5.2.3 复折射率

颗粒的光学性质由复折射率描述：

$$m=n-i\kappa$$

其中：

• $n$：实部，折射率
• $\kappa$：虚部，吸收指数

复折射率与波长和材料有关。对于碳烟，典型值：

• $n\approx 1.5-2.0$
• $\kappa \approx 0.5-1.0$

5.3 颗粒云辐射特性

5.3.1 体积吸收和散射系数

对于颗粒云，定义体积系数：

体积吸收系数：

$${\kappa }_p={\int }_0^{\infty }N(D)\frac{\pi D^2}{4}{Q}_{abs}(D)dD$$

体积散射系数：

$${\sigma }_{sp}={\int }_0^{\infty }N(D)\frac{\pi D^2}{4}{Q}_{sca}(D)dD$$

其中$N(D)$是颗粒数密度分布。

5.3.2 散射相函数

描述散射辐射的方向分布。常用近似：

各向同性散射：

$$\Phi (\Theta )=1$$

线性各向异性散射：

$$\Phi (\Theta )=1+g\cos \Theta$$

其中$g$是不对称因子：

• $g=0$：各向同性
• $g>0$：前向散射为主
• $g<0$：后向散射为主

Henyey-Greenstein相函数：

$$\Phi (\Theta )=\frac{1-g^2}{(1+g^2-2g\cos \Theta {)}^{3/2}}$$

5.3.3 碳烟辐射

碳烟是燃烧中最重要的颗粒辐射体：

碳烟体积分数：

$$f_v=\frac{\text{碳烟质量}}{\text{气体体积}\times {\rho }_{\text{碳烟}}}$$

Rayleigh近似下的吸收系数：

$${\kappa }_{soot}=\frac{36\pi nk}{(n^2-k^2+2{)}^2+4n^2{k}^2}\frac{f_v}{\lambda }$$

对于典型碳烟（$n=1.85$, $k=0.48$）：

$${\kappa }_{soot}\approx 8.7\frac{f_v}{\lambda }$$

这表明碳烟吸收与波长成反比（与黑体不同）。

---

6. 辐射传热数值方法

6.1 离散坐标法(DOM)

6.1.1 方法原理

将连续的方向空间离散为有限个方向：

$${\int }_{4\pi }f(\mathbf{s})d\Omega \approx \sum _{i=1}^M{w}_{i}f({\mathbf{s}}_i)$$

其中${\mathbf{s}}_i$是离散方向，$w_i$是权重。

6.1.2 Sₙ方法

常用的方向离散方案，如S₂、S₄、S₆、S₈等。

S₂方法（2个方向）：

• 方向：$\mu =\pm 1/\sqrt{3}$
• 权重：$w=2\pi$

S₄方法（12个方向）：

• 更精确，但计算量更大

6.1.3 空间离散

沿每个离散方向求解RTE：

$$\frac{d{I}_i}{ds}=-\beta I_i+\beta S_i$$

其中$S_i$是源项。

使用有限差分或有限体积法离散空间。

6.1.4 伪散射问题

DOM的主要缺点：

• 射线效应：离散方向导致的数值误差
• 假散射：数值扩散

改进方法：

• 增加离散方向数
• 使用有限体积法（FVM）
• 使用其他方法如FVM、MC等

6.2 有限体积法(FVM)

6.2.1 方法原理

与DOM类似，但在控制体积内对方向进行积分：

$${\int }_{\Delta V}{\int }_{\Delta \Omega }\frac{dI}{ds}d\Omega dV=\cdots$$

6.2.2 优点

1. 守恒性好
2. 射线效应较弱
3. 与CFD方法一致，便于耦合

6.2.3 实施步骤

1. 将空间域划分为控制体积
2. 将方向空间划分为立体角
3. 对每个控制体积和立体角，建立离散方程
4. 迭代求解

6.3 蒙特卡洛法(MC)

6.3.1 方法原理

通过模拟大量光子的随机传播来求解辐射传递：

1. 从发射点发射光子
2. 随机决定光子的传播方向
3. 计算自由程，确定吸收或散射位置
4. 如果散射，随机决定新方向
5. 重复直到光子被吸收或离开系统

6.3.2 随机抽样

发射位置抽样：

根据发射功率分布：

$$P(\mathbf{r})=\frac{\kappa (\mathbf{r})I_b(\mathbf{r})}{\int \kappa I_{b}dV}$$

方向抽样：

对于各向同性发射：

$$\theta ={\cos }^{-1}(1-2R_1),\varphi =2\pi R_2$$

其中$R_1$, $R_2$是[0,1]区间均匀分布的随机数。

自由程抽样：

根据Beer定律：

$$s=-\frac{1}{\beta }\ln (1-R)$$

6.3.3 优点和缺点

优点：

• 可以处理复杂几何
• 可以处理各向异性散射
• 可以处理光谱特性
• 并行性好

缺点：

• 统计噪声
• 计算量大
• 收敛慢

6.4 其他方法

6.4.1 球谐函数法(Pₙ)

将辐射强度展开为球谐函数的级数：

$$I(\mathbf{r},\mathbf{s})=\sum _{l=0}^N\sum _{m=-l}^l{I}_l^m(\mathbf{r})Y_l^m(\mathbf{s})$$

对于一维问题，简化为P₁近似：

$$I=\frac{1}{4\pi }(G+3\mathbf{q}\cdot \mathbf{s})$$

其中$G$是入射辐射，$\mathbf{q}$是辐射热流。

6.4.2 扩散近似

光学厚介质的近似：

$$\mathbf{q}=-\frac{1}{3\beta }\nabla G$$

类似于Fourier导热定律。

6.4.3 光学薄近似

光学薄介质中，忽略自吸收：

$$\nabla \cdot \mathbf{q}=4\pi \kappa I_b$$

---

7. 燃烧中的辐射传热

7.1 火焰辐射特性

7.1.1 非发光火焰

• 主要辐射源：CO₂和H₂O气体
• 蓝色（来自CH*化学发光）
• 辐射较弱

7.1.2 发光火焰

• 含有碳烟颗粒
• 黄/橙色
• 辐射强（碳烟是强辐射体）

碳烟辐射的重要性：

• 增强传热
• 影响火焰温度分布
• 影响燃烧效率

7.1.3 煤粉火焰

• 挥发分燃烧：气体辐射为主
• 焦炭燃烧：颗粒辐射为主
• 灰分：液态灰滴辐射

7.2 辐射与能量方程耦合

7.2.1 能量方程

包含辐射源项的能量方程：

$$\rho c_p\frac{DT}{Dt}=\nabla \cdot (k\nabla T)-\nabla \cdot {\mathbf{q}}_r+{\dot{Q}}_{\text{chem}}$$

其中辐射源项：

$$-\nabla \cdot {\mathbf{q}}_r=\kappa (G-4\pi I_b)$$

7.2.2 耦合求解

顺序求解：

1. 求解流动和能量方程（无辐射）
2. 求解RTE得到辐射热流
3. 更新温度场
4. 重复直到收敛

耦合求解：

• 同时求解所有方程
• 更稳定但计算量大

7.3 辐射传热计算实例

7.3.1 一维平行平板

两个无限大平行平板之间的辐射传热：

纯辐射介质：

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

其中$\tau$是光学厚度。

7.3.2 圆柱形燃烧室

考虑气体辐射的圆柱形燃烧室：

平均光束长度：

$$L_m=3.6\frac{V}{A}$$

其中$V$是体积，$A$是表面积。

总发射率法：

使用Hottel图表或关联式计算气体发射率，然后：

$$q={\epsilon }_g\sigma T_g^4-{\alpha }_g\sigma T_w^4$$

---

8. Python仿真实现

8.1 黑体辐射计算

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

# 物理常数
h = 6.626e-34      # Planck常数 (J·s)
c = 2.998e8        # 光速 (m/s)
k_B = 1.381e-23    # Boltzmann常数 (J/K)
sigma = 5.67e-8    # Stefan-Boltzmann常数 (W/(m²·K⁴))

print("=" * 70)
print("主题040：辐射传热模型 - Python仿真")
print("=" * 70)

# ============================================================================
# 实验1: 黑体辐射特性
# ============================================================================
print("\n" + "=" * 70)
print("实验1: 黑体辐射特性")
print("=" * 70)

def planck_law(lambda_m, T):
    """
    Planck黑体辐射定律
    lambda_m: 波长 (m)
    T: 温度 (K)
    返回: 光谱辐射力 (W/(m²·m))
    """
    C1 = 2 * np.pi * h * c**2  # 第一辐射常数
    C2 = h * c / k_B           # 第二辐射常数 (hc/k_B)
    
    E_lambda = C1 / (lambda_m**5 * (np.exp(C2 / (lambda_m * T)) - 1))
    return E_lambda

def wien_displacement(T):
    """
    Wien位移定律
    返回峰值波长 (μm)
    """
    return 2898 / T

# 温度范围
T_values = [1000, 1500, 2000, 2500, 3000]  # K

# 波长范围 (0.1 - 20 μm)
wavelengths = np.linspace(0.1, 20, 1000) * 1e-6  # m

fig, axes = plt.subplots(2, 2, figsize=(12, 10))

colors = plt.cm.jet(np.linspace(0, 1, len(T_values)))

# 绘制不同温度下的黑体辐射光谱
for i, T in enumerate(T_values):
    E_lambda = planck_law(wavelengths, T)
    axes[0, 0].plot(wavelengths*1e6, E_lambda/1e6, color=colors[i], 
                    label=f'T={T}K', linewidth=2)
    
    # 标记峰值
    lambda_max = wien_displacement(T)
    E_max = planck_law(lambda_max*1e-6, T)
    axes[0, 0].plot(lambda_max, E_max/1e6, 'o', color=colors[i], markersize=8)

axes[0, 0].set_xlabel('Wavelength (μm)', fontsize=11)
axes[0, 0].set_ylabel('Spectral Emissive Power (W/(m²·μm))', fontsize=11)
axes[0, 0].set_title('Blackbody Radiation Spectrum', fontsize=12)
axes[0, 0].legend(fontsize=9)
axes[0, 0].grid(True, alpha=0.3)
axes[0, 0].set_xlim([0, 20])

# Wien位移定律验证
T_range = np.linspace(500, 4000, 100)
lambda_max_values = [wien_displacement(T) for T in T_range]

axes[0, 1].plot(T_range, lambda_max_values, 'b-', linewidth=2)
axes[0, 1].set_xlabel('Temperature (K)', fontsize=11)
axes[0, 1].set_ylabel('Peak Wavelength (μm)', fontsize=11)
axes[0, 1].set_title('Wien Displacement Law', fontsize=12)
axes[0, 1].grid(True, alpha=0.3)

# Stefan-Boltzmann定律验证
T_range_sb = np.linspace(300, 3000, 100)
E_total = [sigma * T**4 for T in T_range_sb]

axes[1, 0].plot(T_range_sb, E_total, 'r-', linewidth=2, label='E = σT⁴')
axes[1, 0].set_xlabel('Temperature (K)', fontsize=11)
axes[1, 0].set_ylabel('Total Emissive Power (W/m²)', fontsize=11)
axes[1, 0].set_title('Stefan-Boltzmann Law', fontsize=12)
axes[1, 0].legend(fontsize=10)
axes[1, 0].grid(True, alpha=0.3)

# 温度对辐射的影响
T_demo = np.linspace(300, 2500, 50)
E_demo = [sigma * T**4 for T in T_demo]

axes[1, 1].semilogy(T_demo, E_demo, 'g-', linewidth=2)
axes[1, 1].set_xlabel('Temperature (K)', fontsize=11)
axes[1, 1].set_ylabel('Total Emissive Power (W/m²)', fontsize=11)
axes[1, 1].set_title('Temperature Effect on Radiation (Log Scale)', fontsize=12)
axes[1, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('exp1_blackbody_radiation.png', dpi=150, bbox_inches='tight')
print("\n结果已保存至: exp1_blackbody_radiation.png")

print("\n黑体辐射特性:")
for T in T_values:
    lambda_max = wien_displacement(T)
    E_total = sigma * T**4
    print(f"  T = {T}K: λ_max = {lambda_max:.2f} μm, E_total = {E_total/1e6:.2e} MW/m²")

plt.close()

8.2 气体辐射特性计算

# ============================================================================
# 实验2: 气体辐射特性
# ============================================================================
print("\n" + "=" * 70)
print("实验2: 气体辐射特性")
print("=" * 70)

def co2_emissivity(T, P_CO2, L):
    """
    CO₂发射率简化模型 (Leckner模型简化版)
    T: 温度 (K)
    P_CO2: CO₂分压 (atm)
    L: 射线长度 (m)
    """
    # 简化关联式
    a0 = -0.35
    a1 = 0.75
    a2 = -0.15
    
    PL = P_CO2 * L * 100  # atm·cm
    
    # 温度修正
    T_factor = (T / 1000)**a1
    
    # 发射率
    epsilon = 1 - np.exp(-(a0 + a2 * np.log(PL)) * T_factor * np.sqrt(PL))
    
    return max(0, min(1, epsilon))

def h2o_emissivity(T, P_H2O, L):
    """
    H₂O发射率简化模型
    """
    # 简化关联式
    b0 = 0.5
    b1 = 0.6
    b2 = -0.1
    
    PL = P_H2O * L * 100  # atm·cm
    
    T_factor = (T / 1000)**b1
    
    epsilon = 1 - np.exp(-(b0 + b2 * np.log(PL)) * T_factor * np.sqrt(PL))
    
    return max(0, min(1, epsilon))

def gas_mixture_emissivity(T, P_CO2, P_H2O, L):
    """
    CO₂-H₂O混合气体发射率
    """
    eps_CO2 = co2_emissivity(T, P_CO2, L)
    eps_H2O = h2o_emissivity(T, P_H2O, L)
    
    # 重叠修正 (简化)
    delta_eps = 0.1 * eps_CO2 * eps_H2O
    
    eps_total = eps_CO2 + eps_H2O - delta_eps
    
    return min(1, eps_total)

fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# 温度对发射率的影响
T_range = np.linspace(500, 2500, 100)
P_CO2 = 0.1  # atm
P_H2O = 0.2  # atm
L = 1.0      # m

eps_CO2 = [co2_emissivity(T, P_CO2, L) for T in T_range]
eps_H2O = [h2o_emissivity(T, P_H2O, L) for T in T_range]
eps_total = [gas_mixture_emissivity(T, P_CO2, P_H2O, L) for T in T_range]

axes[0, 0].plot(T_range, eps_CO2, 'b-', linewidth=2, label='CO₂')
axes[0, 0].plot(T_range, eps_H2O, 'r-', linewidth=2, label='H₂O')
axes[0, 0].plot(T_range, eps_total, 'g--', linewidth=2, label='Mixture')
axes[0, 0].set_xlabel('Temperature (K)', fontsize=11)
axes[0, 0].set_ylabel('Emissivity', fontsize=11)
axes[0, 0].set_title('Gas Emissivity vs Temperature', fontsize=12)
axes[0, 0].legend(fontsize=10)
axes[0, 0].grid(True, alpha=0.3)

# 射线长度对发射率的影响
L_range = np.linspace(0.1, 5.0, 100)
T_fixed = 1500  # K

eps_CO2_L = [co2_emissivity(T_fixed, P_CO2, L) for L in L_range]
eps_H2O_L = [h2o_emissivity(T_fixed, P_H2O, L) for L in L_range]

axes[0, 1].plot(L_range, eps_CO2_L, 'b-', linewidth=2, label='CO₂')
axes[0, 1].plot(L_range, eps_H2O_L, 'r-', linewidth=2, label='H₂O')
axes[0, 1].set_xlabel('Beam Length (m)', fontsize=11)
axes[0, 1].set_ylabel('Emissivity', fontsize=11)
axes[0, 1].set_title(f'Gas Emissivity vs Beam Length (T={T_fixed}K)', fontsize=12)
axes[0, 1].legend(fontsize=10)
axes[0, 1].grid(True, alpha=0.3)

# 分压对发射率的影响
P_range = np.linspace(0.01, 0.5, 100)
L_fixed = 1.0  # m

eps_CO2_P = [co2_emissivity(1500, P, L_fixed) for P in P_range]
eps_H2O_P = [h2o_emissivity(1500, P, L_fixed) for P in P_range]

axes[1, 0].plot(P_range, eps_CO2_P, 'b-', linewidth=2, label='CO₂')
axes[1, 0].plot(P_range, eps_H2O_P, 'r-', linewidth=2, label='H₂O')
axes[1, 0].set_xlabel('Partial Pressure (atm)', fontsize=11)
axes[1, 0].set_ylabel('Emissivity', fontsize=11)
axes[1, 0].set_title('Gas Emissivity vs Partial Pressure', fontsize=12)
axes[1, 0].legend(fontsize=10)
axes[1, 0].grid(True, alpha=0.3)

# 燃烧产物发射率
T_comb = np.linspace(1000, 2500, 100)
# 典型燃烧产物：CO₂=0.1atm, H₂O=0.2atm
eps_comb = [gas_mixture_emissivity(T, 0.1, 0.2, 1.0) for T in T_comb]

axes[1, 1].plot(T_comb, eps_comb, 'm-', linewidth=2)
axes[1, 1].set_xlabel('Temperature (K)', fontsize=11)
axes[1, 1].set_ylabel('Emissivity', fontsize=11)
axes[1, 1].set_title('Combustion Products Emissivity', fontsize=12)
axes[1, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('exp2_gas_radiation.png', dpi=150, bbox_inches='tight')
print("\n结果已保存至: exp2_gas_radiation.png")

print("\n气体辐射特性:")
print(f"  CO₂发射率 (T=1500K, P=0.1atm, L=1m): {co2_emissivity(1500, 0.1, 1.0):.3f}")
print(f"  H₂O发射率 (T=1500K, P=0.2atm, L=1m): {h2o_emissivity(1500, 0.2, 1.0):.3f}")
print(f"  混合气体发射率: {gas_mixture_emissivity(1500, 0.1, 0.2, 1.0):.3f}")

plt.close()

8.3 颗粒辐射特性

# ============================================================================
# 实验3: 颗粒辐射特性
# ============================================================================
print("\n" + "=" * 70)
print("实验3: 颗粒辐射特性")
print("=" * 70)

def mie_efficiency(x, m):
    """
    简化Mie效率计算 (小颗粒和大颗粒近似)
    x: 尺寸参数 πD/λ
    m: 复折射率 n - ik
    """
    n = m.real
    k = m.abs()
    
    if x < 0.1:  # Rayleigh区
        # 吸收效率
        Q_abs = 8 * x * (n**2 * k) / ((n**2 - k**2 + 2)**2 + (2*n*k)**2)
        # 散射效率 (Rayleigh散射 ∝ x⁴)
        Q_sca = (8/3) * x**4 * ((n**2 - 1)/(n**2 + 2))**2
    elif x > 10:  # 几何光学区
        Q_abs = 1.0
        Q_sca = 1.0
    else:  # Mie区 - 简化插值
        f = (x - 0.1) / 9.9
        Q_abs_rayleigh = 8 * 0.1 * (n**2 * k) / ((n**2 - k**2 + 2)**2 + (2*n*k)**2)
        Q_abs = Q_abs_rayleigh * (1-f) + 1.0 * f
        Q_sca_rayleigh = (8/3) * 0.1**4 * ((n**2 - 1)/(n**2 + 2))**2
        Q_sca = Q_sca_rayleigh * (1-f) + 1.0 * f
    
    Q_ext = Q_abs + Q_sca
    
    return Q_ext, Q_abs, Q_sca

def soot_absorption_coefficient(fv, lambda_m, n=1.85, k=0.48):
    """
    碳烟吸收系数 (Rayleigh近似)
    fv: 碳烟体积分数
    lambda_m: 波长 (m)
    """
    # Rayleigh近似
    C_abs = 36 * np.pi * n * k / ((n**2 - k**2 + 2)**2 + 4*n**2*k**2)
    kappa = C_abs * fv / lambda_m
    return kappa

fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# Mie效率随尺寸参数变化
x_range = np.logspace(-2, 2, 100)
m_soot = 1.85 - 0.48j  # 碳烟复折射率

Q_ext_list = []
Q_abs_list = []
Q_sca_list = []

for x in x_range:
    Q_ext, Q_abs, Q_sca = mie_efficiency(x, m_soot)
    Q_ext_list.append(Q_ext)
    Q_abs_list.append(Q_abs)
    Q_sca_list.append(Q_sca)

axes[0, 0].loglog(x_range, Q_ext_list, 'b-', linewidth=2, label='Q_ext')
axes[0, 0].loglog(x_range, Q_abs_list, 'r-', linewidth=2, label='Q_abs')
axes[0, 0].loglog(x_range, Q_sca_list, 'g-', linewidth=2, label='Q_sca')
axes[0, 0].axvline(x=0.1, color='k', linestyle='--', alpha=0.5, label='Rayleigh limit')
axes[0, 0].axvline(x=10, color='k', linestyle='--', alpha=0.5, label='Geometric limit')
axes[0, 0].set_xlabel('Size Parameter x = πD/λ', fontsize=11)
axes[0, 0].set_ylabel('Efficiency Factor', fontsize=11)
axes[0, 0].set_title('Mie Efficiency Factors', fontsize=12)
axes[0, 0].legend(fontsize=9)
axes[0, 0].grid(True, alpha=0.3)

# 碳烟吸收系数随波长变化
wavelengths_soot = np.linspace(0.5, 10, 100) * 1e-6  # m
fv_values = [1e-6, 5e-6, 1e-5, 5e-5]  # 碳烟体积分数
colors_soot = plt.cm.plasma(np.linspace(0, 1, len(fv_values)))

for i, fv in enumerate(fv_values):
    kappa_soot = [soot_absorption_coefficient(fv, lam) for lam in wavelengths_soot]
    axes[0, 1].plot(wavelengths_soot*1e6, kappa_soot, color=colors_soot[i], 
                    label=f'fv={fv:.0e}', linewidth=2)

axes[0, 1].set_xlabel('Wavelength (μm)', fontsize=11)
axes[0, 1].set_ylabel('Absorption Coefficient (1/m)', fontsize=11)
axes[0, 1].set_title('Soot Absorption Coefficient', fontsize=12)
axes[0, 1].legend(fontsize=9)
axes[0, 1].grid(True, alpha=0.3)

# 颗粒直径对辐射特性的影响
D_range = np.linspace(0.01, 10, 100) * 1e-6  # m
lambda_fixed = 1e-6  # 1 μm

x_values = np.pi * D_range / lambda_fixed
Q_abs_diam = []

for x in x_values:
    _, Q_abs, _ = mie_efficiency(x, m_soot)
    Q_abs_diam.append(Q_abs)

axes[1, 0].semilogy(D_range*1e6, Q_abs_diam, 'b-', linewidth=2)
axes[1, 0].set_xlabel('Particle Diameter (μm)', fontsize=11)
axes[1, 0].set_ylabel('Absorption Efficiency', fontsize=11)
axes[1, 0].set_title(f'Absorption vs Diameter (λ={lambda_fixed*1e6:.1f}μm)', fontsize=12)
axes[1, 0].grid(True, alpha=0.3)

# 碳烟发射率与温度
T_soot = np.linspace(1000, 2500, 100)
# 假设碳烟层厚度1cm，体积分数1e-6
fv_soot = 1e-6
L_soot = 0.01  # m

# 使用平均波长 (Wien定律)
eps_soot = []
for T in T_soot:
    lambda_avg = wien_displacement(T) * 1e-6  # m
    kappa = soot_absorption_coefficient(fv_soot, lambda_avg)
    tau = kappa * L_soot
    eps = 1 - np.exp(-tau)
    eps_soot.append(eps)

axes[1, 1].plot(T_soot, eps_soot, 'r-', linewidth=2)
axes[1, 1].set_xlabel('Temperature (K)', fontsize=11)
axes[1, 1].set_ylabel('Soot Layer Emissivity', fontsize=11)
axes[1, 1].set_title('Soot Emissivity vs Temperature', fontsize=12)
axes[1, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('exp3_particle_radiation.png', dpi=150, bbox_inches='tight')
print("\n结果已保存至: exp3_particle_radiation.png")

print("\n颗粒辐射特性:")
print(f"  碳烟复折射率: m = {m_soot.real} - {abs(m_soot.imag)}i")
print(f"  小颗粒(x=0.01)吸收效率: {mie_efficiency(0.01, m_soot)[1]:.4f}")
print(f"  大颗粒(x=100)吸收效率: {mie_efficiency(100, m_soot)[1]:.4f}")

plt.close()

8.4 辐射传递方程求解

# ============================================================================
# 实验4: 一维辐射传递方程求解
# ============================================================================
print("\n" + "=" * 70)
print("实验4: 一维辐射传递方程求解")
print("=" * 70)

def solve_rte_1d(N, tau_total, omega, epsilon_wall, T_gas, T_wall_hot, T_wall_cold):
    """
    求解一维平板几何的RTE (DOM方法)
    
    参数:
    N: 空间网格数
    tau_total: 总光学厚度
    omega: 单次散射反照率
    epsilon_wall: 壁面发射率
    T_gas: 气体温度 (K)
    T_wall_hot: 热壁温度 (K)
    T_wall_cold: 冷壁温度 (K)
    
    返回:
    tau: 光学厚度网格
    G: 入射辐射
    q: 辐射热流
    """
    # 空间离散
    tau = np.linspace(0, tau_total, N)
    dtau = tau[1] - tau[0]
    
    # S4方向离散 (简化版，使用2个方向)
    mu = [1/np.sqrt(3), -1/np.sqrt(3)]  # 方向余弦
    w = [np.pi, np.pi]  # 权重
    n_dir = len(mu)
    
    # 黑体辐射力
    E_b_gas = sigma * T_gas**4 / np.pi
    E_b_hot = sigma * T_wall_hot**4 / np.pi
    E_b_cold = sigma * T_wall_cold**4 / np.pi
    
    # 初始化辐射强度
    I = np.zeros((N, n_dir))
    
    # 源项
    S = (1 - omega) * E_b_gas * np.ones(N)
    
    # 迭代求解
    max_iter = 1000
    tol = 1e-6
    
    for iteration in range(max_iter):
        I_old = I.copy()
        
        # 正向传播 (mu > 0)
        for d in range(n_dir):
            if mu[d] > 0:
                # 边界条件 (冷壁)
                I[0, d] = epsilon_wall * E_b_cold + (1 - epsilon_wall) * I[0, 1-d] if n_dir > 1 else epsilon_wall * E_b_cold
                
                for i in range(1, N):
                    # 离散RTE
                    I[i, d] = (I[i-1, d] * (mu[d]/dtau - 0.5) + S[i]) / (mu[d]/dtau + 0.5)
        
        # 反向传播 (mu < 0)
        for d in range(n_dir):
            if mu[d] < 0:
                # 边界条件 (热壁)
                I[N-1, d] = epsilon_wall * E_b_hot + (1 - epsilon_wall) * I[N-1, 1-d] if n_dir > 1 else epsilon_wall * E_b_hot
                
                for i in range(N-2, -1, -1):
                    I[i, d] = (I[i+1, d] * (-mu[d]/dtau - 0.5) + S[i]) / (-mu[d]/dtau + 0.5)
        
        # 检查收敛
        if np.max(np.abs(I - I_old)) < tol:
            break
    
    # 计算入射辐射和热流
    G = np.zeros(N)
    q = np.zeros(N)
    
    for i in range(N):
        for d in range(n_dir):
            G[i] += w[d] * I[i, d]
            q[i] += w[d] * mu[d] * I[i, d]
    
    return tau, G, q

fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# 案例1: 纯吸收介质 (omega = 0)
tau, G1, q1 = solve_rte_1d(100, 5.0, 0.0, 0.8, 1500, 2000, 1000)

axes[0, 0].plot(tau, G1*1e-6, 'b-', linewidth=2, label='Incident Radiation G')
ax_twin = axes[0, 0].twinx()
ax_twin.plot(tau, q1/1e3, 'r--', linewidth=2, label='Radiative Heat Flux q')
axes[0, 0].set_xlabel('Optical Thickness τ', fontsize=11)
axes[0, 0].set_ylabel('G (MW/(m²·sr))', color='b', fontsize=11)
ax_twin.set_ylabel('q (kW/m²)', color='r', fontsize=11)
axes[0, 0].set_title('Pure Absorption (ω=0)', fontsize=12)
axes[0, 0].grid(True, alpha=0.3)

# 案例2: 散射介质 (omega = 0.5)
tau, G2, q2 = solve_rte_1d(100, 5.0, 0.5, 0.8, 1500, 2000, 1000)

axes[0, 1].plot(tau, G2*1e-6, 'b-', linewidth=2, label='Incident Radiation G')
ax_twin2 = axes[0, 1].twinx()
ax_twin2.plot(tau, q2/1e3, 'r--', linewidth=2, label='Radiative Heat Flux q')
axes[0, 1].set_xlabel('Optical Thickness τ', fontsize=11)
axes[0, 1].set_ylabel('G (MW/(m²·sr))', color='b', fontsize=11)
ax_twin2.set_ylabel('q (kW/m²)', color='r', fontsize=11)
axes[0, 1].set_title('Scattering Medium (ω=0.5)', fontsize=12)
axes[0, 1].grid(True, alpha=0.3)

# 光学厚度影响
tau_values = [0.1, 0.5, 1.0, 2.0, 5.0, 10.0]
q_hot_wall = []

for tau_val in tau_values:
    tau_grid, G, q = solve_rte_1d(50, tau_val, 0.0, 0.8, 1500, 2000, 1000)
    q_hot_wall.append(q[-1])  # 热壁热流

axes[1, 0].semilogx(tau_values, np.array(q_hot_wall)/1e3, 'bo-', linewidth=2, markersize=8)
axes[1, 0].set_xlabel('Optical Thickness τ', fontsize=11)
axes[1, 0].set_ylabel('Heat Flux at Hot Wall (kW/m²)', fontsize=11)
axes[1, 0].set_title('Effect of Optical Thickness', fontsize=12)
axes[1, 0].grid(True, alpha=0.3)

# 散射反照率影响
omega_values = np.linspace(0, 0.9, 10)
q_hot_wall_omega = []

for omega in omega_values:
    tau_grid, G, q = solve_rte_1d(50, 5.0, omega, 0.8, 1500, 2000, 1000)
    q_hot_wall_omega.append(q[-1])

axes[1, 1].plot(omega_values, np.array(q_hot_wall_omega)/1e3, 'rs-', linewidth=2, markersize=8)
axes[1, 1].set_xlabel('Single Scattering Albedo ω', fontsize=11)
axes[1, 1].set_ylabel('Heat Flux at Hot Wall (kW/m²)', fontsize=11)
axes[1, 1].set_title('Effect of Scattering', fontsize=12)
axes[1, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('exp4_rte_solution.png', dpi=150, bbox_inches='tight')
print("\n结果已保存至: exp4_rte_solution.png")

print("\nRTE求解结果:")
print(f"  纯吸收介质，热壁热流: {q1[-1]/1e3:.2f} kW/m²")
print(f"  散射介质(ω=0.5)，热壁热流: {q2[-1]/1e3:.2f} kW/m²")
print(f"  散射使热流降低: {(1-q2[-1]/q1[-1])*100:.1f}%")

plt.close()

8.5 视角因子计算

# ============================================================================
# 实验5: 视角因子计算
# ============================================================================
print("\n" + "=" * 70)
print("实验5: 视角因子计算")
print("=" * 70)

def view_factor_parallel_rectangles(W, H, L):
    """
    两个平行矩形之间的视角因子
    W: 宽度
    H: 高度
    L: 间距
    """
    X = W / L
    Y = H / L
    
    # 使用Hottel交叉线法或查表
    # 这里使用简化公式
    F = (2 / (np.pi * X * Y)) * (
        np.log(np.sqrt((1 + X**2) * (1 + Y**2) / (1 + X**2 + Y**2))) +
        X * np.sqrt(1 + Y**2) * np.arctan(X / np.sqrt(1 + Y**2)) +
        Y * np.sqrt(1 + X**2) * np.arctan(Y / np.sqrt(1 + X**2)) -
        X * np.arctan(X) - Y * np.arctan(Y)
    )
    
    return F

def view_factor_perpendicular_rectangles(W, H, L):
    """
    两个垂直矩形之间的视角因子
    """
    X = W / L
    Y = H / L
    
    # 简化计算
    F = (1 / (np.pi * Y)) * (
        Y * np.arctan(1/Y) + X * np.arctan(1/X) - 
        np.sqrt(X**2 + Y**2) * np.arctan(1/np.sqrt(X**2 + Y**2)) +
        0.25 * np.log((1 + X**2) * (1 + Y**2) / (1 + X**2 + Y**2))
    )
    
    return F

fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# 平行平板视角因子
W_range = np.linspace(0.1, 5, 50)
L_fixed = 1.0
H_fixed = 1.0

F_parallel = [view_factor_parallel_rectangles(W, H_fixed, L_fixed) for W in W_range]

axes[0, 0].plot(W_range, F_parallel, 'b-', linewidth=2)
axes[0, 0].set_xlabel('W/L Ratio', fontsize=11)
axes[0, 0].set_ylabel('View Factor F₁₂', fontsize=11)
axes[0, 0].set_title(f'Parallel Rectangles (H/L={H_fixed/L_fixed})', fontsize=12)
axes[0, 0].grid(True, alpha=0.3)

# 间距对视角因子的影响
L_range = np.linspace(0.1, 5, 50)
W_fixed = 1.0
H_fixed = 1.0

F_vs_L = [view_factor_parallel_rectangles(W_fixed, H_fixed, L) for L in L_range]

axes[0, 1].plot(L_range, F_vs_L, 'r-', linewidth=2)
axes[0, 1].set_xlabel('Spacing L (m)', fontsize=11)
axes[0, 1].set_ylabel('View Factor F₁₂', fontsize=11)
axes[0, 1].set_title(f'Effect of Spacing (W=H={W_fixed}m)', fontsize=12)
axes[0, 1].grid(True, alpha=0.3)

# 垂直平板视角因子
F_perp = [view_factor_perpendicular_rectangles(W, H_fixed, L_fixed) for W in W_range]

axes[1, 0].plot(W_range, F_perp, 'g-', linewidth=2)
axes[1, 0].set_xlabel('W/L Ratio', fontsize=11)
axes[1, 0].set_ylabel('View Factor F₁₂', fontsize=11)
axes[1, 0].set_title(f'Perpendicular Rectangles (H/L={H_fixed/L_fixed})', fontsize=12)
axes[1, 0].grid(True, alpha=0.3)

# 不同几何形状的视角因子比较
geometries = ['Parallel\n(Inf)', 'Parallel\n(W=H=L)', 'Perpendicular\n(W=H=L)', 
              'Concentric\nSpheres', 'Concentric\nCylinders']
F_values = [1.0, view_factor_parallel_rectangles(1, 1, 1),
            view_factor_perpendicular_rectangles(1, 1, 1), 1.0, 1.0]

colors_geo = ['blue', 'green', 'orange', 'red', 'purple']
axes[1, 1].bar(range(len(geometries)), F_values, color=colors_geo, alpha=0.7)
axes[1, 1].set_xticks(range(len(geometries)))
axes[1, 1].set_xticklabels(geometries, fontsize=9)
axes[1, 1].set_ylabel('View Factor F₁₂', fontsize=11)
axes[1, 1].set_title('View Factors for Different Geometries', fontsize=12)
axes[1, 1].set_ylim([0, 1.2])
axes[1, 1].grid(True, alpha=0.3, axis='y')

plt.tight_layout()
plt.savefig('exp5_view_factors.png', dpi=150, bbox_inches='tight')
print("\n结果已保存至: exp5_view_factors.png")

print("\n视角因子计算结果:")
print(f"  无限大平行平板: F₁₂ = 1.0")
print(f"  平行矩形(W=H=L): F₁₂ = {view_factor_parallel_rectangles(1, 1, 1):.3f}")
print(f"  垂直矩形(W=H=L): F₁₂ = {view_factor_perpendicular_rectangles(1, 1, 1):.3f}")

plt.close()