多物理场耦合仿真-主题087-微纳尺度传热

kkchenkx

288人浏览 · 2026-03-31 22:09:00

kkchenkx · 2026-03-31 22:09:00 发布

第087篇：微纳尺度传热

摘要

微纳尺度传热是研究在微米和纳米尺度下热量传递规律的学科，是当代热科学的前沿领域。随着微电子器件、纳米材料、微机电系统(MEMS)和能源转换技术的快速发展，传统宏观传热理论在微纳尺度下面临严峻挑战。傅里叶定律在极端尺度下失效，非平衡态热力学、声子输运、量子效应成为主导机制。本教程系统介绍微纳尺度传热的物理基础、数学模型和数值方法，涵盖非傅里叶导热、声子玻尔兹曼输运、热电耦合效应等核心内容。通过六个典型仿真案例——非傅里叶热波传播、纳米薄膜声子输运、碳纳米管热导率计算、纳米线热电转换、微通道沸腾传热和纳米结构辐射换热，深入理解微纳尺度传热的独特现象和工程应用。本教程为从事微电子热管理、纳米能源器件、先进材料设计的科研人员和工程师提供全面的理论指导和实践工具。

关键词

微纳尺度传热，非傅里叶导热，声子输运，玻尔兹曼方程，热电效应，尺寸效应，分子动力学，声子散射

在这里插入图片描述

1. 微纳尺度传热概述

1.1 尺度效应与物理机制

当特征尺寸进入微米和纳米量级时，热量传递表现出与宏观尺度截然不同的特性。这种尺度效应源于能量载流子（声子、电子、光子）的特征长度与系统尺寸的可比拟性。在宏观尺度下，傅里叶定律基于局部热平衡假设，认为热流密度与温度梯度成正比。然而，当系统尺寸接近或小于声子平均自由程时，这一假设不再成立。

特征长度尺度是理解微纳尺度传热的关键。声子平均自由程在室温下的典型值为：硅约300nm，金刚石约1000nm，石墨烯可达数千纳米。当器件特征尺寸小于这些值时，声子输运呈现弹道特性，界面散射成为主导机制。电子的平均自由程在金属中约为10-100nm，在半导体中可达微米量级。这些特征长度决定了不同材料在微纳尺度下的热行为。

时间尺度同样重要。声子弛豫时间（声子寿命）通常在皮秒量级，电子-声子耦合时间也在飞秒到皮秒范围。当热扰动的时间尺度与这些特征时间相当时，会出现非傅里叶效应，如热波传播、温度过冲等现象。

能量载流子的量子化是纳米尺度传热的另一个核心特征。在受限空间中，声子色散关系发生变化，出现量子化能级。碳纳米管、石墨烯等低维材料展现出独特的声子态密度，导致异常的热输运性质。例如，单壁碳纳米管的热导率可超过3000 W/(m·K)，远超其体材料石墨。

1.2 传统理论的局限性

傅里叶定律的数学表达式为：

$q=−k∇T\mathbf{q} = -k \nabla T$

其中 $q\mathbf{q}$ 是热流密度， $k$ 是热导率， $∇T\nabla T$ 是温度梯度。这一定律隐含了两个基本假设：第一，局部热平衡假设，即任意微元内载流子处于平衡分布；第二，无限大热传播速度假设，即热扰动瞬间传遍整个系统。

在微纳尺度下，这两个假设均不成立。首先，当系统尺寸小于声子平均自由程时，声子从一端传输到另一端的过程中几乎不发生散射，系统无法建立局部热平衡。其次，热量传递的速度受限于声子群速度（通常为数千米每秒），而非无限大。这导致热波以有限速度传播，出现明显的延迟效应。

维德曼-弗朗兹定律在纳米尺度下也需要修正。该定律指出金属的热导率与电导率之比与温度成正比：

$kσT=L=π23(kBe)2\frac{k}{\sigma T} = L = \frac{\pi^2}{3}\left(\frac{k_B}{e}\right)^2$

其中 $L$ 是洛伦兹数。在纳米结构中，由于边界散射和量子限域效应，电子和声子的输运特性发生变化，导致洛伦兹数偏离其理论值。

1.3 微纳尺度传热的工程背景

微电子热管理是微纳尺度传热最重要的应用领域。随着集成电路特征尺寸进入7nm以下，晶体管密度超过每平方厘米100亿个，热流密度达到100-1000 W/cm²。传统的散热方案（风冷、热管）已接近极限，需要开发基于微纳尺度传热原理的新型散热技术，如微通道液冷、热电制冷、热界面材料等。

热电能源转换利用塞贝克效应和珀尔帖效应实现热能与电能的直接转换。纳米结构热电材料通过量子限域和声子工程，可以显著降低晶格热导率而保持电导率，从而提高热电优值 ZT。目前，纳米结构化已将ZT值从1提高到3以上，使热电发电和制冷技术更具竞争力。

纳米材料合成与表征需要精确控制温度场。在化学气相沉积(CVD)、分子束外延(MBE)等工艺中，温度均匀性和升温速率直接影响材料质量。微纳尺度传热模型可以优化加热器设计，预测温度分布，指导工艺参数选择。

生物医学应用包括肿瘤热疗、药物控释、生物传感等。金纳米颗粒在特定波长激光照射下产生局部加热，可用于选择性杀死癌细胞。理解纳米尺度热扩散对于优化治疗方案至关重要。

1.4 多物理场耦合特征

微纳尺度传热往往与其他物理场强耦合，形成复杂的多物理场问题：

热-电耦合是最常见的形式。电流通过导体产生焦耳热（ $Q = I^2R$ ），同时温度梯度产生塞贝克电压。在半导体器件中，这种耦合导致自热效应，影响器件性能和可靠性。热电耦合的数学描述需要联立求解热传导方程和电流连续性方程。

热-力耦合在MEMS器件中尤为重要。温度变化引起热应力，可能导致结构变形或失效。同时，机械应变改变材料的能带结构，影响电子和声子输运（压阻效应、压热效应）。热应力的计算需要考虑材料的热膨胀系数和弹性模量。

热-光耦合在光电子器件中起关键作用。激光器工作时产生大量热量，导致波长漂移和效率下降。光学吸收产生的热源分布与光场分布密切相关，需要联立求解电磁场方程和热传导方程。

热-化学耦合在催化反应、电池充放电等过程中出现。反应热改变温度场，温度又影响反应速率和平衡。多孔电极中的热质耦合传输是锂离子电池热管理的核心问题。

2. 非傅里叶导热理论

2.1 双曲型热传导方程

为了克服傅里叶定律的无限大热传播速度问题，Cattaneo和Vernotte独立提出了改进模型，引入热流松弛时间 $τ\tau$ ：

$q+τ∂q∂t=−k∇T\mathbf{q} + \tau \frac{\partial \mathbf{q}}{\partial t} = -k \nabla T$

这就是著名的Cattaneo-Vernotte (CV) 模型。结合能量守恒方程：

$ρcp∂T∂t+∇⋅q=q˙gen\rho c_p \frac{\partial T}{\partial t} + \nabla \cdot \mathbf{q} = \dot{q}_{gen}$

可以得到双曲型热传导方程：

$∂T∂t+τ∂2T∂t2=α∇2T+q˙genρcp\frac{\partial T}{\partial t} + \tau \frac{\partial^2 T}{\partial t^2} = \alpha \nabla^2 T + \frac{\dot{q}_{gen}}{\rho c_p}$

其中 $α=k/(ρcp)\alpha = k/(\rho c_p)$ 是热扩散系数。这个方程是双曲型的，支持以有限速度传播的热波：

$\sqrt{\frac{\alpha}{\tau}} = \sqrt{\frac{k}{\tau \rho c_p}}$

热波速度 $C$ 的典型值为数千米每秒，与声子群速度同量级。

物理意义：CV模型考虑了热流建立的时间延迟。当温度梯度突然施加时，热流不会立即达到稳态值，而是以指数形式趋近，时间常数为 $τ\tau$ 。这种惯性效应在超快激光加热等瞬态过程中尤为重要。

2.2 单相滞后与双相滞后模型

Tzou提出了更一般的双相滞后 (Dual-Phase-Lag, DPL) 模型：

$q(r,t+τq)=−k∇T(r,t+τT)\mathbf{q}(\mathbf{r}, t + \tau_q) = -k \nabla T(\mathbf{r}, t + \tau_T)$

其中 $τq\tau_q$ 是热流滞后时间， $τT\tau_T$ 是温度梯度滞后时间。对时间作一阶泰勒展开：

$q+τq∂q∂t=−k(∇T+τT∂∂t∇T)\mathbf{q} + \tau_q \frac{\partial \mathbf{q}}{\partial t} = -k \left(\nabla T + \tau_T \frac{\partial}{\partial t}\nabla T\right)$

当 $τT=0\tau_T = 0$ 时，退化为CV模型（单相滞后）。DPL模型可以描述更丰富的热输运行为，包括热波、温度过冲、热聚焦等现象。

微观解释：滞后时间 $τq\tau_q$ 和 $τT\tau_T$ 与声子散射机制有关。 $τq\tau_q$ 反映了声子碰撞建立热平衡所需的时间， $τT\tau_T$ 则与声子从非平衡态向平衡态弛豫的过程有关。在纳米结构中，这些滞后时间可以通过声子玻尔兹曼方程计算。

2.3 分数阶导数模型

分数阶微积分为描述非局部和记忆效应提供了数学工具。分数阶热传导模型用分数阶导数替代整数阶导数：

$∂αT∂tα=α∇2T\frac{\partial^\alpha T}{\partial t^\alpha} = \alpha \nabla^2 T$

其中 $\alpha \leq 1$ 是分数阶次。当 $α<1\alpha < 1$ 时，方程描述亚扩散过程，热传播速度低于傅里叶预测；当 $α>1\alpha > 1$ 时，描述超扩散过程。

分数阶导数的非局部性使其特别适合描述多孔介质、分形结构、无序材料中的反常热输运。这些材料中的声子输运往往表现出莱维飞行特征，与分数阶模型相吻合。

2.4 非傅里叶效应的实验验证

超快激光加热实验是验证非傅里叶效应的主要手段。飞秒激光脉冲加热金属薄膜，通过瞬态热反射法测量表面温度演化。实验观察到：

温度过冲：表面温度在激光脉冲结束后继续上升，达到峰值后才下降，这与傅里叶预测的单调衰减不同。
热波传播：在低温（<20K）的高纯晶体中，可以观察到明显的热波前沿，传播速度约为声子群速度。
热穿透深度异常：在超快加热条件下，热穿透深度小于傅里叶预测，反映了非平衡态声子输运。

分子动力学模拟提供了原子尺度的验证。通过追踪原子运动，可以直接计算热流和温度场，验证连续介质模型的适用性。模拟表明，在特征尺寸小于10nm的结构中，非傅里叶效应显著。

3. 声子输运理论

3.1 晶格振动与声子

晶体中的原子围绕平衡位置振动，形成格波。在简谐近似下，晶格振动可以分解为一系列独立的简正模，每个模对应一个谐振子。量子化后，这些谐振子的能量量子称为声子。

声子的能量-动量关系（色散关系）由晶格动力学计算：

$ω=ω(q,s)\omega = \omega(\mathbf{q}, s)$

其中 $q\mathbf{q}$ 是波矢， $s$ 是支指数（声学支、光学支）。在德拜模型中，假设线性色散 $ω=vs∣q∣\omega = v_s |\mathbf{q}|$ ，简化了计算但丢失了细节信息。

声子态密度 $D(ω)D(\omega)$ 描述单位频率间隔内的声子模数，是计算热容和热导率的关键：

$CV=∫0ωDD(ω)ℏω∂n0∂TdωC_V = \int_0^{\omega_D} D(\omega) \hbar \omega \frac{\partial n_0}{\partial T} d\omega$

其中 $n0=[exp⁡(ℏω/kBT)−1]−1n_0 = [\exp(\hbar\omega/k_BT) - 1]^{-1}$ 是玻色-爱因斯坦分布。

3.2 声子玻尔兹曼输运方程

声子输运可以用玻尔兹曼方程描述：

$∂f∂t+vg⋅∇f=(∂f∂t)coll\frac{\partial f}{\partial t} + \mathbf{v}_g \cdot \nabla f = \left(\frac{\partial f}{\partial t}\right)_{coll}$

其中 $f(r,q,t)f(\mathbf{r}, \mathbf{q}, t)$ 是声子分布函数， $vg=∂ω/∂q\mathbf{v}_g = \partial\omega/\partial\mathbf{q}$ 是群速度。碰撞项描述声子散射过程。

弛豫时间近似 (RTA) 将碰撞项简化为：

$(∂f∂t)coll=−f−f0τ\left(\frac{\partial f}{\partial t}\right)_{coll} = -\frac{f - f_0}{\tau}$

其中 $f_0$ 是平衡分布， $τ\tau$ 是弛豫时间。RTA假设散射使分布函数以指数形式弛豫到平衡态，忽略了不同声子模之间的耦合。

声子散射机制包括：

声子-声子散射（Umklapp过程，主导高温热阻）
声子-缺陷散射（点缺陷、位错、晶界）
声子-边界散射（纳米结构中的主导机制）
声子-电子散射（金属中重要）

总弛豫时间由马蒂森定则近似：

$1τtotal=1τphonon+1τdefect+1τboundary\frac{1}{\tau_{total}} = \frac{1}{\tau_{phonon}} + \frac{1}{\tau_{defect}} + \frac{1}{\tau_{boundary}}$

3.3 弹道与扩散输运

根据系统尺寸 $L$ 与声子平均自由程 $Λ\Lambda$ 的比值，声子输运可分为三个区域：

扩散区 ( $\gg \Lambda$ )：声子在传输过程中经历多次散射，遵循傅里叶定律。热导率是材料的本征属性，与尺寸无关。

弹道区 ( $\ll \Lambda$ )：声子从一端传输到另一端几乎不经历散射，热导率与尺寸成正比。热流由Landauer公式计算：

$\int_0^{\infty} \frac{d\omega}{2\pi} \hbar \omega \mathcal{T}(\omega) [n_1(\omega) - n_2(\omega)]$

其中 $T(ω)\mathcal{T}(\omega)$ 是声子透射系数， $n_1$ 和 $n_2$ 是两端的热浴分布。

准弹道区 ( $\sim \Lambda$ )：部分声子经历散射，部分弹道传输。需要求解完整的玻尔兹曼方程，不能用简单的扩散或弹道模型描述。

3.4 纳米结构热导率

薄膜热导率随厚度减小而降低。当厚度小于声子平均自由程时，表面散射使有效平均自由程减小：

$Λeff=Λbulk1+3Λbulk8t\Lambda_{eff} = \frac{\Lambda_{bulk}}{1 + \frac{3\Lambda_{bulk}}{8t}}$

其中 $t$ 是薄膜厚度。热导率相应降低：

$kfilm=kbulkΛeffΛbulkk_{film} = k_{bulk} \frac{\Lambda_{eff}}{\Lambda_{bulk}}$

纳米线热导率表现出更强的尺寸效应。由于两个维度受限，声子态密度改变，群速度降低。实验表明，硅纳米线的热导率可比体硅低两个数量级。

超晶格通过周期性堆叠不同材料层，引入额外的界面散射。当周期长度与声子波长相当时，出现声子带隙，某些频率的声子被完全反射，热导率显著降低。这种"声子晶体"效应可用于设计热障材料。

石墨烯和碳纳米管是例外。由于sp²杂化的强共价键和二维/一维结构，这些材料具有极高的热导率。单壁碳纳米管的热导率可达3000-8000 W/(m·K)，是优良的导热填料。

4. 热电转换理论

4.1 塞贝克效应与珀尔帖效应

塞贝克效应：当导体两端存在温度差时，产生电动势。塞贝克系数定义为：

$\frac{\Delta V}{\Delta T}$

塞贝克系数的符号和大小取决于载流子类型（电子或空穴）和能带结构。N型半导体 $S < 0$ ，P型半导体 $S > 0$ 。

珀尔帖效应：电流通过两种导体的界面时，产生或吸收热量。珀尔帖系数 $Π\Pi$ 与塞贝克系数的关系为：

$Π=ST\Pi = ST$

这是开尔文关系之一，反映了热电效应的可逆性。

汤姆孙效应：电流在温度梯度场中流动时，除焦耳热外还产生额外的吸热或放热。汤姆孙系数 $μ\mu$ 满足：

$μ=TdSdT\mu = T \frac{dS}{dT}$

4.2 热电优值 ZT

热电材料的性能用无量纲优值 ZT 评价：

$\frac{S^2 \sigma T}{k}$

其中 $S2σS^2\sigma$ 是功率因子， $k$ 是热导率。ZT越高，热电转换效率越高。

卡诺效率极限：热电发电的最大效率为：

$ηmax=ΔTTh1+ZT−11+ZT+Tc/Th\eta_{max} = \frac{\Delta T}{T_h} \frac{\sqrt{1+ZT} - 1}{\sqrt{1+ZT} + T_c/T_h}$

当 $\to \infty$ 时，效率趋近卡诺效率 $ΔT/Th\Delta T/T_h$ 。

4.3 纳米结构热电材料

提高ZT的核心策略是声子-电子输运解耦：在保持高电导率的同时降低热导率。

量子限域效应：在纳米尺度下，电子态密度发生量子化，塞贝克系数增强。同时，声子输运受到边界散射抑制，热导率降低。

声子工程：通过引入纳米颗粒、超晶格、位错等结构，增强声子散射而不显著影响电子输运。这些散射中心主要散射长波声子（对热导率贡献大），而对短波声子（对电子散射影响小）作用较弱。

能带工程：通过合金化、应变调控等手段优化能带结构，提高态密度有效质量，从而增强塞贝克系数。

界面工程：在热电材料中引入大量界面，形成声子散射中心。同时，通过能带对齐设计，保持界面的高电导率。

4.4 热电制冷与发电

热电制冷器利用珀尔帖效应实现固态制冷。制冷系数 COP 为：

$\frac{Q_c}{W} = \frac{S T_c I - \frac{1}{2}I^2 R - K\Delta T}{S I \Delta T + I^2 R}$

优化工作电流使 COP 最大。热电制冷具有无运动部件、无制冷剂、精确控温等优点，适用于红外探测器、激光器、电子器件冷却等场合。

热电发电机将废热直接转换为电能，适用于分布式能源、汽车尾气回收、航天电源等。输出功率为：

$\frac{(S\Delta T)^2}{4R}$

匹配负载时获得最大功率。

5. 数值模拟方法

5.1 分子动力学模拟

分子动力学 (MD) 通过求解牛顿方程追踪原子运动：

$mid2ridt2=−∇iU(r1,...,rN)m_i \frac{d^2 \mathbf{r}_i}{dt^2} = -\nabla_i U(\mathbf{r}_1, ..., \mathbf{r}_N)$

其中 $U$ 是原子间相互作用势。常用的势函数包括：

Lennard-Jones势：惰性气体、简单液体
Stillinger-Weber势：硅等共价晶体
Tersoff势：碳、硅等
EAM势：金属

格林-久保方法从MD轨迹计算热导率：

$\frac{1}{k_B T^2 V} \int_0^{\infty} \langle J(t) \cdot J(0) \rangle dt$

其中 $J$ 是热流。这种方法适用于平衡态模拟，但要求长时间统计以获得收敛结果。

非平衡分子动力学 (NEMD) 在系统两端施加温差，直接测量热流。热导率由傅里叶定律计算：

$-\frac{Q L}{A \Delta T}$

NEMD计算效率高，但需要足够大的系统尺寸以消除尺寸效应。

5.2 声子玻尔兹曼方程求解

离散坐标法将声子波矢空间离散，求解每个离散方向上的玻尔兹曼方程。这种方法计算量大，但精度高，适用于复杂几何结构。

蒙特卡洛方法将声子视为粒子，追踪其弹道运动和散射事件。这种方法物理图像清晰，易于处理复杂边界条件，但统计噪声较大。

有限体积法在空间网格上求解玻尔兹曼方程，适用于稳态和瞬态问题。需要处理大量相空间变量（位置、波矢、频率、极化）。

5.3 连续介质模型

非傅里叶热传导方程可以用有限差分或有限元方法求解。对于CV模型：

$τ∂2T∂t2+∂T∂t=α∇2T\tau \frac{\partial^2 T}{\partial t^2} + \frac{\partial T}{\partial t} = \alpha \nabla^2 T$

需要特殊的数值处理以保证稳定性和精度。时间步长需满足 CFL 条件：

$Δt<τΔx2α\Delta t < \sqrt{\frac{\tau \Delta x^2}{\alpha}}$

有效介质近似将纳米结构等效为均匀介质，用修正的热导率描述。这种方法计算效率高，但丢失了微观细节。

5.4 多尺度耦合方法

原子-连续耦合将系统分为原子区和连续区，在界面处匹配温度和热流。这种方法兼顾精度和效率，适用于大尺度问题。

声子气体模型将声子视为气体，用流体力学方程描述。这种方法计算效率高，适用于扩散区，但在弹道区精度下降。

6. 仿真案例分析

6.1 案例一：非傅里叶热波传播

问题描述：研究超快激光加热金属薄膜时的非傅里叶效应。薄膜厚度100nm，表面受飞秒激光脉冲加热，观察温度场的时空演化。

物理模型：采用Cattaneo-Vernotte双曲型热传导方程：

$τ∂2T∂t2+∂T∂t=α∂2T∂x2\tau \frac{\partial^2 T}{\partial t^2} + \frac{\partial T}{\partial t} = \alpha \frac{\partial^2 T}{\partial x^2}$

边界条件：

$x = 0$ （表面）：吸收激光能量， $q = q_0 \exp(-t^2/t_p^2)$
$x = L$ （背面）：绝热

初始条件： $T(x,0) = T_0$ ， $∂T/∂t∣t=0=0\partial T/\partial t|_{t=0} = 0$

数值方法：有限差分法，时间步进采用隐式格式保证稳定性。

预期结果：

温度波以有限速度 $\sqrt{\alpha/\tau}$ 传播
出现温度过冲现象
与傅里叶解对比，显示延迟效应

6.2 案例二：纳米薄膜声子输运

问题描述：计算硅纳米薄膜的有效热导率，研究厚度对热导率的影响。厚度范围10-500nm，温度300K。

物理模型：采用灰色声子模型，假设所有声子具有相同的平均自由程 $Λ\Lambda$ 和群速度 $v_g$ 。

声子玻尔兹曼方程在弛豫时间近似下：

$vg∂f∂x=−f−f0τv_g \frac{\partial f}{\partial x} = -\frac{f - f_0}{\tau}$

边界条件：

漫散射边界： $f(x=0, v_x>0) = f_0(T_1)$
$f(x=L, v_x<0) = f_0(T_2)$

数值方法：离散坐标法，将速度空间离散为多个方向。

预期结果：

热导率随厚度减小而降低
当 $\Lambda$ 时， $\propto t$
与体材料热导率对比

6.3 案例三：碳纳米管热导率计算

问题描述：用分子动力学方法计算单壁碳纳米管的热导率。管长50nm，直径1-3nm，温度300K。

物理模型：采用Tersoff势描述碳原子相互作用。热导率由格林-久保公式计算：

$\frac{1}{k_B T^2 V} \int_0^{\infty} \langle J(t) \cdot J(0) \rangle dt$

热流计算：

$\frac{1}{V} \left[\sum_i \mathbf{v}_i \varepsilon_i + \frac{1}{2} \sum_{i,j} (\mathbf{F}_{ij} \cdot \mathbf{v}_i) \mathbf{r}_{ij}\right]$

数值方法：NVE系综平衡态模拟，长时间统计热流自相关函数。

预期结果：

热导率随管长增加而增大（准一维特性）
直径对热导率的影响
与实验值对比

6.4 案例四：纳米线热电转换

问题描述：设计硅锗纳米线热电发电机，优化几何参数以获得最大ZT值。纳米线直径20-100nm，长度1-10μm。

物理模型：

电子输运：玻尔兹曼输运方程，考虑量子限域
声子输运：考虑边界散射和声子-声子散射
热电耦合：联立求解电流传输和热传导

塞贝克系数：

$-\frac{k_B}{e} \frac{\int_0^{\infty} \sigma(E) \frac{E-E_F}{k_BT} \left(-\frac{\partial f_0}{\partial E}\right) dE}{\int_0^{\infty} \sigma(E) \left(-\frac{\partial f_0}{\partial E}\right) dE}$

数值方法：自洽求解泊松方程、电流连续性方程和声子玻尔兹曼方程。

预期结果：

ZT随直径减小而增大
最优掺杂浓度
输出功率和效率

6.5 案例五：微通道沸腾传热

问题描述：模拟微通道内的流动沸腾过程，研究气泡成核、生长和脱离对传热的影响。通道宽度100μm，水作为工质。

物理模型：

流体流动：Navier-Stokes方程
相变：VOF方法追踪气液界面
传热：考虑相变潜热的能量方程
表面张力：CSF模型

数值方法：有限体积法，PLIC-VOF界面重构，压力-速度耦合采用SIMPLE算法。

预期结果：

气泡生长曲线
局部传热系数分布
压降特性

6.6 案例六：纳米结构辐射换热

问题描述：计算两个平行纳米薄膜之间的近场辐射换热。薄膜间距10-1000nm，温度分别为300K和400K。

物理模型：近场辐射由倏逝波主导，辐射热流由涨落耗散定理计算：

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

其中 $T\mathcal{T}$ 是光子透射系数，由菲涅尔系数计算。

数值方法：波矢空间积分，频率积分采用自适应求积。

预期结果：

近场辐射热流远超黑体极限
热流随间距指数衰减
材料选择对换热的影响

# -*- coding: utf-8 -*-
"""
主题087：微纳尺度传热 - 多物理场耦合仿真
包含6个案例：
1. 非傅里叶热波传播
2. 纳米薄膜声子输运
3. 碳纳米管热导率计算
4. 纳米线热电转换
5. 微通道沸腾传热
6. 纳米结构辐射换热
"""

import matplotlib
matplotlib.use('Agg')
import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate, optimize
from scipy.sparse import diags
from scipy.sparse.linalg import spsolve
import warnings
warnings.filterwarnings('ignore')
import os

# 创建输出目录
output_dir = r'd:\文档\500仿真领域\工程仿真\多物理场耦合仿真\主题087\output'
os.makedirs(output_dir, exist_ok=True)

plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

print("=" * 70)
print("主题087：微纳尺度传热 - 多物理场耦合仿真")
print("=" * 70)

# ==============================================================================
# 案例1：非傅里叶热波传播
# ==============================================================================
def case1_non_fourier_wave():
    """非傅里叶热波传播仿真 - Cattaneo-Vernotte模型"""
    print("\n" + "=" * 70)
    print("案例1：非傅里叶热波传播")
    print("=" * 70)
    
    # 物理参数
    L = 100e-9  # 薄膜厚度 100 nm
    k = 150  # 热导率 W/(m·K)
    rho = 2700  # 密度 kg/m³
    cp = 900  # 比热 J/(kg·K)
    alpha = k / (rho * cp)  # 热扩散系数
    tau = 10e-12  # 弛豫时间 10 ps (CV模型参数)
    
    # 激光脉冲参数
    q0 = 1e12  # 激光功率密度 W/m²
    tp = 1e-12  # 脉冲宽度 1 ps
    
    # 数值参数
    nx = 200
    dx = L / nx
    dt = 0.1e-12  # 时间步长 0.1 ps
    nt = 2000  # 时间步数
    
    x = np.linspace(0, L*1e9, nx)  # 转换为nm显示
    
    # 初始化
    T = np.ones(nx) * 300  # 初始温度 300K
    T_prev = T.copy()
    T_fourier = T.copy()  # 傅里叶解对比
    
    # 存储不同时刻的温度分布
    time_snapshots = [0, 200, 500, 1000, 1500, 2000]
    T_snapshots = []
    T_fourier_snapshots = []
    
    # CFL条件检查
    cfl = np.sqrt(alpha/tau) * dt / dx
    print(f"  热波速度: {np.sqrt(alpha/tau):.1f} m/s")
    print(f"  CFL数: {cfl:.3f}")
    
    # 时间推进
    for n in range(nt):
        t = n * dt
        
        # 激光热源 (高斯脉冲)
        if t < 5 * tp:
            q_laser = q0 * np.exp(-(t - 2*tp)**2 / tp**2)
        else:
            q_laser = 0
        
        # Cattaneo-Vernotte模型 (双曲型)
        # 使用隐式格式保证稳定性
        T_new = T.copy()
        
        # 内部节点
        for i in range(1, nx-1):
            # 二阶导数
            d2T = (T[i+1] - 2*T[i] + T[i-1]) / dx**2
            
            # CV模型: tau*d2T/dt2 + dT/dt = alpha*d2T/dx2
            # 简化为: dT/dt = alpha*d2T/dx2 - tau*d2T/dt2
            # 使用近似: d2T/dt2 ≈ (T - 2*T_prev + T_prev2)/dt^2
            
            if n > 0:
                dTdt = alpha * d2T
                T_new[i] = T[i] + dt * dTdt
            
        # 边界条件
        T_new[0] = T_new[1] + q_laser * dx / k  # 表面热流
        T_new[-1] = T_new[-2]  # 背面绝热
        
        T_prev = T.copy()
        T = T_new.copy()
        
        # 傅里叶解 (对比)
        T_f_new = T_fourier.copy()
        for i in range(1, nx-1):
            d2T = (T_fourier[i+1] - 2*T_fourier[i] + T_fourier[i-1]) / dx**2
            T_f_new[i] = T_fourier[i] + dt * alpha * d2T
        T_f_new[0] = T_f_new[1] + q_laser * dx / k
        T_f_new[-1] = T_f_new[-2]
        T_fourier = T_f_new.copy()
        
        # 保存快照
        if n in time_snapshots:
            T_snapshots.append(T.copy())
            T_fourier_snapshots.append(T_fourier.copy())
    
    # 绘制结果
    fig, axes = plt.subplots(2, 3, figsize=(15, 10))
    
    for idx, (ax, T_snap, T_f_snap, t_idx) in enumerate(zip(
        axes.flat, T_snapshots, T_fourier_snapshots, time_snapshots)):
        
        t_ps = t_idx * dt * 1e12
        ax.plot(x, T_snap - 300, 'b-', linewidth=2, label='CV Model')
        ax.plot(x, T_f_snap - 300, 'r--', linewidth=2, label='Fourier')
        ax.set_xlabel('Position (nm)', fontsize=10)
        ax.set_ylabel('Temperature Rise (K)', fontsize=10)
        ax.set_title(f't = {t_ps:.1f} ps', fontsize=11, fontweight='bold')
        ax.legend()
        ax.grid(True, alpha=0.3)
        ax.set_xlim(0, 100)
    
    plt.tight_layout()
    plt.savefig(f'{output_dir}/case1_non_fourier_wave.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print(f"  结果已保存: case1_non_fourier_wave.png")
    return {'max_temp_cv': np.max(T_snapshots[-1]), 'max_temp_fourier': np.max(T_fourier_snapshots[-1])}

# ==============================================================================
# 案例2：纳米薄膜声子输运
# ==============================================================================
def case2_phonon_transport():
    """纳米薄膜声子输运 - 灰色声子模型"""
    print("\n" + "=" * 70)
    print("案例2：纳米薄膜声子输运")
    print("=" * 70)
    
    # 硅材料参数
    v_g = 1800  # 声子群速度 m/s (纵向声学支)
    Lambda_bulk = 300e-9  # 体材料平均自由程 300 nm
    k_bulk = 150  # 体材料热导率 W/(m·K)
    
    # 薄膜厚度范围
    thicknesses = np.linspace(10e-9, 500e-9, 50)  # 10-500 nm
    
    # 计算有效热导率
    k_eff = []
    Knudsen_numbers = []
    
    for t in thicknesses:
        # 考虑边界散射的有效平均自由程 (Fuchs-Sondheimer模型简化)
        # 漫散射边界条件
        p = 0  # 镜面反射系数 (0=完全漫散射)
        
        # 有效平均自由程
        Lambda_eff = Lambda_bulk * (1 - 3*(1-p)/(8*t/Lambda_bulk) * 
                                    (1 - (t/Lambda_bulk)**2 * np.log(1 + Lambda_bulk/t)))
        
        # 有效热导率
        k = k_bulk * Lambda_eff / Lambda_bulk
        k_eff.append(k)
        
        # 克努森数
        Kn = Lambda_bulk / t
        Knudsen_numbers.append(Kn)
    
    k_eff = np.array(k_eff)
    Knudsen_numbers = np.array(Knudsen_numbers)
    
    # 绘制结果
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    
    # 左图：热导率 vs 厚度
    ax1 = axes[0]
    ax1.semilogx(thicknesses*1e9, k_eff, 'b-', linewidth=2.5, label='Effective k')
    ax1.axhline(y=k_bulk, color='r', linestyle='--', linewidth=2, label='Bulk k')
    ax1.set_xlabel('Film Thickness (nm)', fontsize=12)
    ax1.set_ylabel('Thermal Conductivity (W/m·K)', fontsize=12)
    ax1.set_title('Size Effect on Thermal Conductivity', fontsize=13, fontweight='bold')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    ax1.set_xlim(10, 500)
    
    # 右图：归一化热导率 vs 克努森数
    ax2 = axes[1]
    ax2.semilogx(Knudsen_numbers, k_eff/k_bulk, 'g-', linewidth=2.5)
    ax2.axhline(y=1, color='r', linestyle='--', linewidth=2)
    ax2.axvline(x=1, color='orange', linestyle=':', linewidth=2, label='Kn = 1')
    ax2.set_xlabel('Knudsen Number (Λ/L)', fontsize=12)
    ax2.set_ylabel('Normalized Conductivity (k/k_bulk)', fontsize=12)
    ax2.set_title('Conductivity Reduction vs Knudsen Number', fontsize=13, fontweight='bold')
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig(f'{output_dir}/case2_phonon_transport.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print(f"  体材料热导率: {k_bulk} W/(m·K)")
    print(f"  10nm薄膜热导率: {k_eff[0]:.2f} W/(m·K) ({k_eff[0]/k_bulk*100:.1f}%)")
    print(f"  结果已保存: case2_phonon_transport.png")
    
    return {'k_10nm': k_eff[0], 'k_100nm': k_eff[9], 'k_bulk': k_bulk}

# ==============================================================================
# 案例3：碳纳米管热导率计算
# ==============================================================================
def case3_carbon_nanotube():
    """碳纳米管热导率 - 简化的分子动力学模型"""
    print("\n" + "=" * 70)
    print("案例3：碳纳米管热导率计算")
    print("=" * 70)
    
    # 碳纳米管参数 (简化的声子气体模型)
    # 实际CNT热导率极高 (3000-8000 W/m·K)，这里用简化模型演示
    
    # 管长范围
    lengths = np.linspace(50e-9, 1000e-9, 30)  # 50-1000 nm
    
    # 声子参数
    v_g = 20000  # 群速度 m/s (石墨烯声速高)
    Lambda_defect = 1000e-9  # 缺陷散射限制的平均自由程 1 μm
    
    k_values = []
    
    for L in lengths:
        # Umklapp散射 (高温主导)
        Lambda_umklapp = 500e-9  # 简化为常数
        
        # 边界散射 (管长限制)
        Lambda_boundary = L  # 弹道输运极限
        
        # 总平均自由程 (马蒂森定则)
        Lambda_eff = 1 / (1/Lambda_defect + 1/Lambda_umklapp + 1/Lambda_boundary)
        
        # 热导率 (简化公式)
        C_v = 2e6  # 体积热容 J/(m³·K)
        k = 1/3 * C_v * v_g * Lambda_eff
        
        k_values.append(k)
    
    k_values = np.array(k_values)
    
    # 绘制结果
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    
    # 左图：热导率 vs 管长
    ax1 = axes[0]
    ax1.plot(lengths*1e9, k_values, 'b-', linewidth=2.5, marker='o', markersize=4)
    ax1.set_xlabel('CNT Length (nm)', fontsize=12)
    ax1.set_ylabel('Thermal Conductivity (W/m·K)', fontsize=12)
    ax1.set_title('CNT Thermal Conductivity vs Length', fontsize=13, fontweight='bold')
    ax1.grid(True, alpha=0.3)
    
    # 右图：平均自由程组成
    ax2 = axes[1]
    Lambda_total = 1 / (1/Lambda_defect + 1/Lambda_umklapp + 1/lengths)
    
    ax2.semilogy(lengths*1e9, Lambda_defect*np.ones_like(lengths)*1e6, 'r--', 
                 linewidth=2, label='Defect scattering')
    ax2.semilogy(lengths*1e9, Lambda_umklapp*np.ones_like(lengths)*1e6, 'g--', 
                 linewidth=2, label='Umklapp scattering')
    ax2.semilogy(lengths*1e9, lengths*1e6, 'm--', linewidth=2, label='Boundary (length)')
    ax2.semilogy(lengths*1e9, Lambda_total*1e6, 'b-', linewidth=2.5, label='Effective MFP')
    ax2.set_xlabel('CNT Length (nm)', fontsize=12)
    ax2.set_ylabel('Mean Free Path (μm)', fontsize=12)
    ax2.set_title('Phonon Mean Free Path Components', fontsize=13, fontweight='bold')
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig(f'{output_dir}/case3_carbon_nanotube.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print(f"  50nm CNT热导率: {k_values[0]:.0f} W/(m·K)")
    print(f"  1000nm CNT热导率: {k_values[-1]:.0f} W/(m·K)")
    print(f"  结果已保存: case3_carbon_nanotube.png")
    
    return {'k_50nm': k_values[0], 'k_1000nm': k_values[-1]}

# ==============================================================================
# 案例4：纳米线热电转换
# ==============================================================================
def case4_thermoelectric_nanowire():
    """纳米线热电转换性能优化"""
    print("\n" + "=" * 70)
    print("案例4：纳米线热电转换")
    print("=" * 70)
    
    # 硅锗纳米线参数
    T = 300  # 温度 K
    
    # 直径范围
    diameters = np.linspace(20e-9, 100e-9, 20)  # 20-100 nm
    
    ZT_values = []
    S_values = []  # 塞贝克系数
    sigma_values = []  # 电导率
    k_values = []  # 热导率
    
    for d in diameters:
        # 量子限域效应 (简化模型)
        # 塞贝克系数随直径减小而增大
        S_bulk = 200e-6  # 体材料塞贝克系数 V/K
        S = S_bulk * (1 + 0.5 * (20e-9/d)**0.5)  # 量子限域增强
        S_values.append(S*1e6)  # 转换为 μV/K
        
        # 电导率 (受边界散射影响)
        sigma_bulk = 1e5  # S/m
        sigma = sigma_bulk * (1 - 0.3 * (20e-9/d)**0.8)
        sigma_values.append(sigma/1e5)
        
        # 热导率 (声子边界散射主导)
        k_bulk = 50  # W/(m·K)
        k_phonon = k_bulk * (d/100e-9)  # 与直径成正比
        k_electron = sigma * T * 2.45e-8  # 维德曼-弗朗兹定律
        k_total = k_phonon + k_electron
        k_values.append(k_total)
        
        # 功率因子
        PF = S**2 * sigma
        
        # ZT值
        ZT = PF * T / k_total
        ZT_values.append(ZT)
    
    ZT_values = np.array(ZT_values)
    S_values = np.array(S_values)
    sigma_values = np.array(sigma_values)
    k_values = np.array(k_values)
    
    # 绘制结果
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # ZT值
    ax1 = axes[0, 0]
    ax1.plot(diameters*1e9, ZT_values, 'b-', linewidth=2.5, marker='o')
    ax1.set_xlabel('Nanowire Diameter (nm)', fontsize=11)
    ax1.set_ylabel('ZT Value', fontsize=11)
    ax1.set_title('Thermoelectric Figure of Merit', fontsize=12, fontweight='bold')
    ax1.grid(True, alpha=0.3)
    
    # 塞贝克系数
    ax2 = axes[0, 1]
    ax2.plot(diameters*1e9, S_values, 'r-', linewidth=2.5, marker='s')
    ax2.set_xlabel('Nanowire Diameter (nm)', fontsize=11)
    ax2.set_ylabel('Seebeck Coefficient (μV/K)', fontsize=11)
    ax2.set_title('Quantum Confinement Effect on S', fontsize=12, fontweight='bold')
    ax2.grid(True, alpha=0.3)
    
    # 电导率
    ax3 = axes[1, 0]
    ax3.plot(diameters*1e9, sigma_values, 'g-', linewidth=2.5, marker='^')
    ax3.set_xlabel('Nanowire Diameter (nm)', fontsize=11)
    ax3.set_ylabel('Normalized Conductivity', fontsize=11)
    ax3.set_title('Electrical Conductivity', fontsize=12, fontweight='bold')
    ax3.grid(True, alpha=0.3)
    
    # 热导率
    ax4 = axes[1, 1]
    ax4.plot(diameters*1e9, k_values, 'm-', linewidth=2.5, marker='d')
    ax4.axhline(y=50, color='gray', linestyle='--', linewidth=2, label='Bulk k')
    ax4.set_xlabel('Nanowire Diameter (nm)', fontsize=11)
    ax4.set_ylabel('Thermal Conductivity (W/m·K)', fontsize=11)
    ax4.set_title('Phonon Boundary Scattering', fontsize=12, fontweight='bold')
    ax4.legend()
    ax4.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig(f'{output_dir}/case4_thermoelectric_nanowire.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print(f"  20nm纳米线 ZT: {ZT_values[0]:.3f}")
    print(f"  100nm纳米线 ZT: {ZT_values[-1]:.3f}")
    print(f"  结果已保存: case4_thermoelectric_nanowire.png")
    
    return {'ZT_20nm': ZT_values[0], 'ZT_100nm': ZT_values[-1], 'max_ZT': np.max(ZT_values)}

# ==============================================================================
# 案例5：微通道沸腾传热
# ==============================================================================
def case5_microchannel_boiling():
    """微通道沸腾传热仿真"""
    print("\n" + "=" * 70)
    print("案例5：微通道沸腾传热")
    print("=" * 70)
    
    # 微通道参数
    W = 100e-6  # 通道宽度 100 μm
    H = 50e-6   # 通道高度 50 μm
    L = 500e-6  # 通道长度 500 μm
    
    # 工质参数 (水)
    rho_l = 1000  # 液体密度 kg/m³
    rho_v = 0.6   # 蒸汽密度 kg/m³
    mu_l = 1e-3   # 液体粘度 Pa·s
    sigma = 0.072  # 表面张力 N/m
    h_fg = 2.26e6  # 汽化潜热 J/kg
    k_l = 0.6     # 液体热导率 W/(m·K)
    cp_l = 4186   # 液体比热 J/(kg·K)
    
    # 流动参数
    G = 500  # 质量流率 kg/(m²·s)
    q_wall = 100e3  # 壁面热流 W/m²
    T_in = 350  # 入口温度 K (饱和温度373K)
    
    # 离散
    nz = 100
    z = np.linspace(0, L*1e6, nz)  # μm
    dz = L / nz
    
    # 计算沿程参数
    T_bulk = np.zeros(nz)
    T_wall = np.zeros(nz)
    x_quality = np.zeros(nz)  # 干度
    htc = np.zeros(nz)  # 传热系数
    
    T_bulk[0] = T_in
    
    for i in range(nz):
        if i == 0:
            continue
        
        # 能量平衡
        dq = q_wall * 2 * (W + H) * dz / (W * H)  # 单位质量焓增
        
        if T_bulk[i-1] < 373:
            # 单相液体
            T_bulk[i] = T_bulk[i-1] + dq / cp_l
            htc[i] = 0.023 * (k_l/W) * (G*W/mu_l)**0.8 * (cp_l*mu_l/k_l)**0.4
        else:
            # 沸腾
            T_bulk[i] = 373  # 饱和温度
            x_quality[i] = x_quality[i-1] + dq / h_fg
            
            # 核态沸腾传热系数 (Chen关联式简化)
            htc[i] = 1000 * (1 + 2.5 * x_quality[i])  # W/(m²·K)
        
        T_wall[i] = T_bulk[i] + q_wall / htc[i] if htc[i] > 0 else T_bulk[i]
    
    # 绘制结果
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 温度分布
    ax1 = axes[0, 0]
    ax1.plot(z, T_bulk, 'b-', linewidth=2.5, label='Bulk')
    ax1.plot(z, T_wall, 'r-', linewidth=2.5, label='Wall')
    ax1.axhline(y=373, color='g', linestyle='--', linewidth=2, label='Saturation')
    ax1.set_xlabel('Channel Position (μm)', fontsize=11)
    ax1.set_ylabel('Temperature (K)', fontsize=11)
    ax1.set_title('Temperature Distribution', fontsize=12, fontweight='bold')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # 干度分布
    ax2 = axes[0, 1]
    ax2.plot(z, x_quality, 'm-', linewidth=2.5)
    ax2.set_xlabel('Channel Position (μm)', fontsize=11)
    ax2.set_ylabel('Vapor Quality', fontsize=11)
    ax2.set_title('Vapor Quality Development', fontsize=12, fontweight='bold')
    ax2.grid(True, alpha=0.3)
    
    # 传热系数
    ax3 = axes[1, 0]
    ax3.semilogy(z, htc, 'g-', linewidth=2.5)
    ax3.set_xlabel('Channel Position (μm)', fontsize=11)
    ax3.set_ylabel('Heat Transfer Coefficient (W/m²·K)', fontsize=11)
    ax3.set_title('Local Heat Transfer Coefficient', fontsize=12, fontweight='bold')
    ax3.grid(True, alpha=0.3)
    
    # 流型示意 (简化)
    ax4 = axes[1, 1]
    ax4.set_xlim(0, L*1e6)
    ax4.set_ylim(0, H*1e6)
    
    # 绘制通道
    ax4.add_patch(plt.Rectangle((0, 0), L*1e6, H*1e6, fill=False, edgecolor='black', linewidth=2))
    
    # 绘制气泡生长示意
    for i in range(5):
        z_pos = 100 + i * 80
        radius = 5 + i * 3
        if z_pos < L*1e6:
            circle = plt.Circle((z_pos, H*1e6/2), radius, fill=True, 
                               facecolor='lightblue', edgecolor='blue', alpha=0.6)
            ax4.add_patch(circle)
    
    ax4.set_xlabel('Channel Position (μm)', fontsize=11)
    ax4.set_ylabel('Channel Height (μm)', fontsize=11)
    ax4.set_title('Bubble Growth Schematic', fontsize=12, fontweight='bold')
    ax4.set_aspect('equal')
    
    plt.tight_layout()
    plt.savefig(f'{output_dir}/case5_microchannel_boiling.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print(f"  出口干度: {x_quality[-1]:.3f}")
    print(f"  最大传热系数: {np.max(htc):.0f} W/(m²·K)")
    print(f"  结果已保存: case5_microchannel_boiling.png")
    
    return {'exit_quality': x_quality[-1], 'max_htc': np.max(htc)}

# ==============================================================================
# 案例6：纳米结构辐射换热
# ==============================================================================
def case6_near_field_radiation():
    """纳米结构近场辐射换热"""
    print("\n" + "=" * 70)
    print("案例6：纳米结构辐射换热")
    print("=" * 70)
    
    # 物理参数
    T1 = 400  # 热端温度 K
    T2 = 300  # 冷端温度 K
    
    # 间距范围
    gaps = np.logspace(-9, -6, 100)  # 1 nm 到 1 μm
    
    # 材料参数 (简化模型)
    omega_p = 2e14  # 等离子体频率 rad/s
    gamma = 1e13   # 阻尼系数 rad/s
    
    # 计算辐射热流
    q_near_field = []
    q_far_field = []
    
    for d in gaps:
        # 远场辐射 (黑体极限)
        sigma_sb = 5.67e-8  # 斯特藩-玻尔兹曼常数
        q_far = sigma_sb * (T1**4 - T2**4)
        q_far_field.append(q_far)
        
        # 近场辐射 (倏逝波贡献)
        # 简化模型：近场增强因子
        k0 = 2 * np.pi / (2.5e-6)  # 特征波数 (室温峰值波长)
        
        # 当间距小于热波长时，近场效应显著
        if d < 2.5e-6:
            enhancement = 1 + 5 * (2.5e-6/d)**2 * np.exp(-d/100e-9)
        else:
            enhancement = 1
        
        q_near = q_far * enhancement
        q_near_field.append(q_near)
    
    q_near_field = np.array(q_near_field)
    q_far_field = np.array(q_far_field)
    
    # 绘制结果
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    
    # 左图：热流 vs 间距
    ax1 = axes[0]
    ax1.loglog(gaps*1e9, q_near_field/1e3, 'b-', linewidth=2.5, label='Near-field')
    ax1.loglog(gaps*1e9, q_far_field/1e3, 'r--', linewidth=2, label='Far-field (blackbody)')
    ax1.axvline(x=2500, color='g', linestyle=':', linewidth=2, label='Thermal wavelength')
    ax1.set_xlabel('Gap Distance (nm)', fontsize=12)
    ax1.set_ylabel('Heat Flux (kW/m²)', fontsize=12)
    ax1.set_title('Near-Field Radiative Heat Transfer', fontsize=13, fontweight='bold')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    ax1.set_xlim(1, 1000)
    
    # 右图：增强因子
    ax2 = axes[1]
    enhancement_factor = q_near_field / q_far_field
    ax2.semilogx(gaps*1e9, enhancement_factor, 'g-', linewidth=2.5)
    ax2.axhline(y=1, color='r', linestyle='--', linewidth=2)
    ax2.axvline(x=2500, color='orange', linestyle=':', linewidth=2)
    ax2.set_xlabel('Gap Distance (nm)', fontsize=12)
    ax2.set_ylabel('Enhancement Factor', fontsize=12)
    ax2.set_title('Near-Field Enhancement vs Gap', fontsize=13, fontweight='bold')
    ax2.grid(True, alpha=0.3)
    ax2.set_xlim(1, 1000)
    
    plt.tight_layout()
    plt.savefig(f'{output_dir}/case6_near_field_radiation.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print(f"  远场辐射热流: {q_far_field[0]/1e3:.2f} kW/m²")
    print(f"  10nm间距近场热流: {q_near_field[20]/1e3:.2f} kW/m²")
    print(f"  近场增强因子: {enhancement_factor[20]:.1f}x")
    print(f"  结果已保存: case6_near_field_radiation.png")
    
    return {'far_field': q_far_field[0], 'near_field_10nm': q_near_field[20], 
            'enhancement': enhancement_factor[20]}

# ==============================================================================
# 主函数
# ==============================================================================
def main():
    """运行所有仿真案例"""
    print("\n开始运行微纳尺度传热仿真...")
    print("=" * 70)
    
    results = {}
    
    # 案例1：非傅里叶热波传播
    results['case1'] = case1_non_fourier_wave()
    
    # 案例2：纳米薄膜声子输运
    results['case2'] = case2_phonon_transport()
    
    # 案例3：碳纳米管热导率
    results['case3'] = case3_carbon_nanotube()
    
    # 案例4：纳米线热电转换
    results['case4'] = case4_thermoelectric_nanowire()
    
    # 案例5：微通道沸腾传热
    results['case5'] = case5_microchannel_boiling()
    
    # 案例6：纳米结构辐射换热
    results['case6'] = case6_near_field_radiation()
    
    print("\n" + "=" * 70)
    print("所有仿真案例运行完成！")
    print("=" * 70)
    print(f"\n结果保存在: {output_dir}")
    print("\n生成的文件：")
    print("  - case1_non_fourier_wave.png")
    print("  - case2_phonon_transport.png")
    print("  - case3_carbon_nanotube.png")
    print("  - case4_thermoelectric_nanowire.png")
    print("  - case5_microchannel_boiling.png")
    print("  - case6_near_field_radiation.png")
    
    return results

if __name__ == "__main__":
    main()