主题076_辐射传热极限

kkchenjj

21人浏览 · 2026-04-28 22:00:54

kkchenjj · 2026-04-28 22:00:54 发布

主题076：辐射传热极限

摘要

辐射传热极限是热辐射领域的基础理论问题，涉及远场黑体辐射极限和近场超普朗克辐射极限两个核心概念。本主题深入分析斯蒂芬-玻尔兹曼定律描述的黑体辐射极限，探讨近场辐射中由于光子隧穿效应产生的超普朗克热流，建立纳米间隙辐射传热的理论模型，并通过Python仿真展示不同尺度下的传热极限特性。

关键词

传热极限、黑体极限、超普朗克、纳米间隙、热流密度、光子隧穿、近场辐射、远场辐射

在这里插入图片描述

1. 引言

辐射传热极限研究的是两个物体之间通过热辐射能够传递的最大热流密度。这一极限取决于物体之间的距离、温度、材料特性以及周围介质。理解辐射传热极限对于热管理、能源转换和纳米技术应用具有重要意义。

1.1 研究背景

传统热辐射理论基于远场近似，认为热辐射传热存在由斯蒂芬-玻尔兹曼定律决定的黑体极限。然而，当两个物体之间的距离减小到纳米尺度时，由于近场效应和光子隧穿，热流密度可以远超黑体极限，达到所谓的"超普朗克"状态。这一现象为高效热管理和能量转换提供了新的可能性。

1.2 研究意义

辐射传热极限的研究具有重要的科学价值和工程应用前景：

热管理技术：为电子设备散热、热光伏系统等提供理论指导
能源转换：提高热辐射能量转换效率，开发新型热电转换器件
纳米技术：为纳米尺度热传导和能量操控提供理论基础
空间应用：航天器热控系统设计需要精确理解辐射传热极限

2. 理论基础

2.1 黑体辐射极限

2.1.1 斯蒂芬-玻尔兹曼定律

黑体辐射的总辐射功率由斯蒂芬-玻尔兹曼定律给出：

$\sigma A T^4$

其中：

$P$ 为辐射功率（W）
$σ=5.67×10−8\sigma = 5.67 \times 10^{-8}$ W/(m²·K⁴) 为斯蒂芬-玻尔兹曼常数
$A$ 为表面积（m²）
$T$ 为绝对温度（K）

两个黑体表面之间的净辐射换热为：

$qbb=σ(T14−T24)q_{bb} = \sigma (T_1^4 - T_2^4)$

这就是远场黑体辐射传热的理论极限。

2.1.2 普朗克定律与波谱分布

黑体辐射的光谱辐射出射度由普朗克定律描述：

$Mλ(λ,T)=2πhc2λ51ehc/(λkBT)−1M_{\lambda}(\lambda, T) = \frac{2\pi h c^2}{\lambda^5} \frac{1}{e^{hc/(\lambda k_B T)} - 1}$

其中：

$\times 10^{-34}$ J·s 为普朗克常数
$\times 10^8$ m/s 为光速
$kB=1.381×10−23k_B = 1.381 \times 10^{-23}$ J/K 为玻尔兹曼常数
$λ\lambda$ 为波长（m）

2.1.3 维恩位移定律

黑体辐射的峰值波长由维恩位移定律给出：

$λmax=bT\lambda_{max} = \frac{b}{T}$

其中 $\times 10^{-3}$ m·K 为维恩位移常数。

2.2 近场辐射传热理论

2.2.1 近场效应的物理机制

当两个物体之间的距离 $d$ 远小于热辐射的特征波长 $λth=hc/(kBT)\lambda_{th} = hc/(k_B T)$ 时，近场效应变得显著。在室温下（T = 300 K），特征波长约为 10 μm。

近场辐射传热增强的主要机制包括：

光子隧穿（Photon Tunneling）：倏逝波可以穿过纳米间隙，贡献额外的传热通道
表面等离激元耦合（Surface Plasmon Polariton Coupling）：金属表面的集体电子振荡增强近场传热
表面声子极化激元（Surface Phonon Polariton）：极性介质中的晶格振动增强近场耦合

2.2.2 涨落耗散定理

近场辐射传热的理论基础是涨落耗散定理（Fluctuation-Dissipation Theorem）。根据该定理，热辐射源于介质内部电荷密度的热涨落，这些涨落产生电磁场，通过麦克斯韦方程组传播。

介电张量的虚部与能量耗散相关：

$⟨Ji(r,ω)Jj∗(r′,ω′)⟩=ωε0πεij′′(r,ω)Θ(ω,T)δ(r−r′)δ(ω−ω′)\langle J_i(\mathbf{r}, \omega) J_j^*(\mathbf{r}', \omega') \rangle = \frac{\omega \varepsilon_0}{\pi} \varepsilon''_{ij}(\mathbf{r}, \omega) \Theta(\omega, T) \delta(\mathbf{r} - \mathbf{r}') \delta(\omega - \omega')$

其中：

$J$ 为电流密度
$ε′′\varepsilon''$ 为介电张量的虚部
$Θ(ω,T)=ℏω/(eℏω/kBT−1)\Theta(\omega, T) = \hbar \omega / (e^{\hbar \omega / k_B T} - 1)$ 为平均光子能量

2.2.3 近场热流密度公式

两个半无限大介质之间的近场辐射热流密度可以表示为：

$\int_0^{\infty} \frac{d\omega}{2\pi} [\Theta(\omega, T_1) - \Theta(\omega, T_2)] \int \frac{d^2\mathbf{k}}{(2\pi)^2} \mathcal{T}(\omega, \mathbf{k})$

其中 $T(ω,k)\mathcal{T}(\omega, \mathbf{k})$ 是能量透射系数，对于平行平面几何，可以分解为TE和TM模式：

$T(ω,k)={(1−∣r01TE∣2)(1−∣r02TE∣2)∣1−r01TEr02TEe2ikzd∣2TE模式(1−∣r01TM∣2)(1−∣r02TM∣2)∣1−r01TMr02TMe2ikzd∣2TM模式\mathcal{T}(\omega, \mathbf{k}) = \begin{cases} \frac{(1 - |r_{01}^{TE}|^2)(1 - |r_{02}^{TE}|^2)}{|1 - r_{01}^{TE} r_{02}^{TE} e^{2ik_z d}|^2} & \text{TE模式} \\ \frac{(1 - |r_{01}^{TM}|^2)(1 - |r_{02}^{TM}|^2)}{|1 - r_{01}^{TM} r_{02}^{TM} e^{2ik_z d}|^2} & \text{TM模式} \end{cases}$

其中 $r_{01}$ 和 $r_{02}$ 是菲涅尔反射系数， $kz=ω2/c2−k2k_z = \sqrt{\omega^2/c^2 - k^2}$ 是波矢的垂直分量。

2.3 超普朗克辐射

2.3.1 超普朗克因子的定义

超普朗克因子定义为近场热流密度与黑体极限的比值：

$ηSH=qnearqbb\eta_{SH} = \frac{q_{near}}{q_{bb}}$

当 $ηSH>1\eta_{SH} > 1$ 时，称为超普朗克辐射。在纳米间隙中，超普朗克因子可以达到数百甚至数千。

2.3.2 超普朗克辐射的条件

实现超普朗克辐射需要满足以下条件：

纳米级间隙：间隙距离 $d$ 应小于热辐射的特征波长
合适的材料：支持表面等离激元或声子极化激元的材料（如金、银、SiC、SiO₂等）
温度差：足够的温度梯度驱动热流

2.3.3 材料选择对超普朗克辐射的影响

不同材料的近场辐射特性差异显著：

材料类型	代表材料	主要机制	超普朗克因子
金属	Au, Ag, Al	表面等离激元	10-100
极性介质	SiC, SiO₂, BN	表面声子极化激元	100-1000
双曲材料	hBN, 超材料	双曲色散	1000+

3. 远场辐射传热极限

3.1 黑体辐射极限的推导

3.1.1 从普朗克定律到斯蒂芬-玻尔兹曼定律

对普朗克定律在所有波长上积分，得到总辐射出射度：

$\int_0^{\infty} M_{\lambda}(\lambda, T) d\lambda = \sigma T^4$

其中：

$W/(m2⋅K4)\sigma = \frac{2\pi^5 k_B^4}{15 h^3 c^2} = 5.670374 \times 10^{-8} \text{ W/(m}^2\text{·K}^4\text{)}$

3.1.2 黑体极限的物理意义

黑体辐射极限代表了理想情况下热辐射传热的上限。实际物体的辐射功率总是小于或等于同温度黑体的辐射功率，由发射率 $ε\varepsilon$ 修正：

$Mreal=εσT4,0≤ε≤1M_{real} = \varepsilon \sigma T^4, \quad 0 \leq \varepsilon \leq 1$

3.2 灰体辐射极限

3.2.1 灰体假设

灰体假设认为材料的发射率与波长无关。在此假设下，两个灰体表面之间的辐射换热为：

$\frac{\sigma (T_1^4 - T_2^4)}{\frac{1}{\varepsilon_1} + \frac{1}{\varepsilon_2} - 1}$

3.2.2 发射率对传热极限的影响

发射率对辐射传热极限有显著影响：

当两个表面都是理想黑体（ $ε1=ε2=1\varepsilon_1 = \varepsilon_2 = 1$ ）时，达到黑体极限
当发射率降低时，传热极限相应降低
选择性表面（波长相关发射率）可以实现对特定波段的辐射控制

3.3 视角因子对传热极限的影响

3.3.1 有限几何尺寸的修正

对于有限尺寸的表面，需要考虑视角因子 $F_{12}$ ：

$F_{12} \sigma (T_1^4 - T_2^4)$

视角因子取决于表面的几何形状和相对位置，取值范围为 0 到 1。

3.3.2 平行平板间的辐射换热

两个无限大平行平板间的视角因子为 1，辐射换热达到最大。对于有限尺寸的平行平板，视角因子小于 1，传热极限相应降低。

4. 近场辐射传热极限

4.1 纳米间隙辐射传热

4.1.1 间隙距离的影响

近场辐射热流密度强烈依赖于间隙距离 $d$ 。当 $d$ 减小时：

倏逝波贡献增加：更多高波矢的倏逝波能够隧穿间隙
共振耦合增强：表面等离激元或声子极化激元的耦合增强
热流密度增大：可以远超黑体极限

间隙距离与热流密度的关系通常近似为：

$\propto \frac{1}{d^2}$

在极小的间隙（< 10 nm）下，这种关系可能转变为 $1/ d$ 或其他形式。

4.1.2 温度对近场传热的影响

温度对近场辐射传热的影响体现在：

热辐射谱的变化：温度升高，辐射峰值向短波方向移动
材料光学特性的温度依赖性：介电函数随温度变化
非线性效应：高温下可能出现非线性热辐射现象

4.2 表面等离激元增强传热

4.2.1 金属-介质界面的表面等离激元

在金属-介质界面，自由电子的集体振荡形成表面等离激元（Surface Plasmon Polariton, SPP）。SPP的色散关系为：

$kspp=ωcεmεdεm+εdk_{spp} = \frac{\omega}{c} \sqrt{\frac{\varepsilon_m \varepsilon_d}{\varepsilon_m + \varepsilon_d}}$

其中 $εm\varepsilon_m$ 和 $εd\varepsilon_d$ 分别是金属和介质的介电常数。

4.2.2 表面等离激元耦合

当两个金属表面靠近时，各自的SPP发生耦合，形成对称和反对称模式。这种耦合可以显著增强近场传热。

4.3 表面声子极化激元增强传热

4.3.1 极性介质中的声子极化激元

在极性介质（如SiC、SiO₂）中，晶格振动与电磁场耦合形成表面声子极化激元（Surface Phonon Polariton, SPhP）。SPhP在特定的频率范围内（Reststrahlen带）具有强耦合特性。

4.3.2 SiC的近场辐射特性

碳化硅（SiC）是研究近场辐射传热的典型材料。在 10-13 μm 波长范围内，SiC支持强SPhP，可以实现极高的超普朗克因子。

5. 超普朗克辐射的理论分析

5.1 超普朗克因子的计算

5.1.1 理论模型

对于两个半无限大SiC平板，超普朗克因子可以通过数值积分计算：

$ηSH=∫0∞dω[Θ(ω,T1)−Θ(ω,T2)]∫0∞dkkT(ω,k)σ(T14−T24)\eta_{SH} = \frac{\int_0^{\infty} d\omega [\Theta(\omega, T_1) - \Theta(\omega, T_2)] \int_0^{\infty} dk k \mathcal{T}(\omega, k)}{\sigma (T_1^4 - T_2^4)}$

其中波矢积分包括传播波（ $\omega/c$ ）和倏逝波（ $\omega/c$ ）的贡献。

5.1.2 数值计算方法

计算超普朗克因子需要：

材料的介电函数：通常使用Drude-Lorentz模型或实验数据
菲涅尔反射系数：根据介电函数计算
双重数值积分：对频率和波矢进行积分

5.2 超普朗克辐射的物理限制

5.2.1 能量守恒约束

虽然近场辐射可以超过黑体极限，但仍受能量守恒约束。超普朗克辐射并不意味着违反热力学定律，而是反映了近场条件下额外的传热通道。

5.2.2 量子极限

在极小的间隙（< 1 nm）下，量子效应变得重要。量子极限下的近场传热需要考虑：

非局域效应：介电响应的非局域性
量子隧穿：电子隧穿对传热的贡献
卡西米尔效应：量子真空涨落的影响

5.3 双曲材料中的超普朗克辐射

5.3.1 双曲色散特性

双曲材料（如六方氮化硼hBN）在特定频率范围内具有双曲型色散关系：

$kx2+ky2ε⊥+kz2ε∥=ω2c2\frac{k_x^2 + k_y^2}{\varepsilon_{\perp}} + \frac{k_z^2}{\varepsilon_{\parallel}} = \frac{\omega^2}{c^2}$

当 $ε⊥\varepsilon_{\perp}$ 和 $ε∥\varepsilon_{\parallel}$ 符号相反时，等频面为双曲面，支持高波矢模式。

5.3.2 双曲材料的高超普朗克因子

双曲材料可以支持无限多的高波矢模式，理论上可以实现任意大的超普朗克因子。实际中，材料损耗和几何因素限制了最大热流密度。

6. 工程应用

6.1 近场热光伏系统

6.1.1 系统原理

近场热光伏（Near-Field Thermophotovoltaic, NF-TPV）系统利用近场辐射增强热辐射与光伏电池的耦合，提高能量转换效率。

系统效率可以表示为：

$η=PoutPin=qnear⋅ηPVqnear+qloss\eta = \frac{P_{out}}{P_{in}} = \frac{q_{near} \cdot \eta_{PV}}{q_{near} + q_{loss}}$

其中 $ηPV\eta_{PV}$ 是光伏电池的效率， $q_{loss}$ 是热损失。

6.1.2 性能优化

优化NF-TPV系统的关键参数：

间隙距离：减小间隙以增强近场传热
发射器温度：提高温度以增加辐射功率
材料选择：选择支持强近场耦合的材料
光伏电池带隙：匹配热辐射谱与电池响应

6.2 纳米尺度热管理

6.2.1 纳米器件散热

纳米电子器件的散热是一个重大挑战。近场辐射可以提供额外的散热通道，帮助解决热管理问题。

6.2.2 热整流与热晶体管

利用近场辐射的非对称性，可以实现：

热整流器：单向热流控制
热晶体管：热流的放大和开关控制
热逻辑器件：基于热流的逻辑运算

6.3 热辐射调控

6.3.1 选择性热辐射器

通过纳米结构设计，可以实现对热辐射谱的精确调控：

窄带辐射器：特定波段的强辐射
宽带吸收器：全波段的强吸收
角度选择性辐射器：特定方向的辐射增强

6.3.2 热伪装与隐身

利用近场辐射特性，可以实现：

热伪装：改变物体的热辐射特征
热隐身：抑制特定方向的热辐射
热错觉：制造虚假的热信号

7. 仿真案例分析

7.1 案例1：黑体辐射极限验证

7.1.1 仿真目标

验证斯蒂芬-玻尔兹曼定律，计算不同温度下黑体辐射的总功率和光谱分布。

7.1.2 理论计算

根据普朗克定律，计算光谱辐射出射度：

import numpy as np

def planck_law(wavelength, T):
    """普朗克定律计算光谱辐射出射度"""
    h = 6.626e-34  # 普朗克常数
    c = 3e8        # 光速
    k_B = 1.381e-23  # 玻尔兹曼常数
    
    M_lambda = (2 * np.pi * h * c**2) / (wavelength**5) / \
               (np.exp(h * c / (wavelength * k_B * T)) - 1)
    return M_lambda

7.1.3 仿真结果分析

通过数值积分验证斯蒂芬-玻尔兹曼常数，分析不同温度下的辐射谱分布和峰值波长变化。

7.2 案例2：间隙距离对近场传热的影响

7.2.1 仿真目标

研究两个SiC平板之间的近场辐射传热随间隙距离的变化规律。

7.2.2 理论模型

使用近场辐射传热公式，考虑SiC的介电函数：

def sic_dielectric(omega):
    """SiC的介电函数（Drude-Lorentz模型）"""
    omega_LO = 1.827e14  # 纵向光学声子频率
    omega_TO = 1.495e14  # 横向光学声子频率
    Gamma = 8.97e11      # 阻尼系数
    eps_inf = 6.7        # 高频介电常数
    
    epsilon = eps_inf * (omega_LO**2 - omega**2 - 1j*Gamma*omega) / \
              (omega_TO**2 - omega**2 - 1j*Gamma*omega)
    return epsilon

7.2.3 仿真结果分析

计算不同间隙距离（1 nm - 10 μm）下的热流密度，分析超普朗克因子随距离的变化规律。

7.3 案例3：材料选择对超普朗克辐射的影响

7.3.1 仿真目标

比较不同材料（金属Au、极性介质SiC、双曲材料hBN）的近场辐射特性。

7.3.2 理论模型

建立不同材料的介电函数模型，计算各自的近场热流密度。

7.3.3 仿真结果分析

对比不同材料的超普朗克因子，分析材料特性对近场传热的影响。

7.4 案例4：温度对近场传热的影响

7.4.1 仿真目标

研究温度差对近场辐射传热的影响，分析热流密度与温度的关系。

7.4.2 理论模型

考虑温度对介电函数的影响，计算不同温度组合下的热流密度。

7.4.3 仿真结果分析

分析热流密度随温度差的变化，验证近场传热的非线性特性。

7.5 案例5：近场热光伏系统性能分析

7.5.1 仿真目标

设计并优化近场热光伏系统，分析系统效率与关键参数的关系。

7.5.2 理论模型

建立NF-TPV系统的能量平衡模型，计算系统效率。

7.5.3 仿真结果分析

优化发射器温度、间隙距离和光伏电池带隙，获得最大系统效率。

7.6 案例6：超普朗克辐射动态演示

7.6.1 仿真目标

创建动画展示超普朗克辐射随间隙距离变化的动态过程。

7.6.2 可视化设计

设计动态图表，展示热流密度、超普朗克因子和辐射谱随间隙距离的变化。

7.6.3 动画效果

生成GIF动画，直观展示近场辐射的增强效应。

附录：Python仿真代码

完整的Python仿真代码见同目录下的 run_simulation.py 文件。代码包含以下主要功能：

黑体辐射极限计算：验证斯蒂芬-玻尔兹曼定律
近场传热仿真：计算不同间隙距离下的热流密度
材料特性分析：比较不同材料的近场辐射特性
温度影响分析：研究温度对近场传热的影响
热光伏系统优化：设计和优化近场热光伏系统
动态可视化：生成超普朗克辐射的动态演示动画

运行代码将生成6个案例的仿真结果和可视化图表。

"""
主题076：辐射传热极限仿真程序
=====================================
本程序包含以下仿真案例：
1. 黑体辐射极限验证
2. 间隙距离对近场传热的影响
3. 材料选择对超普朗克辐射的影响
4. 温度对近场传热的影响
5. 近场热光伏系统性能分析
6. 超普朗克辐射动态演示动画
"""

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle, FancyBboxPatch, Circle, FancyArrowPatch
from matplotlib.collections import LineCollection
import matplotlib.animation as animation
from scipy import integrate
from scipy.special import exp1
import warnings
warnings.filterwarnings('ignore')

# 物理常数
h = 6.626e-34  # 普朗克常数 (J·s)
c = 2.998e8    # 光速 (m/s)
k_B = 1.381e-23  # 玻尔兹曼常数 (J/K)
sigma = 5.670374e-8  # 斯蒂芬-玻尔兹曼常数 (W/(m²·K⁴))


def planck_law(wavelength, T):
    """
    普朗克定律：计算黑体光谱辐射出射度
    
    参数:
        wavelength: 波长 (m)
        T: 温度 (K)
    
    返回:
        M_lambda: 光谱辐射出射度 (W/(m²·m))
    """
    M_lambda = (2 * np.pi * h * c**2) / (wavelength**5) / \
               (np.exp(h * c / (wavelength * k_B * T)) - 1)
    return M_lambda


def wien_displacement(T):
    """
    维恩位移定律：计算峰值波长
    
    参数:
        T: 温度 (K)
    
    返回:
        lambda_max: 峰值波长 (m)
    """
    b = 2.898e-3  # 维恩位移常数 (m·K)
    return b / T


def stefan_boltzmann(T):
    """
    斯蒂芬-玻尔兹曼定律：计算黑体总辐射功率密度
    
    参数:
        T: 温度 (K)
    
    返回:
        M: 辐射出射度 (W/m²)
    """
    return sigma * T**4


def sic_dielectric(omega):
    """
    SiC的介电函数（Drude-Lorentz模型）
    
    参数:
        omega: 角频率 (rad/s)
    
    返回:
        epsilon: 复介电常数
    """
    omega_LO = 1.827e14  # 纵向光学声子频率 (rad/s)
    omega_TO = 1.495e14  # 横向光学声子频率 (rad/s)
    Gamma = 8.97e11      # 阻尼系数 (rad/s)
    eps_inf = 6.7        # 高频介电常数
    
    epsilon = eps_inf * (omega_LO**2 - omega**2 - 1j*Gamma*omega) / \
              (omega_TO**2 - omega**2 - 1j*Gamma*omega)
    return epsilon


def au_dielectric(omega):
    """
    金的介电函数（Drude模型）
    
    参数:
        omega: 角频率 (rad/s)
    
    返回:
        epsilon: 复介电常数
    """
    omega_p = 1.37e16  # 等离子体频率 (rad/s)
    gamma = 4.05e13    # 阻尼系数 (rad/s)
    
    epsilon = 1 - omega_p**2 / (omega**2 + 1j*gamma*omega)
    return epsilon


def sio2_dielectric(omega):
    """
    SiO2的介电函数（简化模型）
    
    参数:
        omega: 角频率 (rad/s)
    
    返回:
        epsilon: 复介电常数
    """
    # 简化模型参数
    omega_TO = 8.9e13   # 横向光学声子频率
    omega_LO = 1.03e14  # 纵向光学声子频率
    Gamma = 5e12        # 阻尼系数
    eps_inf = 2.1       # 高频介电常数
    
    epsilon = eps_inf * (omega_LO**2 - omega**2 - 1j*Gamma*omega) / \
              (omega_TO**2 - omega**2 - 1j*Gamma*omega)
    return epsilon


def fresnel_coefficients(k, omega, epsilon1, epsilon2, d):
    """
    计算菲涅尔反射系数和近场透射系数
    
    参数:
        k: 平行波矢 (1/m)
        omega: 角频率 (rad/s)
        epsilon1, epsilon2: 两种材料的介电常数
        d: 间隙距离 (m)
    
    返回:
        T_te, T_tm: TE和TM模式的透射系数
    """
    k0 = omega / c
    kz0 = np.sqrt(k0**2 - k**2 + 0j)
    kz1 = np.sqrt(k0**2 * epsilon1 - k**2 + 0j)
    kz2 = np.sqrt(k0**2 * epsilon2 - k**2 + 0j)
    
    # TE模式菲涅尔反射系数
    r01_te = (kz0 - kz1) / (kz0 + kz1)
    r02_te = (kz0 - kz2) / (kz0 + kz2)
    
    # TM模式菲涅尔反射系数
    r01_tm = (epsilon1 * kz0 - kz1) / (epsilon1 * kz0 + kz1)
    r02_tm = (epsilon2 * kz0 - kz2) / (epsilon2 * kz0 + kz2)
    
    # 透射系数
    if np.isreal(kz0):
        # 传播波
        T_te = 0
        T_tm = 0
    else:
        # 倏逝波
        T_te = np.abs((1 - np.abs(r01_te)**2) * (1 - np.abs(r02_te)**2) / \
                      np.abs(1 - r01_te * r02_te * np.exp(2j * kz0 * d))**2)
        T_tm = np.abs((1 - np.abs(r01_tm)**2) * (1 - np.abs(r02_tm)**2) / \
                      np.abs(1 - r01_tm * r02_tm * np.exp(2j * kz0 * d))**2)
    
    return T_te, T_tm


def near_field_heat_flux(d, T1, T2, material='SiC'):
    """
    计算近场辐射热流密度（简化模型）
    
    参数:
        d: 间隙距离 (m)
        T1, T2: 两个表面的温度 (K)
        material: 材料类型 ('SiC', 'Au', 'SiO2')
    
    返回:
        q: 热流密度 (W/m²)
        eta_sh: 超普朗克因子
    """
    # 黑体极限
    q_bb = sigma * (T1**4 - T2**4)
    
    # 根据材料选择介电函数
    if material == 'SiC':
        eps_func = sic_dielectric
        enhancement_factor = 500  # SiC的典型增强因子
    elif material == 'Au':
        eps_func = au_dielectric
        enhancement_factor = 50   # 金属的典型增强因子
    elif material == 'SiO2':
        eps_func = sio2_dielectric
        enhancement_factor = 200  # SiO2的典型增强因子
    else:
        enhancement_factor = 1
    
    # 简化的近场增强模型
    # 实际计算需要复杂的数值积分
    d_nm = d * 1e9  # 转换为nm
    
    if d_nm < 10:
        # 极小间隙：强近场效应
        near_field_factor = enhancement_factor * (10 / (d_nm + 1))**2
    elif d_nm < 100:
        # 小间隙：中等近场效应
        near_field_factor = enhancement_factor * (100 / (d_nm + 10))**1.5
    elif d_nm < 1000:
        # 亚微米间隙：弱近场效应
        near_field_factor = 1 + (enhancement_factor - 1) * (1000 - d_nm) / 900
    else:
        # 远场：无近场效应
        near_field_factor = 1
    
    q_near = q_bb * near_field_factor
    eta_sh = near_field_factor
    
    return q_near, eta_sh


def case1_blackbody_limit():
    """案例1：黑体辐射极限验证"""
    print("\n" + "=" * 70)
    print("案例1：黑体辐射极限验证")
    print("=" * 70)
    
    # 温度范围
    temperatures = np.linspace(273, 1273, 100)  # 0°C 到 1000°C
    
    # 计算辐射功率
    radiation_power = [stefan_boltzmann(T) for T in temperatures]
    
    # 验证斯蒂芬-玻尔兹曼常数
    T_test = 300  # K
    M_numerical = integrate.quad(lambda lam: planck_law(lam, T_test), 
                                  1e-8, 1e-3)[0]
    M_theoretical = stefan_boltzmann(T_test)
    
    print(f"\n斯蒂芬-玻尔兹曼常数验证:")
    print(f"  数值积分结果: {M_numerical:.6e} W/m²")
    print(f"  理论公式结果: {M_theoretical:.6e} W/m²")
    print(f"  相对误差: {abs(M_numerical - M_theoretical) / M_theoretical * 100:.4f}%")
    
    # 不同温度下的辐射特性
    print(f"\n不同温度下的黑体辐射特性:")
    print("-" * 60)
    print(f"{'温度(K)':<12} {'温度(°C)':<12} {'辐射功率(W/m²)':<20} {'峰值波长(μm)':<15}")
    print("-" * 60)
    
    for T in [300, 500, 700, 1000, 1273]:
        M = stefan_boltzmann(T)
        lambda_max = wien_displacement(T) * 1e6  # 转换为μm
        print(f"{T:<12} {T-273.15:<12.1f} {M:<20.2e} {lambda_max:<15.2f}")
    
    # 绘图
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 图1：斯蒂芬-玻尔兹曼定律
    ax1 = axes[0, 0]
    ax1.plot(temperatures, radiation_power, 'b-', linewidth=2, label='Blackbody')
    ax1.set_xlabel('Temperature (K)', fontsize=11)
    ax1.set_ylabel('Radiation Power (W/m²)', fontsize=11)
    ax1.set_title('Stefan-Boltzmann Law', fontsize=12, fontweight='bold')
    ax1.grid(True, alpha=0.3)
    ax1.legend(fontsize=10)
    
    # 添加T^4拟合
    T_fit = np.linspace(273, 1273, 100)
    power_fit = sigma * T_fit**4
    ax1.plot(T_fit, power_fit, 'r--', linewidth=1.5, alpha=0.7, label=r'$\sigma T^4$')
    ax1.legend(fontsize=10)
    
    # 图2：普朗克定律（不同温度）
    ax2 = axes[0, 1]
    wavelengths = np.linspace(0.1, 20, 500) * 1e-6  # 0.1 - 20 μm
    colors = plt.cm.jet(np.linspace(0, 1, 5))
    
    for i, T in enumerate([300, 500, 700, 1000, 1273]):
        spectral_power = [planck_law(lam, T) for lam in wavelengths]
        ax2.plot(wavelengths * 1e6, np.array(spectral_power) / 1e9, 
                 color=colors[i], linewidth=2, label=f'T = {T} K')
        
        # 标记峰值
        lambda_max = wien_displacement(T) * 1e6
        peak_idx = np.argmin(np.abs(wavelengths * 1e6 - lambda_max))
        ax2.plot(lambda_max, spectral_power[peak_idx] / 1e9, 'o', 
                 color=colors[i], markersize=8)
    
    ax2.set_xlabel('Wavelength (μm)', fontsize=11)
    ax2.set_ylabel('Spectral Radiance (GW/(m²·μm))', fontsize=11)
    ax2.set_title('Planck Law at Different Temperatures', fontsize=12, fontweight='bold')
    ax2.legend(fontsize=9)
    ax2.grid(True, alpha=0.3)
    ax2.set_xlim(0, 20)
    
    # 图3：维恩位移定律验证
    ax3 = axes[1, 0]
    T_range = np.linspace(200, 2000, 100)
    lambda_max_range = [wien_displacement(T) * 1e6 for T in T_range]  # μm
    
    ax3.plot(T_range, lambda_max_range, 'g-', linewidth=2)
    ax3.set_xlabel('Temperature (K)', fontsize=11)
    ax3.set_ylabel('Peak Wavelength (μm)', fontsize=11)
    ax3.set_title('Wien Displacement Law', fontsize=12, fontweight='bold')
    ax3.grid(True, alpha=0.3)
    
    # 添加理论曲线
    b = 2.898e3  # μm·K
    ax3.plot(T_range, b / T_range, 'r--', linewidth=1.5, alpha=0.7, label=r'$\lambda_{max} = b/T$')
    ax3.legend(fontsize=10)
    
    # 图4：温度与辐射功率关系（对数坐标）
    ax4 = axes[1, 1]
    ax4.loglog(temperatures, radiation_power, 'b-', linewidth=2, marker='o', 
               markersize=4, markevery=10)
    ax4.set_xlabel('Temperature (K)', fontsize=11)
    ax4.set_ylabel('Radiation Power (W/m²)', fontsize=11)
    ax4.set_title('Radiation Power vs Temperature (Log-Log)', fontsize=12, fontweight='bold')
    ax4.grid(True, alpha=0.3, which='both')
    
    # 添加斜率标注
    ax4.text(400, 1e4, r'$Slope = 4$', fontsize=12, 
             bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))
    
    plt.tight_layout()
    plt.savefig('case1_blackbody_limit.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print("\n结果已保存至: case1_blackbody_limit.png")
    print("\n案例1完成！")


def case2_gap_distance_effect():
    """案例2：间隙距离对近场传热的影响"""
    print("\n" + "=" * 70)
    print("案例2：间隙距离对近场传热的影响")
    print("=" * 70)
    
    # 温度设置
    T1 = 700  # K
    T2 = 300  # K
    
    # 间隙距离范围
    gaps = np.logspace(-9, -5, 100)  # 1 nm 到 100 μm
    
    # 计算不同材料的近场热流
    materials = ['SiC', 'Au', 'SiO2']
    colors = ['blue', 'gold', 'cyan']
    
    results = {}
    
    for material in materials:
        q_values = []
        eta_values = []
        
        for d in gaps:
            q, eta = near_field_heat_flux(d, T1, T2, material)
            q_values.append(q)
            eta_values.append(eta)
        
        results[material] = {
            'q': np.array(q_values),
            'eta': np.array(eta_values)
        }
    
    # 黑体极限
    q_bb = sigma * (T1**4 - T2**4)
    
    print(f"\n仿真参数:")
    print(f"  高温表面温度: {T1} K ({T1-273.15:.1f}°C)")
    print(f"  低温表面温度: {T2} K ({T2-273.15:.1f}°C)")
    print(f"  黑体极限热流: {q_bb:.2e} W/m²")
    
    print(f"\n不同间隙距离下的近场热流密度 (SiC材料):")
    print("-" * 60)
    print(f"{'间隙距离':<15} {'热流密度(W/m²)':<20} {'超普朗克因子':<15}")
    print("-" * 60)
    
    for d_val in [1e-9, 10e-9, 100e-9, 1e-6, 10e-6, 100e-6]:
        idx = np.argmin(np.abs(gaps - d_val))
        q_val = results['SiC']['q'][idx]
        eta_val = results['SiC']['eta'][idx]
        print(f"{d_val*1e9:.1f} nm{'':<8} {q_val:<20.2e} {eta_val:<15.1f}")
    
    # 绘图
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 图1：热流密度vs间隙距离
    ax1 = axes[0, 0]
    for i, material in enumerate(materials):
        ax1.loglog(gaps * 1e9, results[material]['q'], 
                   color=colors[i], linewidth=2, label=material)
    
    ax1.axhline(y=q_bb, color='red', linestyle='--', linewidth=2, 
                label='Blackbody Limit')
    ax1.set_xlabel('Gap Distance (nm)', fontsize=11)
    ax1.set_ylabel('Heat Flux (W/m²)', fontsize=11)
    ax1.set_title('Near-Field Heat Flux vs Gap Distance', fontsize=12, fontweight='bold')
    ax1.legend(fontsize=10)
    ax1.grid(True, alpha=0.3, which='both')
    
    # 图2：超普朗克因子vs间隙距离
    ax2 = axes[0, 1]
    for i, material in enumerate(materials):
        ax2.semilogx(gaps * 1e9, results[material]['eta'], 
                     color=colors[i], linewidth=2, label=material)
    
    ax2.axhline(y=1, color='red', linestyle='--', linewidth=2, 
                label='Blackbody Limit (η = 1)')
    ax2.set_xlabel('Gap Distance (nm)', fontsize=11)
    ax2.set_ylabel('Super-Planck Factor', fontsize=11)
    ax2.set_title('Super-Planck Factor vs Gap Distance', fontsize=12, fontweight='bold')
    ax2.legend(fontsize=10)
    ax2.grid(True, alpha=0.3)
    
    # 图3：不同间隙下的热流对比（柱状图）
    ax3 = axes[1, 0]
    gap_labels = ['1 nm', '10 nm', '100 nm', '1 μm', '10 μm']
    gap_indices = [0, 20, 40, 60, 80]
    x = np.arange(len(gap_labels))
    width = 0.25
    
    for i, material in enumerate(materials):
        q_vals = [results[material]['q'][idx] / q_bb for idx in gap_indices]
        ax3.bar(x + i * width, q_vals, width, label=material, 
                color=colors[i], alpha=0.7, edgecolor='black')
    
    ax3.set_xlabel('Gap Distance', fontsize=11)
    ax3.set_ylabel('Heat Flux / Blackbody Limit', fontsize=11)
    ax3.set_title('Heat Flux Enhancement at Different Gaps', fontsize=12, fontweight='bold')
    ax3.set_xticks(x + width)
    ax3.set_xticklabels(gap_labels)
    ax3.legend(fontsize=10)
    ax3.grid(True, alpha=0.3, axis='y')
    ax3.axhline(y=1, color='red', linestyle='--', linewidth=1.5, alpha=0.7)
    
    # 图4：近场效应物理机制示意图
    ax4 = axes[1, 1]
    ax4.set_xlim(0, 10)
    ax4.set_ylim(0, 10)
    ax4.axis('off')
    ax4.set_title('Near-Field Radiation Mechanisms', fontsize=12, fontweight='bold')
    
    # 绘制近场效应示意图
    # 表面1
    rect1 = Rectangle((1, 6), 3, 0.5, facecolor='lightblue', edgecolor='blue', linewidth=2)
    ax4.add_patch(rect1)
    ax4.text(2.5, 6.75, 'Surface 1 (T₁)', ha='center', fontsize=10, fontweight='bold')
    
    # 表面2
    rect2 = Rectangle((1, 3.5), 3, 0.5, facecolor='lightcoral', edgecolor='red', linewidth=2)
    ax4.add_patch(rect2)
    ax4.text(2.5, 3.25, 'Surface 2 (T₂)', ha='center', fontsize=10, fontweight='bold')
    
    # 间隙标注
    ax4.annotate('', xy=(4.5, 6), xytext=(4.5, 4),
                 arrowprops=dict(arrowstyle='<->', color='green', lw=2))
    ax4.text(5, 5, 'Gap d', fontsize=10, color='green', fontweight='bold')
    
    # 倏逝波
    for i in range(5):
        x_pos = 1.5 + i * 0.5
        # 从表面1到表面2的波
        ax4.plot([x_pos, x_pos], [6, 4], 'purple', linewidth=1.5, alpha=0.6)
        ax4.plot(x_pos, 5, 'o', color='purple', markersize=4)
    
    ax4.text(2.5, 5, 'Evanescent\nWaves', ha='center', fontsize=9, 
             color='purple', style='italic')
    
    # 说明文字
    ax4.text(5, 8, 'Photon Tunneling Effect', fontsize=11, fontweight='bold',
             bbox=dict(boxstyle='round', facecolor='yellow', alpha=0.3))
    
    mechanisms = [
        '• Evanescent wave coupling',
        '• Surface plasmon polaritons',
        '• Surface phonon polaritons',
        '• Enhanced heat transfer at nanoscale'
    ]
    
    for i, mech in enumerate(mechanisms):
        ax4.text(0.5, 2.5 - i * 0.4, mech, fontsize=9)
    
    plt.tight_layout()
    plt.savefig('case2_gap_distance_effect.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print("\n结果已保存至: case2_gap_distance_effect.png")
    print("\n案例2完成！")


def case3_material_comparison():
    """案例3：材料选择对超普朗克辐射的影响"""
    print("\n" + "=" * 70)
    print("案例3：材料选择对超普朗克辐射的影响")
    print("=" * 70)
    
    # 温度设置
    T1 = 700  # K
    T2 = 300  # K
    
    # 固定间隙距离
    d = 10e-9  # 10 nm
    
    # 材料参数
    materials = {
        'SiC': {
            'type': 'Polar Dielectric',
            'mechanism': 'Surface Phonon Polaritons',
            'color': 'blue',
            'typical_eta': 500
        },
        'SiO2': {
            'type': 'Polar Dielectric',
            'mechanism': 'Surface Phonon Polaritons',
            'color': 'cyan',
            'typical_eta': 200
        },
        'Au': {
            'type': 'Metal',
            'mechanism': 'Surface Plasmon Polaritons',
            'color': 'gold',
            'typical_eta': 50
        },
        'Ag': {
            'type': 'Metal',
            'mechanism': 'Surface Plasmon Polaritons',
            'color': 'silver',
            'typical_eta': 60
        },
        'hBN': {
            'type': 'Hyperbolic Material',
            'mechanism': 'Hyperbolic Phonon Polaritons',
            'color': 'purple',
            'typical_eta': 1000
        }
    }
    
    # 计算各材料的近场热流
    print(f"\n仿真参数:")
    print(f"  间隙距离: {d*1e9:.0f} nm")
    print(f"  高温表面温度: {T1} K")
    print(f"  低温表面温度: {T2} K")
    
    print(f"\n不同材料的近场辐射特性:")
    print("-" * 80)
    print(f"{'Material':<10} {'Type':<20} {'Mechanism':<30} {'η_SH':<10}")
    print("-" * 80)
    
    results = {}
    for material, props in materials.items():
        q, eta = near_field_heat_flux(d, T1, T2, material)
        results[material] = {
            'q': q,
            'eta': eta,
            **props
        }
        print(f"{material:<10} {props['type']:<20} {props['mechanism']:<30} {eta:<10.0f}")
    
    # 绘图
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 图1：超普朗克因子对比（柱状图）
    ax1 = axes[0, 0]
    material_names = list(materials.keys())
    eta_values = [results[m]['eta'] for m in material_names]
    colors_list = [materials[m]['color'] for m in material_names]
    
    bars = ax1.bar(material_names, eta_values, color=colors_list, 
                   alpha=0.7, edgecolor='black', linewidth=1.5)
    ax1.axhline(y=1, color='red', linestyle='--', linewidth=2, 
                label='Blackbody Limit')
    ax1.set_ylabel('Super-Planck Factor', fontsize=11)
    ax1.set_title('Super-Planck Factor by Material', fontsize=12, fontweight='bold')
    ax1.legend(fontsize=10)
    ax1.grid(True, alpha=0.3, axis='y')
    
    # 添加数值标签
    for bar, eta in zip(bars, eta_values):
        height = bar.get_height()
        ax1.text(bar.get_x() + bar.get_width()/2., height,
                 f'{eta:.0f}', ha='center', va='bottom', fontsize=9, fontweight='bold')
    
    # 图2：热流密度对比
    ax2 = axes[0, 1]
    q_bb = sigma * (T1**4 - T2**4)
    q_values = [results[m]['q'] / q_bb for m in material_names]
    
    bars2 = ax2.bar(material_names, q_values, color=colors_list, 
                    alpha=0.7, edgecolor='black', linewidth=1.5)
    ax2.axhline(y=1, color='red', linestyle='--', linewidth=2, 
                label='Blackbody Limit')
    ax2.set_ylabel('Heat Flux / Blackbody Limit', fontsize=11)
    ax2.set_title('Heat Flux Enhancement by Material', fontsize=12, fontweight='bold')
    ax2.legend(fontsize=10)
    ax2.grid(True, alpha=0.3, axis='y')
    
    # 添加数值标签
    for bar, q in zip(bars2, q_values):
        height = bar.get_height()
        ax2.text(bar.get_x() + bar.get_width()/2., height,
                 f'{q:.0f}×', ha='center', va='bottom', fontsize=9, fontweight='bold')
    
    # 图3：材料分类与特性
    ax3 = axes[1, 0]
    ax3.axis('off')
    ax3.set_title('Material Classification for Near-Field Radiation', 
                  fontsize=12, fontweight='bold')
    
    # 创建材料分类表
    table_data = []
    for material, props in materials.items():
        table_data.append([material, props['type'], props['mechanism'], f"{props['typical_eta']:.0f}"])
    
    table = ax3.table(cellText=table_data,
                      colLabels=['Material', 'Type', 'Mechanism', 'Typical η_SH'],
                      cellLoc='center',
                      loc='center',
                      bbox=[0.1, 0.3, 0.8, 0.6])
    table.auto_set_font_size(False)
    table.set_fontsize(9)
    table.scale(1, 2)
    
    # 设置表头样式
    for i in range(4):
        table[(0, i)].set_facecolor('#4CAF50')
        table[(0, i)].set_text_props(weight='bold', color='white')
    
    # 图4：物理机制示意图
    ax4 = axes[1, 1]
    ax4.set_xlim(0, 10)
    ax4.set_ylim(0, 10)
    ax4.axis('off')
    ax4.set_title('Near-Field Enhancement Mechanisms', fontsize=12, fontweight='bold')
    
    # 绘制三种机制的示意图
    mechanisms = [
        ('Surface Plasmon\nPolaritons (SPP)', 'Metal', 8, 'gold'),
        ('Surface Phonon\nPolaritons (SPhP)', 'SiC/SiO₂', 5, 'lightblue'),
        ('Hyperbolic Phonon\nPolaritons (HPP)', 'hBN', 2, 'plum')
    ]
    
    for i, (mechanism, material, y_pos, color) in enumerate(mechanisms):
        # 绘制材料块
        rect = FancyBboxPatch((1, y_pos - 0.5), 2, 1, 
                              boxstyle="round,pad=0.1", 
                              facecolor=color, edgecolor='black', linewidth=2)
        ax4.add_patch(rect)
        ax4.text(2, y_pos, material, ha='center', va='center', 
                fontsize=9, fontweight='bold')
        
        # 绘制说明
        ax4.text(4, y_pos, mechanism, ha='left', va='center', fontsize=9)
        
        # 绘制增强因子指示
        eta_vals = [50, 200, 1000]
        ax4.text(7.5, y_pos, f'η_SH ≈ {eta_vals[i]}', ha='center', va='center',
                fontsize=9, fontweight='bold',
                bbox=dict(boxstyle='round', facecolor='yellow', alpha=0.3))
    
    plt.tight_layout()
    plt.savefig('case3_material_comparison.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print("\n结果已保存至: case3_material_comparison.png")
    print("\n案例3完成！")


def case4_temperature_effect():
    """案例4：温度对近场传热的影响"""
    print("\n" + "=" * 70)
    print("案例4：温度对近场传热的影响")
    print("=" * 70)
    
    # 固定间隙距离和材料
    d = 10e-9  # 10 nm
    material = 'SiC'
    
    # 温度范围
    T2 = 300  # K (固定低温)
    T1_range = np.linspace(350, 1000, 50)  # 高温范围
    
    # 计算不同温度下的热流
    q_values = []
    eta_values = []
    q_bb_values = []
    
    for T1 in T1_range:
        q, eta = near_field_heat_flux(d, T1, T2, material)
        q_bb = sigma * (T1**4 - T2**4)
        
        q_values.append(q)
        eta_values.append(eta)
        q_bb_values.append(q_bb)
    
    q_values = np.array(q_values)
    eta_values = np.array(eta_values)
    q_bb_values = np.array(q_bb_values)
    
    print(f"\n仿真参数:")
    print(f"  间隙距离: {d*1e9:.0f} nm")
    print(f"  材料: {material}")
    print(f"  低温表面温度: {T2} K")
    
    print(f"\n不同温度差下的近场热流:")
    print("-" * 60)
    print(f"{'T₁ (K)':<12} {'ΔT (K)':<12} {'q_near (W/m²)':<20} {'η_SH':<10}")
    print("-" * 60)
    
    for i in [0, 12, 25, 37, 49]:
        T1 = T1_range[i]
        print(f"{T1:<12.0f} {T1-T2:<12.0f} {q_values[i]:<20.2e} {eta_values[i]:<10.0f}")
    
    # 绘图
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 图1：热流密度vs温度
    ax1 = axes[0, 0]
    ax1.plot(T1_range, q_values / 1e6, 'b-', linewidth=2, label='Near-Field')
    ax1.plot(T1_range, q_bb_values / 1e6, 'r--', linewidth=2, label='Blackbody Limit')
    ax1.set_xlabel('Temperature T₁ (K)', fontsize=11)
    ax1.set_ylabel('Heat Flux (MW/m²)', fontsize=11)
    ax1.set_title('Heat Flux vs Temperature', fontsize=12, fontweight='bold')
    ax1.legend(fontsize=10)
    ax1.grid(True, alpha=0.3)
    
    # 图2：超普朗克因子vs温度
    ax2 = axes[0, 1]
    ax2.plot(T1_range, eta_values, 'g-', linewidth=2)
    ax2.axhline(y=1, color='red', linestyle='--', linewidth=2, 
                label='Blackbody Limit')
    ax2.set_xlabel('Temperature T₁ (K)', fontsize=11)
    ax2.set_ylabel('Super-Planck Factor', fontsize=11)
    ax2.set_title('Super-Planck Factor vs Temperature', fontsize=12, fontweight='bold')
    ax2.legend(fontsize=10)
    ax2.grid(True, alpha=0.3)
    
    # 图3：热流增强比vs温度差
    ax3 = axes[1, 0]
    delta_T = T1_range - T2
    enhancement_ratio = q_values / q_bb_values
    
    ax3.plot(delta_T, enhancement_ratio, 'purple', linewidth=2, marker='o', 
             markersize=4, markevery=5)
    ax3.axhline(y=1, color='red', linestyle='--', linewidth=2, 
                label='Blackbody Limit')
    ax3.set_xlabel('Temperature Difference ΔT (K)', fontsize=11)
    ax3.set_ylabel('Heat Flux Enhancement Ratio', fontsize=11)
    ax3.set_title('Heat Flux Enhancement vs Temperature Difference', 
                  fontsize=12, fontweight='bold')
    ax3.legend(fontsize=10)
    ax3.grid(True, alpha=0.3)
    
    # 图4：温度对辐射谱的影响
    ax4 = axes[1, 1]
    wavelengths = np.linspace(1, 20, 200) * 1e-6  # 1 - 20 μm
    
    colors = plt.cm.jet(np.linspace(0, 1, 4))
    for i, T1 in enumerate([400, 600, 800, 1000]):
        spectral_power = [planck_law(lam, T1) for lam in wavelengths]
        ax4.plot(wavelengths * 1e6, np.array(spectral_power) / 1e9, 
                 color=colors[i], linewidth=2, label=f'T₁ = {T1} K')
    
    # 标记SiC的SPhP波长范围
    ax4.axvspan(10, 13, alpha=0.2, color='blue', label='SiC SPhP Band')
    
    ax4.set_xlabel('Wavelength (μm)', fontsize=11)
    ax4.set_ylabel('Spectral Radiance (GW/(m²·μm))', fontsize=11)
    ax4.set_title('Radiation Spectrum at Different Temperatures', 
                  fontsize=12, fontweight='bold')
    ax4.legend(fontsize=9)
    ax4.grid(True, alpha=0.3)
    ax4.set_xlim(1, 20)
    
    plt.tight_layout()
    plt.savefig('case4_temperature_effect.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print("\n结果已保存至: case4_temperature_effect.png")
    print("\n案例4完成！")


def case5_near_field_tpv():
    """案例5：近场热光伏系统性能分析"""
    print("\n" + "=" * 70)
    print("案例5：近场热光伏系统性能分析")
    print("=" * 70)
    
    # 系统参数
    T_emitter = 1500  # 发射器温度 (K)
    T_cell = 300      # 电池温度 (K)
    
    # 光伏电池带隙范围 (eV)
    bandgaps = np.linspace(0.3, 2.0, 50)
    
    # 间隙距离范围
    gaps = [10e-9, 50e-9, 100e-9, 500e-9, 1e-6]  # m
    gap_labels = ['10 nm', '50 nm', '100 nm', '500 nm', '1 μm']
    
    # 计算系统效率
    def calculate_tpv_efficiency(Eg, T_emitter, T_cell, d):
        """计算近场热光伏系统效率"""
        # 简化的效率模型
        # 实际计算需要详细的辐射传递计算
        
        # 近场热流
        q_near, _ = near_field_heat_flux(d, T_emitter, T_cell, 'SiC')
        
        # 光伏电池效率（简化Shockley-Queisser极限）
        # 带隙对应的波长
        h_ev = 4.136e-15  # eV·s
        c = 3e8
        lambda_g = h_ev * c / Eg * 1e6  # μm
        
        # 简化效率公式
        voc_ratio = 0.8  # 开路电压比
        fill_factor = 0.85  # 填充因子
        
        # 带隙效率（高于带隙的光子才能产生电子-空穴对）
        eta_qe = 0.9  # 量子效率
        
        # 电压因子
        v_oc = voc_ratio * Eg
        eta_voltage = v_oc / Eg
        
        # 总体电池效率
        eta_pv = eta_qe * eta_voltage * fill_factor
        
        # 系统效率 = 近场耦合效率 × 电池效率
        # 考虑热损失
        q_loss = 0.1 * q_near  # 10% 热损失
        eta_system = q_near * eta_pv / (q_near + q_loss)
        
        return eta_system * 100  # 转换为百分比
    
    # 计算不同间隙和带隙下的效率
    results = {}
    for i, d in enumerate(gaps):
        efficiencies = []
        for Eg in bandgaps:
            eta = calculate_tpv_efficiency(Eg, T_emitter, T_cell, d)
            efficiencies.append(eta)
        results[gap_labels[i]] = efficiencies
    
    print(f"\n近场热光伏系统参数:")
    print(f"  发射器温度: {T_emitter} K")
    print(f"  电池温度: {T_cell} K")
    print(f"  发射器材料: SiC")
    
    print(f"\n不同间隙距离下的最优效率:")
    print("-" * 50)
    print(f"{'间隙距离':<15} {'最优带隙(eV)':<15} {'最大效率(%)':<15}")
    print("-" * 50)
    
    for label in gap_labels:
        efficiencies = results[label]
        max_eta = max(efficiencies)
        optimal_Eg = bandgaps[efficiencies.index(max_eta)]
        print(f"{label:<15} {optimal_Eg:<15.2f} {max_eta:<15.2f}")
    
    # 绘图
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 图1：效率vs带隙（不同间隙）
    ax1 = axes[0, 0]
    colors = plt.cm.viridis(np.linspace(0, 1, len(gaps)))
    
    for i, label in enumerate(gap_labels):
        ax1.plot(bandgaps, results[label], color=colors[i], linewidth=2, 
                 label=f'Gap = {label}')
    
    # 远场极限（Shockley-Queisser极限约33%）
    ax1.axhline(y=33, color='red', linestyle='--', linewidth=2, 
                label='Shockley-Queisser Limit')
    
    ax1.set_xlabel('Bandgap Energy (eV)', fontsize=11)
    ax1.set_ylabel('System Efficiency (%)', fontsize=11)
    ax1.set_title('TPV Efficiency vs Bandgap', fontsize=12, fontweight='bold')
    ax1.legend(fontsize=9)
    ax1.grid(True, alpha=0.3)
    ax1.set_ylim(0, 40)
    
    # 图2：最大效率vs间隙距离
    ax2 = axes[0, 1]
    max_efficiencies = [max(results[label]) for label in gap_labels]
    gap_values = [10, 50, 100, 500, 1000]  # nm
    
    ax2.semilogx(gap_values, max_efficiencies, 'bo-', linewidth=2, 
                 markersize=8, markerfacecolor='red')
    ax2.axhline(y=33, color='red', linestyle='--', linewidth=2, 
                label='Shockley-Queisser Limit')
    ax2.set_xlabel('Gap Distance (nm)', fontsize=11)
    ax2.set_ylabel('Maximum Efficiency (%)', fontsize=11)
    ax2.set_title('Maximum TPV Efficiency vs Gap Distance', 
                  fontsize=12, fontweight='bold')
    ax2.legend(fontsize=10)
    ax2.grid(True, alpha=0.3)
    
    # 图3：系统示意图
    ax3 = axes[1, 0]
    ax3.set_xlim(0, 10)
    ax3.set_ylim(0, 10)
    ax3.axis('off')
    ax3.set_title('Near-Field TPV System Schematic', fontsize=12, fontweight='bold')
    
    # 发射器
    emitter = Rectangle((1, 6), 3, 1.5, facecolor='red', 
                        edgecolor='darkred', linewidth=2, alpha=0.7)
    ax3.add_patch(emitter)
    ax3.text(2.5, 6.75, f'Emitter\n{T_emitter} K', ha='center', va='center',
            fontsize=10, fontweight='bold', color='white')
    
    # 光伏电池
    cell = Rectangle((1, 2.5), 3, 1, facecolor='blue', 
                     edgecolor='darkblue', linewidth=2, alpha=0.7)
    ax3.add_patch(cell)
    ax3.text(2.5, 3, f'PV Cell\n{T_cell} K', ha='center', va='center',
            fontsize=10, fontweight='bold', color='white')
    
    # 间隙
    ax3.annotate('', xy=(4.5, 6), xytext=(4.5, 3.5),
                 arrowprops=dict(arrowstyle='<->', color='green', lw=2))
    ax3.text(5.2, 4.75, 'Gap d\n(10-100 nm)', fontsize=9, color='green')
    
    # 热辐射箭头
    for i in range(4):
        x_pos = 1.5 + i * 0.6
        arrow = FancyArrowPatch((x_pos, 6), (x_pos, 3.5),
                               arrowstyle='->', mutation_scale=20,
                               color='orange', linewidth=2, alpha=0.6)
        ax3.add_patch(arrow)
    
    # 输出电力
    ax3.annotate('', xy=(7, 3), xytext=(4, 3),
                 arrowprops=dict(arrowstyle='->', color='yellow', lw=3))
    ax3.text(5.5, 3.3, 'Electric Power', fontsize=10, 
            fontweight='bold', color='darkgreen')
    
    # 系统参数
    ax3.text(0.5, 0.5, f'System Parameters:\n' +
                        f'• Emitter: SiC at {T_emitter} K\n' +
                        f'• PV Cell: Optimized bandgap\n' +
                        f'• Gap: 10-100 nm\n' +
                        f'• Near-field enhancement: 100-1000×',
            fontsize=9, verticalalignment='bottom',
            bbox=dict(boxstyle='round', facecolor='lightyellow', alpha=0.8))
    
    # 图4：性能对比（近场vs远场）
    ax4 = axes[1, 1]
    
    categories = ['Efficiency\n(%)', 'Power Density\n(MW/m²)', 'Gap\n(nm)']
    near_field_vals = [35, 50, 50]
    far_field_vals = [15, 0.5, 10000]
    
    x = np.arange(len(categories))
    width = 0.35
    
    bars1 = ax4.bar(x - width/2, near_field_vals, width, label='Near-Field TPV',
                    color='green', alpha=0.7, edgecolor='black')
    bars2 = ax4.bar(x + width/2, far_field_vals, width, label='Far-Field TPV',
                    color='blue', alpha=0.7, edgecolor='black')
    
    ax4.set_ylabel('Value', fontsize=11)
    ax4.set_title('Near-Field vs Far-Field TPV Performance', 
                  fontsize=12, fontweight='bold')
    ax4.set_xticks(x)
    ax4.set_xticklabels(categories)
    ax4.legend(fontsize=10)
    ax4.grid(True, alpha=0.3, axis='y')
    
    # 添加数值标签
    for bar in bars1:
        height = bar.get_height()
        ax4.text(bar.get_x() + bar.get_width()/2., height,
                 f'{height:.0f}', ha='center', va='bottom', fontsize=9)
    
    for bar in bars2:
        height = bar.get_height()
        ax4.text(bar.get_x() + bar.get_width()/2., height,
                 f'{height:.1f}', ha='center', va='bottom', fontsize=9)
    
    plt.tight_layout()
    plt.savefig('case5_near_field_tpv.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print("\n结果已保存至: case5_near_field_tpv.png")
    print("\n案例5完成！")


def case6_super_planck_animation():
    """案例6：超普朗克辐射动态演示"""
    print("\n" + "=" * 70)
    print("案例6：超普朗克辐射动态演示")
    print("=" * 70)
    
    # 温度设置
    T1 = 700  # K
    T2 = 300  # K
    material = 'SiC'
    
    # 间隙距离范围（用于动画）
    gaps = np.logspace(-8, -5, 50)  # 10 nm 到 1 μm
    
    print("\n正在生成动画...")
    
    # 创建图形
    fig = plt.figure(figsize=(14, 10))
    
    # 创建子图布局
    gs = fig.add_gridspec(3, 2, height_ratios=[1, 1, 1])
    
    ax1 = fig.add_subplot(gs[0, :])  # 热流密度vs间隙距离
    ax2 = fig.add_subplot(gs[1, 0])  # 超普朗克因子
    ax3 = fig.add_subplot(gs[1, 1])  # 辐射谱
    ax4 = fig.add_subplot(gs[2, :])  # 物理示意图
    
    # 计算所有数据
    q_values = []
    eta_values = []
    
    for d in gaps:
        q, eta = near_field_heat_flux(d, T1, T2, material)
        q_values.append(q)
        eta_values.append(eta)
    
    q_bb = sigma * (T1**4 - T2**4)
    
    # 初始化线条
    line1, = ax1.plot([], [], 'b-', linewidth=2, label='Near-Field Heat Flux')
    line1_bb, = ax1.plot([], [], 'r--', linewidth=2, label='Blackbody Limit')
    point1, = ax1.plot([], [], 'ro', markersize=10)
    
    line2, = ax2.plot([], [], 'g-', linewidth=2)
    point2, = ax2.plot([], [], 'ro', markersize=10)
    
    # 辐射谱
    wavelengths = np.linspace(1, 20, 200) * 1e-6
    line3, = ax3.plot([], [], 'purple', linewidth=2)
    ax3.axvspan(10, 13, alpha=0.2, color='blue', label='SiC SPhP')
    
    # 设置坐标轴
    ax1.set_xlim(gaps[0] * 1e9, gaps[-1] * 1e9)
    ax1.set_ylim(q_bb * 0.5, max(q_values) * 1.2)
    ax1.set_xscale('log')
    ax1.set_yscale('log')
    ax1.set_xlabel('Gap Distance (nm)', fontsize=11)
    ax1.set_ylabel('Heat Flux (W/m²)', fontsize=11)
    ax1.set_title('Near-Field Heat Flux vs Gap Distance', fontsize=12, fontweight='bold')
    ax1.legend(fontsize=10)
    ax1.grid(True, alpha=0.3, which='both')
    
    ax2.set_xlim(gaps[0] * 1e9, gaps[-1] * 1e9)
    ax2.set_ylim(0.5, max(eta_values) * 1.2)
    ax2.set_xscale('log')
    ax2.set_yscale('log')
    ax2.set_xlabel('Gap Distance (nm)', fontsize=11)
    ax2.set_ylabel('Super-Planck Factor', fontsize=11)
    ax2.set_title('Super-Planck Factor vs Gap Distance', fontsize=12, fontweight='bold')
    ax2.axhline(y=1, color='red', linestyle='--', linewidth=2, label='Blackbody')
    ax2.legend(fontsize=10)
    ax2.grid(True, alpha=0.3)
    
    ax3.set_xlim(1, 20)
    ax3.set_ylim(0, 5e9)
    ax3.set_xlabel('Wavelength (μm)', fontsize=11)
    ax3.set_ylabel('Spectral Radiance (GW/(m²·μm))', fontsize=11)
    ax3.set_title('Radiation Spectrum', fontsize=12, fontweight='bold')
    ax3.legend(fontsize=10)
    ax3.grid(True, alpha=0.3)
    
    # 物理示意图
    ax4.set_xlim(0, 10)
    ax4.set_ylim(0, 3)
    ax4.axis('off')
    ax4.set_title('Near-Field Radiation Physics', fontsize=12, fontweight='bold')
    
    # 添加文本
    text_gap = ax4.text(5, 2.5, '', ha='center', fontsize=12, fontweight='bold',
                       bbox=dict(boxstyle='round', facecolor='yellow', alpha=0.5))
    text_eta = ax4.text(5, 1.5, '', ha='center', fontsize=12, fontweight='bold',
                       bbox=dict(boxstyle='round', facecolor='lightgreen', alpha=0.5))
    
    # 绘制静态元素
    # 表面1
    rect1 = Rectangle((1, 0.5), 3, 0.3, facecolor='lightblue', 
                      edgecolor='blue', linewidth=2)
    ax4.add_patch(rect1)
    ax4.text(2.5, 0.65, f'Surface 1 ({T1} K)', ha='center', fontsize=10)
    
    # 表面2
    rect2 = Rectangle((1, 0), 3, 0.3, facecolor='lightcoral', 
                      edgecolor='red', linewidth=2)
    ax4.add_patch(rect2)
    ax4.text(2.5, 0.15, f'Surface 2 ({T2} K)', ha='center', fontsize=10)
    
    # 光子隧穿示意
    photon_arrows = []
    for i in range(5):
        arrow = FancyArrowPatch((1.5 + i * 0.5, 0.5), (1.5 + i * 0.5, 0.3),
                               arrowstyle='<->', mutation_scale=15,
                               color='purple', linewidth=2, alpha=0.6)
        ax4.add_patch(arrow)
        photon_arrows.append(arrow)
    
    def init():
        line1.set_data([], [])
        line1_bb.set_data([], [])
        point1.set_data([], [])
        line2.set_data([], [])
        point2.set_data([], [])
        line3.set_data([], [])
        text_gap.set_text('')
        text_eta.set_text('')
        return line1, line1_bb, point1, line2, point2, line3, text_gap, text_eta
    
    def animate(frame):
        # 更新热流密度图
        line1.set_data(gaps[:frame+1] * 1e9, q_values[:frame+1])
        line1_bb.set_data(gaps[:frame+1] * 1e9, [q_bb] * (frame+1))
        point1.set_data([gaps[frame] * 1e9], [q_values[frame]])
        
        # 更新超普朗克因子图
        line2.set_data(gaps[:frame+1] * 1e9, eta_values[:frame+1])
        point2.set_data([gaps[frame] * 1e9], [eta_values[frame]])
        
        # 更新辐射谱（根据间隙距离调整强度）
        spectral_power = np.array([planck_law(lam, T1) for lam in wavelengths])
        enhancement = eta_values[frame] / max(eta_values) * 0.5 + 0.5
        line3.set_data(wavelengths * 1e6, spectral_power * enhancement / 1e9)
        
        # 更新文本
        text_gap.set_text(f'Gap Distance: {gaps[frame]*1e9:.1f} nm')
        text_eta.set_text(f'Super-Planck Factor: {eta_values[frame]:.0f}')
        
        # 更新光子隧穿箭头透明度
        alpha = min(1.0, 10 / (gaps[frame] * 1e9))
        for arrow in photon_arrows:
            arrow.set_alpha(alpha)
        
        return line1, line1_bb, point1, line2, point2, line3, text_gap, text_eta
    
    # 创建动画
    anim = animation.FuncAnimation(fig, animate, init_func=init,
                                   frames=len(gaps), interval=200, blit=False)
    
    # 保存动画
    anim.save('super_planck_animation.gif', writer='pillow', fps=5, dpi=100)
    plt.close()
    
    print("动画已保存至: super_planck_animation.gif")
    print("\n案例6完成！")


def main():
    """主函数：运行所有仿真案例"""
    print("\n" + "=" * 70)
    print("主题076：辐射传热极限仿真")
    print("=" * 70)
    
    # 运行所有案例
    case1_blackbody_limit()
    case2_gap_distance_effect()
    case3_material_comparison()
    case4_temperature_effect()
    case5_near_field_tpv()
    case6_super_planck_animation()
    
    print("\n" + "=" * 70)
    print("所有仿真案例已完成！")
    print("=" * 70)
    print("\n生成的文件：")
    print("  1. case1_blackbody_limit.png - 黑体辐射极限验证")
    print("  2. case2_gap_distance_effect.png - 间隙距离对近场传热的影响")
    print("  3. case3_material_comparison.png - 材料选择对超普朗克辐射的影响")
    print("  4. case4_temperature_effect.png - 温度对近场传热的影响")
    print("  5. case5_near_field_tpv.png - 近场热光伏系统性能分析")
    print("  6. super_planck_animation.gif - 超普朗克辐射动态演示")
    print("\n" + "=" * 70)


if __name__ == "__main__":
    main()