主题095：辐射换热与传质耦合现象

摘要

辐射换热与传质耦合是工程热物理领域的重要研究方向，涉及蒸发冷却、燃烧过程、大气科学、材料加工等众多应用场景。本主题将系统介绍辐射与传质耦合的基本理论、数学模型和数值方法，重点探讨辐射对相变传质的影响、参与性介质中的辐射-扩散耦合、以及多孔介质中的辐射-对流传质问题。通过Python仿真程序，我们将模拟液滴蒸发过程中的辐射效应、湿空气在辐射场中的传热传质、以及工业干燥过程中的耦合现象，为理解这一复杂的多物理场耦合问题提供全面的数值实验平台。

关键词

辐射换热，传质耦合，蒸发冷却，湿空气，参与性介质，多孔介质，相变传热，数值模拟，扩散-辐射耦合

---

1. 引言

1.1 辐射与传质耦合的物理背景

在许多工程应用中，辐射换热与质量传递过程紧密耦合，形成复杂的多物理场相互作用。这种耦合现象普遍存在于自然界和工业过程中：

蒸发冷却过程：太阳辐射加热水面或湿润表面，同时水分蒸发带走潜热，形成辐射-蒸发耦合。这是大气水循环、冷却塔、喷雾冷却等过程的核心机制。

燃烧与气化：燃料颗粒在高温辐射场中加热，同时发生热解、气化等质量传递过程。辐射传热显著影响燃烧速率和污染物生成。

干燥过程：物料在干燥器中吸收辐射热，内部水分向外扩散并蒸发。辐射-传质耦合决定了干燥效率和产品质量。

大气科学：太阳辐射驱动地球表面蒸发，水蒸气向大气扩散，形成云和降水。辐射-对流-扩散耦合是气候模型的核心。

生物系统：植物叶片吸收太阳辐射进行光合作用，同时通过蒸腾作用散失水分。辐射-传质耦合影响生态系统能量平衡。

1.2 耦合机制的分类

根据辐射与传质相互作用的方式，可将耦合现象分为以下几类：

类型I：辐射驱动传质
辐射作为能量输入，驱动相变或化学反应，进而引起质量传递。例如：太阳辐射驱动海水蒸发、红外辐射加热促进溶剂挥发。

类型II：传质影响辐射
质量传递改变介质的辐射特性，进而影响辐射传输。例如：蒸发增加空气湿度，改变大气辐射特性；燃烧产生碳烟，增强火焰辐射。

类型III：双向强耦合
辐射和传质相互强烈影响，形成反馈回路。例如：燃烧过程中，辐射加热促进燃料蒸发，蒸发产物改变火焰辐射特性，进而影响燃烧温度。

1.3 本主题的研究内容

本主题将系统研究辐射换热与传质耦合现象，主要内容包括：

1. 辐射与蒸发传质耦合：建立液滴蒸发模型，分析辐射对蒸发速率的影响，探讨蒸发冷却效应。

2. 湿空气中的辐射-扩散耦合：研究水蒸气在辐射场中的扩散行为，分析湿度对辐射传输的影响。

3. 多孔介质中的辐射-对流传质：建立多孔介质干燥模型，分析辐射加热对内部水分迁移的影响。

4. 参与性介质中的辐射-化学反应耦合：研究辐射对化学反应速率的影响，分析光化学反应机制。

5. 工业应用案例分析：喷雾冷却、喷雾干燥、冷却塔、太阳能蒸馏等典型应用。

6. 数值方法与仿真技术：介绍耦合问题的数值求解方法，包括算子分裂法、全耦合求解等。

---

2. 辐射与蒸发传质耦合基础

2.1 蒸发传质基本理论

蒸发是液体表面分子获得足够能量克服分子间作用力而进入气相的过程。蒸发传质速率由以下因素决定：

传质驱动力：
蒸发速率正比于液面蒸汽压与环境中蒸汽压之差：

$$\dot{m}=h_{m}A({\rho }_{v,s}-{\rho }_{v,\infty })$$

其中$h_m$是传质系数，$A$是表面积，${\rho }_{v,s}$和${\rho }_{v,\infty }$分别是液面和环境的水蒸气密度。

传质系数：
传质系数与流动状态、几何形状有关。对于球体在静止空气中的蒸发：

$$Sh=\frac{h_{m}d}{D}=2+0.6R{e}^{1/2}S{c}^{1/3}$$

其中$Sh$是Sherwood数，$D$是扩散系数，$Re$是Reynolds数，$Sc$是Schmidt数。

蒸汽压与温度的关系：
液面蒸汽压由Clausius-Clapeyron方程给出：

$$p_{v,s}=p_{sat}(T_s)=p_0\exp \left(-\frac{L_v{M}_w}{R}\left(\frac{1}{T_s}-\frac{1}{T_0}\right)\right)$$

其中$L_v$是汽化潜热，$M_w$是水分子摩尔质量。

2.2 辐射对蒸发的影响

辐射加热通过以下机制影响蒸发过程：

表面温度升高：
辐射热流$q_{rad}$加热液面，使其温度$T_s$高于环境温度$T_{\infty }$：

$$q_{rad}+q_{conv}=\dot{m}{L}_v+q_{cond}$$

其中$q_{conv}$是对流热流，$q_{cond}$是向液体内部的导热。

蒸汽压增加：
根据Clausius-Clapeyron关系，表面温度升高导致蒸汽压指数增加，传质驱动力显著增强。

Marangoni效应：
表面温度不均匀引起表面张力梯度，产生Marangoni对流，增强传质。

辐射穿透效应：
对于半透明液体（如水），辐射可穿透一定深度，产生体积加热效应，改变内部温度分布。

2.3 液滴蒸发模型

考虑辐射影响的液滴蒸发模型包括以下控制方程：

质量守恒：
$$\frac{dm}{dt}=-\dot{m}=-h_m\pi d^2({\rho }_{v,s}-{\rho }_{v,\infty })$$

能量守恒：
$$m{c}_p\frac{dT}{dt}=q_{rad}+q_{conv}-\dot{m}{L}_v$$

液滴直径变化：
假设液滴保持球形，质量变化导致直径变化：
$$\frac{dd}{dt}=-\frac{2\dot{m}}{\pi {\rho }_l{d}^2}$$

辐射热流计算：
对于灰体液滴，辐射吸收为：
$$q_{rad}=\alpha \pi d^{2}G-ϵ\pi d^2\sigma T_s^4$$

其中$G$是入射辐射，$\alpha$和$ϵ$分别是吸收率和发射率。

2.4 蒸发冷却效应

蒸发过程伴随潜热吸收，产生冷却效应：

湿球温度：
当辐射热流与蒸发冷却平衡时，液滴达到稳态温度（湿球温度）：

$$T_{wet}=T_{\infty }-\frac{L_v}{c_p}\frac{(p_{v,s}-p_{v,\infty })}{p}\frac{P{r}^{2/3}}{L{e}^{2/3}}$$

其中$Le=\alpha /D$是Lewis数。

辐射对湿球温度的影响：
入射辐射使湿球温度升高：
$$T_{wet,rad}=T_{wet}+\frac{q_{rad}}{h_{eff}}$$

其中$h_{eff}$是有效传热系数。

冷却效率：
蒸发冷却效率定义为：
$${\eta }_{cooling}=\frac{\dot{m}{L}_v}{q_{rad}+q_{conv}}$$

理想情况下，所有输入热量都用于蒸发，${\eta }_{cooling}=1$。

---

3. 湿空气中的辐射-扩散耦合

3.1 湿空气的热物理性质

湿空气是干空气和水蒸气的混合物，其性质随湿度变化：

密度：
$$\rho =\frac{p}{R_{a}T}(1-x_v(1-\frac{R_a}{R_v}))$$

其中$x_v$是水蒸气摩尔分数，$R_a$和$R_v$分别是干空气和水蒸气的气体常数。

比热容：
$$c_p=c_{p,a}(1-\omega )+c_{p,v}\omega$$

其中$\omega$是含湿量（kg水蒸气/kg干空气）。

扩散系数：
水蒸气在空气中的扩散系数随温度变化：
$$D=D_0{\left(\frac{T}{T_0}\right)}^{1.5}\frac{p_0}{p}$$

3.2 湿空气中的辐射传输

水蒸气是红外辐射的强吸收体，对辐射传输有显著影响：

水蒸气的吸收带：

• 2.7 μm带（强吸收）
• 6.3 μm带（强吸收）
• 转动带（>15 μm，连续吸收）

吸收系数模型：
水蒸气吸收系数可用指数宽谱带模型描述：
$${\kappa }_{\lambda }={\kappa }_0\exp \left(-\frac{E}{k_{B}T}\right){\rho }_v$$

辐射传输方程：
在湿空气中，辐射传输方程为：
$$\frac{d{I}_{\lambda }}{ds}=-{\kappa }_{\lambda }{I}_{\lambda }+{\kappa }_{\lambda }{B}_{\lambda }(T)$$

其中${\kappa }_{\lambda }$是湿空气的总吸收系数（干空气+水蒸气）。

3.3 辐射-扩散耦合方程

考虑辐射影响的湿空气传热传质控制方程：

连续性方程：
$$\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho \mathbf{u})=0$$

动量方程：
$$\rho \frac{D\mathbf{u}}{Dt}=-\nabla p+\mu {\nabla }^2\mathbf{u}+\rho \mathbf{g}$$

能量方程：
$$\rho c_p\frac{DT}{Dt}=k{\nabla }^{2}T+{\dot{q}}_{rad}+{\dot{q}}_{latent}$$

其中${\dot{q}}_{rad}$是辐射热源项，${\dot{q}}_{latent}$是相变潜热。

水蒸气扩散方程：
$$\frac{\partial {\rho }_v}{\partial t}+\nabla \cdot ({\rho }_v\mathbf{u})=\nabla \cdot (D\nabla {\rho }_v)+{\dot{m}}_{evap}$$

辐射热源项：
$${\dot{q}}_{rad}=-\nabla \cdot {\mathbf{q}}_{rad}={\int }_0^{\infty }{\kappa }_{\lambda }(4\pi B_{\lambda }-G_{\lambda })d\lambda$$

3.4 边界层内的耦合效应

在壁面边界层内，辐射-扩散耦合效应尤为显著：

温度边界层：
辐射加热改变温度分布，影响自然对流强度：
$${\delta }_T\sim {\left(\frac{kL}{h_{rad}}\right)}^{1/2}$$

其中$h_{rad}=4ϵ\sigma T^3$是辐射传热系数。

浓度边界层：
温度变化影响扩散系数和蒸汽压，进而改变浓度分布。

耦合参数：
定义辐射-扩散耦合参数：
$$N_{RD}=\frac{q_{rad}}{\rho D{L}_v(\partial \omega /\partial y{)}_w}$$

当$N_{RD}\gg 1$时，辐射主导传质过程。

---

4. 多孔介质中的辐射-对流传质

4.1 多孔介质干燥基础

多孔介质干燥是涉及传热、传质和流动的复杂过程：

干燥阶段：

1. 预热阶段：物料温度升高至湿球温度
2. 恒速干燥阶段：表面维持饱和，干燥速率恒定
3. 降速干燥阶段：内部水分迁移控制，干燥速率下降
4. 平衡阶段：物料达到平衡含水率

水分存在形式：

• 自由水：大孔隙中的液态水
• 毛细水：毛细作用力保持的水分
• 结合水：吸附在固体表面的分子层

水分迁移机制：

• 液态扩散：浓度梯度驱动的液态水流动
• 气态扩散：水蒸气在孔隙中的扩散
• 毛细流动：表面张力驱动的流动
• 压力驱动流：内外压差引起的流动

4.2 辐射加热对干燥的影响

辐射加热改变传统对流传热的干燥特性：

体积加热效应：
辐射穿透多孔介质一定深度，产生体积热源：
$${\dot{q}}_{rad}(x)=\alpha G\exp (-\alpha x)$$

其中$\alpha$是多孔介质的有效吸收系数。

温度分布改变：
体积加热使内部温度高于表面温度，形成"逆温度梯度"，促进内部水分向表面迁移。

干燥速率增强：
辐射干燥通常比对流干燥快2-5倍，特别是对于厚层物料。

选择性加热：
不同组分对辐射的吸收不同，可实现选择性加热（如红外干燥中水分选择性吸收）。

4.3 多孔介质辐射-传质模型

考虑辐射的多孔介质干燥数学模型：

能量方程：
$$(\rho c_p{)}_{eff}\frac{\partial T}{\partial t}=k_{eff}{\nabla }^{2}T+{\dot{q}}_{rad}-{\dot{m}}_{evap}{L}_v$$

其中$(\rho c_p{)}_{eff}$和$k_{eff}$是有效热物性参数。

水分传输方程：
$$\frac{\partial X}{\partial t}=\nabla \cdot (D_{eff}\nabla X)-\frac{{\dot{m}}_{evap}}{{\rho }_s}$$

其中$X$是干基含水率，${\rho }_s$是干固体密度。

蒸汽传输方程：
$$\epsilon \frac{\partial {\rho }_v}{\partial t}=\nabla \cdot (D_{v,eff}\nabla {\rho }_v)+{\dot{m}}_{evap}$$

其中$\epsilon$是孔隙率。

辐射传输方程：
对于半透明多孔介质：
$$\frac{dI}{ds}=-\kappa I+\kappa \sigma T^4/\pi$$

4.4 有效辐射特性

多孔介质的辐射特性与组分、孔隙结构有关：

吸收系数：
$$\kappa =(1-\epsilon ){\kappa }_s+\epsilon {\kappa }_g$$

其中${\kappa }_s$和${\kappa }_g$分别是固体骨架和气体的吸收系数。

散射系数：
孔隙结构引起辐射散射：
$${\sigma }_s=\frac{3}{4}\frac{(1-\epsilon )}{d_p}{Q}_s$$

其中$d_p$是特征孔径，$Q_s$是散射效率。

消光系数：
$$\beta =\kappa +{\sigma }_s$$

穿透深度：
辐射穿透深度定义为强度衰减至表面值的1/e时的距离：
$$\delta =1/\beta$$

---

5. 参与性介质中的辐射-化学反应耦合

5.1 光化学反应基础

辐射可以激发化学反应，形成辐射-化学反应耦合：

光化学反应类型：

1. 光解离：分子吸收光子分解为碎片
   $$AB+h\nu \to A+B$$

2. 光异构化：分子吸收光子改变结构
   $$A+h\nu \to A^{\ast }$$

3. 光催化：辐射激活催化剂促进反应
   $$Cat+h\nu \to Ca{t}^{\ast }$$

量子产率：
量子产率定义为：
$$\Phi =\frac{\text{反应分子数}}{\text{吸收光子数}}$$

反应速率：
光化学反应速率与辐射强度成正比：
$$r=\Phi {\int }_{\lambda }{\kappa }_{\lambda }{G}_{\lambda }d\lambda$$

5.2 辐射对热化学反应的影响

除光化学反应外，辐射还通过热效应影响化学反应：

Arrhenius定律：
反应速率常数与温度关系：
$$k=A\exp \left(-\frac{E_a}{RT}\right)$$

辐射加热效应：
辐射提高反应区温度，指数增加反应速率。

热辐射反馈：
放热反应的辐射增强反过来加热反应物，形成正反馈。

5.3 燃烧中的辐射-化学反应耦合

燃烧是最典型的辐射-化学反应耦合过程：

热解过程：
固体燃料在辐射加热下热解：
$$Fuel\stackrel{rad}{\to }Char+Volatiles$$

挥发分燃烧：
挥发分与氧气反应释放热量和辐射：
$$Volatiles+O_2\to C{O}_2+H_{2}O+h\nu$$

碳烟生成与氧化：
碳烟颗粒吸收和发射辐射，影响火焰温度分布：
$$C_{soot}+O_2\stackrel{rad}{\to }C{O}_2$$

辐射熄火/点火：
辐射散热可能导致火焰熄火，或辐射预热促进点火。

5.4 光催化反应工程

光催化是利用辐射驱动化学反应的重要技术：

TiO₂光催化：
紫外光激发TiO₂产生电子-空穴对：
$$Ti{O}_2+h\nu \to e^{-}+h^{+}$$

反应机理：
$$h^{+}+H_{2}O\to O{H}^{\bullet }+H^{+}$$
$$e^{-}+O_2\to O_2^{-\bullet }$$

辐射传输模型：
光催化反应器中辐射分布：
$$\frac{dI}{dx}=-(\kappa +{\sigma }_s)I$$

反应-辐射耦合：
反应速率依赖于局部辐射强度，而反应物浓度影响辐射吸收。

---

6. 工业应用案例分析

6.1 喷雾冷却

喷雾冷却是利用液滴蒸发带走热量的高效冷却技术：

工作原理：

1. 液滴喷射到热表面
2. 液滴吸收表面热量和辐射热
3. 液滴蒸发，产生强烈冷却
4. 蒸汽通过对流扩散离开

辐射增强效应：
在高温环境中，辐射可显著增强喷雾冷却效果：

• 液滴吸收环境辐射，温度升高
• 表面发射辐射被液滴吸收
• 蒸发速率提高2-3倍

关键参数：

• 液滴直径：50-200 μm
• 喷射速度：10-50 m/s
• 热流密度：可达10 MW/m²

6.2 喷雾干燥

喷雾干燥将液体物料雾化成液滴，在热气流中快速干燥：

干燥过程：

1. 料液雾化成10-200 μm液滴
2. 液滴与热空气/辐射接触
3. 水分快速蒸发
4. 干粉收集

辐射辅助干燥：
红外辐射辅助喷雾干燥的优势：

• 干燥时间缩短30-50%
• 产品粒度更均匀
• 热效率提高
• 适用于热敏性物料

数学模型：
单液滴干燥方程：
$$\frac{dX}{dt}=-\frac{h_{m}A}{m_s}(p_{v,s}-p_{v,\infty })$$

$$m_d{c}_p\frac{d{T}_d}{dt}=hA(T_g-T_d)+\alpha AG-\dot{m}{L}_v$$

6.3 冷却塔

冷却塔利用水蒸发冷却循环水，辐射影响其性能：

工作原理：

1. 热水喷洒在填料上
2. 空气与水接触，部分水蒸发
3. 蒸发带走潜热，水温降低
4. 冷却水收集循环使用

辐射效应：

• 太阳辐射加热填料和水膜
• 增加水蒸发速率
• 降低冷却效率（夏季白天）
• 夜间辐射冷却增强性能

性能评估：
冷却效率：
$$\eta =\frac{T_{in}-T_{out}}{T_{in}-T_{wb}}$$

其中$T_{wb}$是湿球温度。

6.4 太阳能蒸馏

太阳能蒸馏利用太阳辐射蒸发海水，产生淡水：

工作原理：

1. 太阳辐射加热海水
2. 水蒸发，蒸汽上升
3. 蒸汽在冷凝面冷凝
4. 淡水收集

辐射-蒸发耦合：

• 辐射是唯一的能量输入
• 蒸发速率与辐射强度成正比
• 玻璃盖板透射太阳辐射但阻挡红外辐射（温室效应）

效率提升措施：

• 使用选择性吸收涂层
• 增加冷凝面积
• 优化倾角和保温
• 多级蒸馏设计

---

7. 数值方法与仿真技术

7.1 耦合问题的数值挑战

辐射-传质耦合问题的数值求解面临以下挑战：

多尺度特性：

• 辐射传输：光学厚度尺度
• 扩散传质：边界层尺度
• 对流：宏观流动尺度

非线性耦合：

• 辐射与温度四次方相关
• 蒸发与温度指数相关
• 物性参数温度依赖

计算复杂性：

• 辐射传输计算量大
• 瞬态过程需要小时间步
• 多物理场迭代收敛慢

7.2 算子分裂法

算子分裂法将耦合问题分解为子问题依次求解：

基本思想：
$$\frac{\partial \varphi }{\partial t}=L_1(\varphi )+L_2(\varphi )$$

分裂为：
$$\frac{\partial \varphi }{\partial t}=L_1(\varphi )\text{(步骤1)}$$
$$\frac{\partial \varphi }{\partial t}=L_2(\varphi )\text{(步骤2)}$$

辐射-传质分裂：

1. 传热步骤：求解能量方程（含辐射源项）
2. 传质步骤：求解质量扩散方程
3. 辐射步骤：求解辐射传输方程
4. 更新步骤：更新物性和源项

时间步进：

• 显式分裂：简单但可能不稳定
• 隐式分裂：稳定但计算量大
• 半隐式：折中方案

7.3 全耦合求解

对于强耦合问题，需要全耦合求解：

Newton-Raphson方法：
$$\mathbf{J}\delta \varphi =-\mathbf{R}$$

其中$\mathbf{J}$是Jacobian矩阵，$\mathbf{R}$是残差向量。

Picard迭代：
$$\mathbf{A}({\varphi }^n){\varphi }^{n+1}=\mathbf{b}({\varphi }^n)$$

收敛判据：
$$\Vert {\varphi }^{n+1}-{\varphi }^n\Vert <ϵ$$

7.4 辐射传输的数值方法

离散坐标法（DOM）：
将立体角离散为若干方向：
$${\int }_{4\pi }Id\Omega \approx \sum _{i=1}^N{w}_i{I}_i$$

有限体积法（FVM）：
在控制体积上积分辐射传输方程：
$$\sum _{faces}{I}_f{A}_f(\mathbf{s}\cdot \mathbf{n})=-\kappa IV+\kappa \sigma T^{4}V/\pi$$

蒙特卡洛法：
追踪大量光子/能束的随机行走：

• 发射：根据表面温度发射光子
• 传输：根据吸收/散射概率决定路径
• 吸收：记录能量沉积

7.5 验证与确认

数值模型的验证与确认至关重要：

验证（Verification）：

• 与解析解对比
• 网格无关性检验
• 时间步长收敛性
• 守恒性检查

确认（Validation）：

• 与实验数据对比
• 文献数据验证
• 敏感性分析
• 不确定性量化

---

8. Python仿真程序设计

8.1 程序架构

本主题的Python仿真程序包含以下模块：

1. 物理模型模块：蒸发模型、湿空气模型、多孔介质模型
2. 辐射计算模块：辐射传输计算、吸收系数模型
3. 数值求解模块：时间推进、空间离散、迭代求解
4. 可视化模块：温度分布、浓度分布、蒸发速率等

8.2 仿真案例设计

程序包含以下仿真案例：

案例1：单液滴蒸发

• 模拟液滴在辐射场中的蒸发过程
• 分析辐射强度对蒸发速率的影响
• 研究蒸发冷却效应

案例2：湿空气辐射传输

• 计算不同湿度条件下的辐射传输
• 分析水蒸气对辐射的吸收
• 研究辐射-扩散耦合

案例3：多孔介质干燥

• 模拟多孔物料的辐射干燥过程
• 分析温度分布和含水率变化
• 研究干燥速率曲线

案例4：喷雾冷却

• 模拟液滴群撞击热表面的冷却过程
• 分析辐射对冷却效率的影响
• 研究液滴动力学

案例5：太阳能蒸馏

• 模拟太阳能蒸馏器的热质传递
• 分析产水率和效率
• 优化设计参数

案例6：辐射-化学反应耦合

• 模拟光催化反应器
• 分析辐射分布对反应速率的影响
• 优化反应器设计

8.3 代码实现要点

数值稳定性：

• 使用隐式时间推进
• 添加人工扩散
• 限制时间步长

计算效率：

• 向量化运算
• 稀疏矩阵求解
• 自适应网格

物理准确性：

• 准确的物性参数
• 适当的边界条件
• 守恒性保证

---

附录：Python代码

完整的Python仿真程序见同目录下的 run_simulation.py 文件。程序包含上述所有案例的实现，可直接运行生成仿真结果。

# -*- coding: utf-8 -*-
"""
主题095：辐射换热与传质耦合现象
Radiative Heat Transfer and Mass Transfer Coupling
"""

import matplotlib
matplotlib.use('Agg')  # 非交互式后端
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle, Rectangle, FancyBboxPatch
from matplotlib.collections import PatchCollection
import warnings
warnings.filterwarnings('ignore')
import os

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

# 创建输出目录
output_dir = r'd:\\文档\\500仿真领域\\工程仿真\\辐射换热仿真\\主题095_辐射换热与传质耦合现象'
os.makedirs(output_dir, exist_ok=True)

print("=" * 70)
print("主题095：辐射换热与传质耦合现象")
print("=" * 70)

# 物理常数
sigma = 5.67e-8  # Stefan-Boltzmann常数 [W/m²K⁴]
R_universal = 8.314  # 通用气体常数 [J/mol·K]
M_w = 0.018  # 水分子摩尔质量 [kg/mol]
M_a = 0.029  # 干空气摩尔质量 [kg/mol]
L_v = 2.45e6  # 水的汽化潜热 [J/kg] (at 20°C)
cp_water = 4186  # 水比热容 [J/kg·K]
cp_air = 1005  # 空气比热容 [J/kg·K]
cp_vapor = 1860  # 水蒸气比热容 [J/kg·K]


# ==================== 案例1：单液滴蒸发 ====================
print("\n" + "=" * 70)
print("案例1：单液滴在辐射场中的蒸发")
print("=" * 70)

class DropletEvaporation:
    """液滴蒸发模型（考虑辐射效应）"""
    
    def __init__(self, d0, T0, T_ambient, RH):
        """
        初始化液滴参数
        d0: 初始直径 [m]
        T0: 初始温度 [K]
        T_ambient: 环境温度 [K]
        RH: 相对湿度 [0-1]
        """
        self.d0 = d0
        self.T = T0
        self.T_ambient = T_ambient
        self.RH = RH
        
        # 物性参数
        self.rho_l = 1000  # 水密度 [kg/m³]
        self.rho_air = 1.2  # 空气密度 [kg/m³]
        
        # 辐射特性（水在红外波段近似黑体）
        self.emissivity = 0.95
        self.absorptivity = 0.95
        
    def saturation_pressure(self, T):
        """计算饱和蒸汽压 [Pa] - Antoine方程"""
        if T < 273.15:
            T = 273.15
        # Antoine方程参数（水，1-100°C）
        A, B, C = 8.07131, 1730.63, 233.426
        # 转换为mmHg，再转换为Pa
        p_mmHg = 10 ** (A - B / (T - 273.15 + C))
        return p_mmHg * 133.322
    
    def vapor_pressure_ambient(self):
        """环境水蒸气分压 [Pa]"""
        return self.RH * self.saturation_pressure(self.T_ambient)
    
    def mass_transfer_coefficient(self, d, velocity=0.1):
        """计算传质系数 [m/s]"""
        # 扩散系数（水蒸气在空气中）
        D = 2.5e-5 * (self.T / 298.15) ** 1.5  # [m²/s]
        
        # Reynolds数
        nu = 1.5e-5  # 运动粘度 [m²/s]
        Re = velocity * d / nu
        
        # Schmidt数
        Sc = nu / D
        
        # Sherwood数（Ranz-Marshall关联式）
        Sh = 2 + 0.6 * Re**0.5 * Sc**(1/3)
        
        # 传质系数
        h_m = Sh * D / d
        return h_m
    
    def heat_transfer_coefficient(self, d, velocity=0.1):
        """计算传热系数 [W/m²K]"""
        k_air = 0.026  # 空气导热系数 [W/mK]
        
        # Reynolds数
        nu = 1.5e-5
        Re = velocity * d / nu
        
        # Prandtl数
        Pr = 0.71
        
        # Nusselt数
        Nu = 2 + 0.6 * Re**0.5 * Pr**(1/3)
        
        h = Nu * k_air / d
        return h
    
    def evaporate(self, t_max, dt, G_rad=0, velocity=0.1):
        """
        模拟蒸发过程
        t_max: 最大时间 [s]
        dt: 时间步长 [s]
        G_rad: 入射辐射强度 [W/m²]
        velocity: 相对速度 [m/s]
        """
        n_steps = int(t_max / dt) + 1
        t = np.linspace(0, t_max, n_steps)
        
        # 初始化数组
        d = np.zeros(n_steps)
        T = np.zeros(n_steps)
        m_dot = np.zeros(n_steps)
        
        d[0] = self.d0
        T[0] = self.T
        
        for i in range(n_steps - 1):
            # 当前状态
            d_i = d[i]
            T_i = T[i]
            
            if d_i < 1e-6:  # 液滴已完全蒸发
                d[i+1:] = 0
                T[i+1:] = self.T_ambient
                break
            
            # 表面积
            A = np.pi * d_i**2
            
            # 传质系数
            h_m = self.mass_transfer_coefficient(d_i, velocity)
            
            # 表面蒸汽压
            p_sat = self.saturation_pressure(T_i)
            rho_v_s = p_sat * M_w / (R_universal * T_i)
            
            # 环境蒸汽密度
            p_v_inf = self.vapor_pressure_ambient()
            rho_v_inf = p_v_inf * M_w / (R_universal * self.T_ambient)
            
            # 蒸发速率 [kg/s]
            m_dot_i = h_m * A * (rho_v_s - rho_v_inf)
            m_dot[i] = m_dot_i
            
            # 传热系数
            h = self.heat_transfer_coefficient(d_i, velocity)
            
            # 对流热流 [W]
            q_conv = h * A * (self.T_ambient - T_i)
            
            # 辐射热流 [W]（吸收 - 发射）
            q_rad = self.absorptivity * A * G_rad - self.emissivity * A * sigma * T_i**4
            
            # 液滴质量
            m = np.pi/6 * self.rho_l * d_i**3
            
            # 能量方程
            if m > 1e-12:
                dT_dt = (q_conv + q_rad - m_dot_i * L_v) / (m * cp_water)
            else:
                dT_dt = 0
            
            # 直径变化率（d²定律）
            if d_i > 1e-6:
                dd_dt = -2 * m_dot_i / (np.pi * self.rho_l * d_i**2)
            else:
                dd_dt = 0
            
            # 时间推进（欧拉法）
            T[i+1] = T_i + dT_dt * dt
            d[i+1] = d_i + dd_dt * dt
            
            # 限制温度
            T[i+1] = np.clip(T[i+1], 273.15, 373.15)
            
        return t, d, T, m_dot


# 案例1仿真
print("\n正在进行案例1仿真：单液滴蒸发...")

# 液滴参数
d0 = 100e-6  # 初始直径 100 μm
T0 = 293.15  # 初始温度 20°C
T_ambient = 303.15  # 环境温度 30°C
RH = 0.5  # 相对湿度 50%

# 不同辐射强度
t_max = 5.0  # 最大时间 [s]
dt = 0.001  # 时间步长 [s]
G_values = [0, 500, 1000, 2000]  # 辐射强度 [W/m²]

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

colors = ['blue', 'green', 'orange', 'red']
labels = ['No Radiation', 'G=500 W/m²', 'G=1000 W/m²', 'G=2000 W/m²']

for idx, G in enumerate(G_values):
    droplet = DropletEvaporation(d0, T0, T_ambient, RH)
    t, d, T, m_dot = droplet.evaporate(t_max, dt, G_rad=G)
    
    # 直径变化
    axes[0, 0].plot(t*1000, d*1e6, color=colors[idx], label=labels[idx], linewidth=2)
    
    # 温度变化
    axes[0, 1].plot(t*1000, T-273.15, color=colors[idx], label=labels[idx], linewidth=2)
    
    # 蒸发速率
    axes[1, 0].plot(t*1000, m_dot*1e6, color=colors[idx], label=labels[idx], linewidth=2)
    
    # 直径平方（d²定律验证）
    d_squared = (d/d0)**2
    axes[1, 1].plot(t*1000, d_squared, color=colors[idx], label=labels[idx], linewidth=2)

axes[0, 0].set_xlabel('Time (ms)', fontsize=11)
axes[0, 0].set_ylabel('Diameter (μm)', fontsize=11)
axes[0, 0].set_title('Droplet Diameter Evolution', fontsize=12, fontweight='bold')
axes[0, 0].legend(fontsize=9)
axes[0, 0].grid(True, alpha=0.3)

axes[0, 1].set_xlabel('Time (ms)', fontsize=11)
axes[0, 1].set_ylabel('Temperature (°C)', fontsize=11)
axes[0, 1].set_title('Droplet Temperature Evolution', fontsize=12, fontweight='bold')
axes[0, 1].legend(fontsize=9)
axes[0, 1].grid(True, alpha=0.3)

axes[1, 0].set_xlabel('Time (ms)', fontsize=11)
axes[1, 0].set_ylabel('Evaporation Rate (mg/s)', fontsize=11)
axes[1, 0].set_title('Evaporation Rate vs Time', fontsize=12, fontweight='bold')
axes[1, 0].legend(fontsize=9)
axes[1, 0].grid(True, alpha=0.3)

axes[1, 1].set_xlabel('Time (ms)', fontsize=11)
axes[1, 1].set_ylabel('(d/d₀)²', fontsize=11)
axes[1, 1].set_title('d² Law Verification', fontsize=12, fontweight='bold')
axes[1, 1].legend(fontsize=9)
axes[1, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/case1_droplet_evaporation.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 案例1图像已保存: case1_droplet_evaporation.png")


# ==================== 案例2：湿空气辐射传输 ====================
print("\n" + "=" * 70)
print("案例2：湿空气中的辐射传输")
print("=" * 70)

class HumidAirRadiation:
    """湿空气辐射传输模型"""
    
    def __init__(self, T, P=101325):
        """
        T: 温度 [K]
        P: 压力 [Pa]
        """
        self.T = T
        self.P = P
        
    def water_vapor_absorption(self, wavelength, RH):
        """
        计算水蒸气吸收系数 [1/m]
        使用简化模型
        """
        # 饱和蒸汽压
        p_sat = 611.2 * np.exp(17.67 * (self.T - 273.15) / (self.T - 29.65))
        p_vapor = RH * p_sat
        
        # 水蒸气密度
        rho_v = p_vapor * M_w / (R_universal * self.T)
        
        # 吸收系数模型（简化）
        # 主要吸收带：2.7 μm, 6.3 μm, 转动带(>15 μm)
        lambda_um = wavelength * 1e6  # 转换为μm
        
        # 基础吸收系数
        kappa_base = 0.1 * rho_v  # [1/m]
        
        # 吸收带增强
        if 2.5 <= lambda_um <= 3.0:
            kappa = kappa_base * 50  # 2.7 μm带
        elif 5.5 <= lambda_um <= 7.5:
            kappa = kappa_base * 100  # 6.3 μm带
        elif lambda_um >= 15:
            kappa = kappa_base * 30  # 转动带
        else:
            kappa = kappa_base * 0.1  # 弱吸收区
        
        return kappa
    
    def radiative_transfer(self, x, I0, wavelength, RH, T_profile=None):
        """
        计算辐射强度分布（Beer-Lambert定律）
        x: 位置数组 [m]
        I0: 入射强度 [W/m²]
        wavelength: 波长 [m]
        RH: 相对湿度
        T_profile: 温度分布（可选）
        """
        kappa = self.water_vapor_absorption(wavelength, RH)
        
        # 辐射衰减
        I = I0 * np.exp(-kappa * x)
        
        # 如果提供温度分布，添加发射项
        if T_profile is not None:
            # 发射项（简化处理）
            emission = sigma * T_profile**4 / np.pi
            # 源项积分（简化）
            for i in range(1, len(x)):
                dx = x[i] - x[i-1]
                I[i] = I[i-1] * np.exp(-kappa * dx) + \
                       kappa * emission[i] * (1 - np.exp(-kappa * dx))
        
        return I
    
    def greenhouse_effect(self, T_surface, T_atm, RH, atm_thickness=10000):
        """
        计算温室效应
        """
        # 地表发射的长波辐射
        I_surface = sigma * T_surface**4
        
        # 大气吸收率（简化）
        p_sat = 611.2 * np.exp(17.67 * (T_atm - 273.15) / (T_atm - 29.65))
        p_vapor = RH * p_sat
        
        # 光学厚度（简化）
        tau = 0.1 * p_vapor / 101325 * atm_thickness / 1000
        
        # 大气吸收
        absorbance = 1 - np.exp(-tau)
        
        # 大气发射回地表
        I_back = absorbance * sigma * T_atm**4
        
        # 净辐射
        I_net = I_surface - I_back
        
        return {
            'I_surface': I_surface,
            'I_back': I_back,
            'I_net': I_net,
            'absorbance': absorbance,
            'tau': tau
        }


# 案例2仿真
print("\n正在进行案例2仿真：湿空气辐射传输...")

humid_air = HumidAirRadiation(T=300)

# 创建图形
fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# 1. 不同湿度下的辐射衰减
x = np.linspace(0, 100, 100)  # 0-100 m
I0 = 1000  # 入射强度 [W/m²]
wavelength = 10e-6  # 10 μm (红外)

RH_values = [0.2, 0.4, 0.6, 0.8]
colors = ['blue', 'green', 'orange', 'red']

for idx, RH in enumerate(RH_values):
    I = humid_air.radiative_transfer(x, I0, wavelength, RH)
    axes[0, 0].plot(x, I, color=colors[idx], label=f'RH={RH*100:.0f}%', linewidth=2)

axes[0, 0].set_xlabel('Distance (m)', fontsize=11)
axes[0, 0].set_ylabel('Radiation Intensity (W/m²)', fontsize=11)
axes[0, 0].set_title('Radiation Attenuation in Humid Air (λ=10μm)', fontsize=12, fontweight='bold')
axes[0, 0].legend(fontsize=9)
axes[0, 0].grid(True, alpha=0.3)

# 2. 不同波长下的吸收
wavelengths = np.linspace(1, 20, 100) * 1e-6  # 1-20 μm
RH = 0.6

kappa_values = []
for wl in wavelengths:
    kappa = humid_air.water_vapor_absorption(wl, RH)
    kappa_values.append(kappa)

axes[0, 1].semilogy(wavelengths*1e6, kappa_values, 'b-', linewidth=2)
axes[0, 1].axvline(x=2.7, color='r', linestyle='--', label='2.7 μm band', alpha=0.7)
axes[0, 1].axvline(x=6.3, color='g', linestyle='--', label='6.3 μm band', alpha=0.7)
axes[0, 1].axvline(x=15, color='m', linestyle='--', label='Rotational band', alpha=0.7)
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('Water Vapor Absorption Spectrum (RH=60%)', fontsize=12, fontweight='bold')
axes[0, 1].legend(fontsize=9)
axes[0, 1].grid(True, alpha=0.3)

# 3. 温室效应分析
T_surface = 288  # 地表温度 15°C
T_atm = 255  # 大气有效温度 -18°C
RH_range = np.linspace(0.1, 0.9, 50)

I_back_values = []
absorbance_values = []

for RH in RH_range:
    result = humid_air.greenhouse_effect(T_surface, T_atm, RH)
    I_back_values.append(result['I_back'])
    absorbance_values.append(result['absorbance'])

axes[1, 0].plot(RH_range*100, np.array(I_back_values), 'r-', linewidth=2)
axes[1, 0].axhline(y=sigma*T_surface**4, color='b', linestyle='--', label='Surface emission', alpha=0.7)
axes[1, 0].set_xlabel('Relative Humidity (%)', fontsize=11)
axes[1, 0].set_ylabel('Back Radiation (W/m²)', fontsize=11)
axes[1, 0].set_title('Greenhouse Effect vs Humidity', fontsize=12, fontweight='bold')
axes[1, 0].legend(fontsize=9)
axes[1, 0].grid(True, alpha=0.3)

# 4. 穿透深度
penetration_depths = []
for RH in RH_range:
    kappa_avg = np.mean([humid_air.water_vapor_absorption(wl * 1e-6, RH) 
                         for wl in [5, 10, 15]])
    penetration_depths.append(1 / (kappa_avg + 1e-10))

axes[1, 1].plot(RH_range*100, np.array(penetration_depths), 'g-', linewidth=2)
axes[1, 1].set_xlabel('Relative Humidity (%)', fontsize=11)
axes[1, 1].set_ylabel('Penetration Depth (m)', fontsize=11)
axes[1, 1].set_title('Radiation Penetration Depth', fontsize=12, fontweight='bold')
axes[1, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/case2_humid_air_radiation.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 案例2图像已保存: case2_humid_air_radiation.png")


# ==================== 案例3：多孔介质干燥 ====================
print("\n" + "=" * 70)
print("案例3：多孔介质辐射干燥")
print("=" * 70)

class PorousMediaDrying:
    """多孔介质辐射干燥模型"""
    
    def __init__(self, L, n_points, porosity=0.5):
        """
        L: 物料厚度 [m]
        n_points: 网格点数
        porosity: 孔隙率
        """
        self.L = L
        self.nx = n_points
        self.porosity = porosity
        self.x = np.linspace(0, L, n_points)
        self.dx = L / (n_points - 1)
        
        # 物性参数
        self.rho_s = 1500  # 固体密度 [kg/m³]
        self.cp_s = 1500  # 固体比热 [J/kgK]
        self.k_s = 0.5  # 固体导热系数 [W/mK]
        
        # 有效物性
        self.rho_eff = (1 - porosity) * self.rho_s
        self.cp_eff = (1 - porosity) * self.cp_s
        self.k_eff = (1 - porosity) * self.k_s
        
        # 辐射特性
        self.kappa = 500  # 吸收系数 [1/m]
        
    def simulate(self, t_max, dt, T_initial, T_air, RH, G_rad, h_conv=10):
        """
        模拟干燥过程
        t_max: 最大时间 [s]
        dt: 时间步长 [s]
        T_initial: 初始温度 [K]
        T_air: 空气温度 [K]
        RH: 相对湿度
        G_rad: 辐射强度 [W/m²]
        h_conv: 对流传热系数 [W/m²K]
        """
        n_steps = int(t_max / dt) + 1
        
        # 初始化场变量
        T = np.ones(self.nx) * T_initial
        X = np.ones(self.nx) * 0.5  # 初始含水率 50% (干基)
        
        # 记录历史
        T_history = [T.copy()]
        X_history = [X.copy()]
        drying_rate_history = []
        
        for step in range(n_steps - 1):
            T_old = T.copy()
            X_old = X.copy()
            
            # 1. 计算辐射热源项
            q_rad = G_rad * self.kappa * np.exp(-self.kappa * self.x)
            
            # 2. 能量方程（隐式求解）
            # 构建三对角矩阵
            alpha = self.k_eff / self.rho_eff / self.cp_eff * dt / self.dx**2
            
            # 内部点
            for i in range(1, self.nx - 1):
                T[i] = T_old[i] + alpha * (T_old[i+1] - 2*T_old[i] + T_old[i-1]) + \
                       q_rad[i] * dt / (self.rho_eff * self.cp_eff)
            
            # 边界条件
            # 表面（x=0）：对流+辐射
            T[0] = T_old[0] + 2*alpha*(T_old[1] - T_old[0]) + \
                   (h_conv * (T_air - T_old[0]) + G_rad) * dt / (self.rho_eff * self.cp_eff)
            
            # 底部（x=L）：绝热
            T[-1] = T_old[-1] + 2*alpha*(T_old[-2] - T_old[-1])
            
            # 3. 传质方程（简化模型）
            # 干燥速率与温度和含水率相关
            D_eff = 1e-9 * np.exp(0.05 * (T - 300))  # 有效扩散系数 [m²/s]
            
            # 表面蒸发
            p_sat_surface = 611.2 * np.exp(17.67 * (T[0] - 273.15) / (T[0] - 29.65))
            p_vapor_air = RH * 611.2 * np.exp(17.67 * (T_air - 273.15) / (T_air - 29.65))
            
            # 表面传质
            h_m = 0.01  # 传质系数 [m/s]
            m_dot_surface = h_m * (p_sat_surface - p_vapor_air) * M_w / (R_universal * T[0])
            
            # 内部水分迁移（扩散）
            for i in range(1, self.nx - 1):
                dX_dt = D_eff[i] * (X_old[i+1] - 2*X_old[i] + X_old[i-1]) / self.dx**2
                X[i] = X_old[i] + dX_dt * dt
            
            # 表面含水率变化
            X[0] = X_old[0] - m_dot_surface * dt / (self.rho_s * self.dx)
            
            # 底部边界
            X[-1] = X_old[-1]
            
            # 限制含水率
            X = np.clip(X, 0.05, 1.0)
            
            # 记录
            if step % 100 == 0:
                T_history.append(T.copy())
                X_history.append(X.copy())
                drying_rate_history.append(m_dot_surface)
        
        return np.array(T_history), np.array(X_history), np.array(drying_rate_history)


# 案例3仿真
print("\n正在进行案例3仿真：多孔介质干燥...")

# 创建模型
porous = PorousMediaDrying(L=0.02, n_points=50, porosity=0.6)

# 仿真参数
t_max = 3600  # 1小时
dt = 1.0  # 1秒

# 两种工况：对流干燥 vs 辐射干燥
T_history_conv, X_history_conv, rate_conv = porous.simulate(
    t_max, dt, T_initial=293, T_air=333, RH=0.3, G_rad=0, h_conv=15)

T_history_rad, X_history_rad, rate_rad = porous.simulate(
    t_max, dt, T_initial=293, T_air=333, RH=0.3, G_rad=2000, h_conv=15)

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

# 1. 温度分布演变
times = np.linspace(0, t_max/60, len(T_history_conv))  # 分钟

# 选择几个时刻
time_indices = [0, len(T_history_conv)//4, len(T_history_conv)//2, -1]
for idx in time_indices:
    axes[0, 0].plot(porous.x*1000, T_history_conv[idx]-273.15, 
                    '--', alpha=0.7, label=f'Conv: {times[idx]:.0f} min')
    axes[0, 0].plot(porous.x*1000, T_history_rad[idx]-273.15, 
                    '-', linewidth=2, label=f'Rad: {times[idx]:.0f} min')

axes[0, 0].set_xlabel('Depth (mm)', fontsize=11)
axes[0, 0].set_ylabel('Temperature (°C)', fontsize=11)
axes[0, 0].set_title('Temperature Distribution Evolution', fontsize=12, fontweight='bold')
axes[0, 0].legend(fontsize=8, ncol=2)
axes[0, 0].grid(True, alpha=0.3)

# 2. 含水率分布演变
for idx in time_indices:
    axes[0, 1].plot(porous.x*1000, X_history_conv[idx]*100, 
                    '--', alpha=0.7, label=f'Conv: {times[idx]:.0f} min')
    axes[0, 1].plot(porous.x*1000, X_history_rad[idx]*100, 
                    '-', linewidth=2, label=f'Rad: {times[idx]:.0f} min')

axes[0, 1].set_xlabel('Depth (mm)', fontsize=11)
axes[0, 1].set_ylabel('Moisture Content (%)', fontsize=11)
axes[0, 1].set_title('Moisture Distribution Evolution', fontsize=12, fontweight='bold')
axes[0, 1].legend(fontsize=8, ncol=2)
axes[0, 1].grid(True, alpha=0.3)

# 3. 干燥速率曲线
rate_time = np.linspace(0, t_max/60, len(rate_conv))
axes[1, 0].plot(rate_time, rate_conv*1e3, 'b--', linewidth=2, label='Convection')
axes[1, 0].plot(rate_time[:len(rate_rad)], rate_rad*1e3, 'r-', linewidth=2, label='Radiation')
axes[1, 0].set_xlabel('Time (min)', fontsize=11)
axes[1, 0].set_ylabel('Drying Rate (g/m²s)', fontsize=11)
axes[1, 0].set_title('Drying Rate vs Time', fontsize=12, fontweight='bold')
axes[1, 0].legend(fontsize=9)
axes[1, 0].grid(True, alpha=0.3)

# 4. 平均含水率变化
avg_X_conv = np.mean(X_history_conv, axis=1) * 100
avg_X_rad = np.mean(X_history_rad, axis=1) * 100

# 确保数组长度一致
min_len = min(len(rate_time), len(avg_X_conv), len(avg_X_rad))
axes[1, 1].plot(rate_time[:min_len], avg_X_conv[:min_len], 'b--', linewidth=2, label='Convection')
axes[1, 1].plot(rate_time[:min_len], avg_X_rad[:min_len], 'r-', linewidth=2, label='Radiation')
axes[1, 1].set_xlabel('Time (min)', fontsize=11)
axes[1, 1].set_ylabel('Average Moisture Content (%)', fontsize=11)
axes[1, 1].set_title('Average Moisture Content vs Time', fontsize=12, fontweight='bold')
axes[1, 1].legend(fontsize=9)
axes[1, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/case3_porous_drying.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 案例3图像已保存: case3_porous_drying.png")


# ==================== 案例4：喷雾冷却 ====================
print("\n" + "=" * 70)
print("案例4：喷雾冷却（辐射增强效应）")
print("=" * 70)

class SprayCooling:
    """喷雾冷却模型"""
    
    def __init__(self, n_droplets, surface_area):
        """
        n_droplets: 液滴数量
        surface_area: 表面积 [m²]
        """
        self.n_droplets = n_droplets
        self.surface_area = surface_area
        
    def simulate_cooling(self, T_surface_initial, T_liquid, G_rad, 
                         droplet_diameter, impact_velocity, t_max, dt):
        """
        模拟喷雾冷却过程
        """
        n_steps = int(t_max / dt) + 1
        t = np.linspace(0, t_max, n_steps)
        
        # 表面温度
        T_surface = np.zeros(n_steps)
        T_surface[0] = T_surface_initial
        
        # 材料热物性（假设为不锈钢）
        rho_s = 8000  # [kg/m³]
        cp_s = 500  # [J/kgK]
        k_s = 16  # [W/mK]
        thickness = 0.005  # 有效厚度 [m]
        
        # 有效热容
        C_eff = rho_s * cp_s * thickness * self.surface_area
        
        # 液滴参数
        V_droplet = np.pi/6 * droplet_diameter**3
        m_droplet = 1000 * V_droplet
        
        for i in range(n_steps - 1):
            # 1. 辐射热流
            q_rad = G_rad * self.surface_area
            
            # 2. 液滴蒸发冷却
            # 简化模型：假设每个液滴带走一部分热量
            # 传热系数
            h_droplet = 10000  # 液滴撞击传热系数 [W/m²K]
            A_droplet = np.pi * droplet_diameter**2
            
            # 每个液滴带走的热量
            q_per_droplet = h_droplet * A_droplet * (T_surface[i] - T_liquid) * dt
            
            # 总冷却功率
            n_impact = self.n_droplets * dt  # 单位时间撞击的液滴数（简化）
            q_cooling = n_impact * q_per_droplet
            
            # 蒸发潜热（假设部分蒸发）
            evap_fraction = 0.3
            q_evap = n_impact * m_droplet * evap_fraction * L_v
            
            # 能量平衡
            dT = (q_rad - q_cooling - q_evap) * dt / C_eff
            T_surface[i+1] = T_surface[i] + dT
            
            # 限制温度
            if T_surface[i+1] < T_liquid:
                T_surface[i+1] = T_liquid
        
        return t, T_surface
    
    def efficiency_analysis(self, T_surface_range, T_liquid, G_rad, 
                           droplet_diameter, impact_velocity):
        """
        分析冷却效率
        """
        efficiencies = []
        heat_fluxes = []
        
        for T_s in T_surface_range:
            # 最大可能冷却（所有热量用于蒸发）
            q_max = G_rad
            
            # 实际冷却
            h_droplet = 10000
            A_droplet = np.pi * droplet_diameter**2
            q_cooling = h_droplet * A_droplet * self.n_droplets * (T_s - T_liquid)
            q_evap = self.n_droplets * (np.pi/6 * droplet_diameter**3 * 1000) * L_v
            
            q_total = q_cooling + q_evap
            
            # 效率
            eta = min(q_total / (q_max + 1e-10), 1.0)
            efficiencies.append(eta)
            heat_fluxes.append(q_total / self.surface_area / 1e6)  # MW/m²
        
        return np.array(efficiencies), np.array(heat_fluxes)


# 案例4仿真
print("\n正在进行案例4仿真：喷雾冷却...")

spray = SprayCooling(n_droplets=10000, surface_area=0.01)

# 仿真参数
t_max = 10.0  # 10秒
dt = 0.001
droplet_diameter = 100e-6  # 100 μm
impact_velocity = 20  # m/s

# 不同辐射强度
G_values = [0, 5000, 10000, 20000]  # W/m²

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

colors = ['blue', 'green', 'orange', 'red']
labels = ['No Radiation', 'G=5 kW/m²', 'G=10 kW/m²', 'G=20 kW/m²']

for idx, G in enumerate(G_values):
    t, T_surface = spray.simulate_cooling(
        T_surface_initial=500,  # 227°C
        T_liquid=293,  # 20°C
        G_rad=G,
        droplet_diameter=droplet_diameter,
        impact_velocity=impact_velocity,
        t_max=t_max,
        dt=dt
    )
    
    axes[0, 0].plot(t*1000, T_surface-273.15, color=colors[idx], 
                    label=labels[idx], linewidth=2)

axes[0, 0].set_xlabel('Time (ms)', fontsize=11)
axes[0, 0].set_ylabel('Surface Temperature (°C)', fontsize=11)
axes[0, 0].set_title('Surface Cooling with Spray', fontsize=12, fontweight='bold')
axes[0, 0].legend(fontsize=9)
axes[0, 0].grid(True, alpha=0.3)

# 效率分析
T_range = np.linspace(300, 600, 50)
efficiencies, heat_fluxes = spray.efficiency_analysis(
    T_range, T_liquid=293, G_rad=10000,
    droplet_diameter=droplet_diameter, impact_velocity=impact_velocity)

axes[0, 1].plot(T_range-273.15, efficiencies*100, 'b-', linewidth=2)
axes[0, 1].set_xlabel('Surface Temperature (°C)', fontsize=11)
axes[0, 1].set_ylabel('Cooling Efficiency (%)', fontsize=11)
axes[0, 1].set_title('Spray Cooling Efficiency', fontsize=12, fontweight='bold')
axes[0, 1].grid(True, alpha=0.3)

# 热流密度
axes[1, 0].plot(T_range-273.15, heat_fluxes, 'r-', linewidth=2)
axes[1, 0].set_xlabel('Surface Temperature (°C)', fontsize=11)
axes[1, 0].set_ylabel('Heat Flux (MW/m²)', fontsize=11)
axes[1, 0].set_title('Heat Removal Rate', fontsize=12, fontweight='bold')
axes[1, 0].grid(True, alpha=0.3)

# 液滴直径影响
diameters = np.linspace(20, 200, 50) * 1e-6
cooling_rates = []

for d in diameters:
    t, T_s = spray.simulate_cooling(
        T_surface_initial=500, T_liquid=293, G_rad=10000,
        droplet_diameter=d, impact_velocity=impact_velocity,
        t_max=0.1, dt=0.0001
    )
    cooling_rate = (T_s[0] - T_s[-1]) / 0.1
    cooling_rates.append(cooling_rate)

axes[1, 1].plot(diameters*1e6, cooling_rates, 'g-', linewidth=2)
axes[1, 1].set_xlabel('Droplet Diameter (μm)', fontsize=11)
axes[1, 1].set_ylabel('Cooling Rate (K/s)', fontsize=11)
axes[1, 1].set_title('Effect of Droplet Size', fontsize=12, fontweight='bold')
axes[1, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/case4_spray_cooling.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 案例4图像已保存: case4_spray_cooling.png")


# ==================== 案例5：太阳能蒸馏 ====================
print("\n" + "=" * 70)
print("案例5：太阳能蒸馏器")
print("=" * 70)

class SolarStill:
    """太阳能蒸馏器模型"""
    
    def __init__(self, area, tilt_angle, glass_thickness=0.004):
        """
        area: 蒸馏器面积 [m²]
        tilt_angle: 倾角 [度]
        glass_thickness: 玻璃厚度 [m]
        """
        self.area = area
        self.tilt_angle = np.radians(tilt_angle)
        self.glass_thickness = glass_thickness
        
        # 玻璃特性
        self.tau_glass = 0.85  # 透射率
        self.alpha_water = 0.95  # 水吸收率
        self.epsilon_water = 0.95  # 水发射率
        
    def daily_simulation(self, day_of_year, latitude, G_max, T_ambient_func, RH_func):
        """
        日仿真
        day_of_year: 一年中的第几天
        latitude: 纬度 [度]
        G_max: 最大太阳辐射 [W/m²]
        T_ambient_func: 环境温度函数
        RH_func: 相对湿度函数
        """
        # 时间数组（每小时）
        hours = np.linspace(0, 24, 100)
        
        # 太阳高度角（简化）
        declination = 23.45 * np.sin(np.radians(360 * (284 + day_of_year) / 365))
        hour_angle = np.radians(15 * (hours - 12))
        
        sin_altitude = np.sin(np.radians(latitude)) * np.sin(np.radians(declination)) + \
                       np.cos(np.radians(latitude)) * np.cos(np.radians(declination)) * np.cos(hour_angle)
        sin_altitude = np.maximum(sin_altitude, 0)
        
        # 入射太阳辐射（简化模型）
        G_solar = G_max * sin_altitude * (sin_altitude > 0)
        
        # 初始化
        T_water = np.zeros_like(hours)
        T_glass = np.zeros_like(hours)
        evaporation_rate = np.zeros_like(hours)
        cumulative_yield = np.zeros_like(hours)
        
        T_water[0] = T_ambient_func(0) + 5
        T_glass[0] = T_ambient_func(0)
        
        dt = (hours[1] - hours[0]) * 3600  # 秒
        
        # 热容
        m_water = 10  # 水量 [kg]
        C_water = m_water * cp_water
        C_glass = 2500 * 840 * self.glass_thickness * self.area
        
        for i in range(len(hours) - 1):
            T_amb = T_ambient_func(hours[i])
            RH = RH_func(hours[i])
            
            # 吸收的太阳能
            Q_absorbed = G_solar[i] * self.area * self.tau_glass * self.alpha_water
            
            # 水向玻璃的传热
            h_wg = 10  # 水-玻璃传热系数 [W/m²K]
            Q_wg = h_wg * self.area * (T_water[i] - T_glass[i])
            
            # 蒸发传热
            p_sat_water = 611.2 * np.exp(17.67 * (T_water[i] - 273.15) / (T_water[i] - 29.65))
            p_sat_glass = 611.2 * np.exp(17.67 * (T_glass[i] - 273.15) / (T_glass[i] - 29.65))
            
            h_m = 0.02  # 传质系数 [m/s]
            m_dot_evap = h_m * self.area * (p_sat_water - p_sat_glass) * M_w / (R_universal * T_water[i])
            Q_evap = m_dot_evap * L_v
            
            # 玻璃向环境的散热
            h_ga = 5 + 3 * np.random.random()  # 对流+辐射 [W/m²K]
            Q_ga = h_ga * self.area * (T_glass[i] - T_amb)
            
            # 能量平衡
            dT_water = (Q_absorbed - Q_wg - Q_evap) * dt / C_water
            dT_glass = (Q_wg + Q_evap - Q_ga) * dt / C_glass
            
            T_water[i+1] = T_water[i] + dT_water
            T_glass[i+1] = T_glass[i] + dT_glass
            
            evaporation_rate[i+1] = m_dot_evap * 3600  # kg/h
            
            if i == 0:
                cumulative_yield[i+1] = m_dot_evap * dt
            else:
                cumulative_yield[i+1] = cumulative_yield[i] + m_dot_evap * dt
        
        return {
            'hours': hours,
            'G_solar': G_solar,
            'T_water': T_water,
            'T_glass': T_glass,
            'evaporation_rate': evaporation_rate,
            'cumulative_yield': cumulative_yield * 1000  # 转换为g
        }


# 案例5仿真
print("\n正在进行案例5仿真：太阳能蒸馏...")

# 创建蒸馏器
still = SolarStill(area=1.0, tilt_angle=15)

# 环境条件函数
def T_ambient(hour):
    return 298 + 5 * np.sin(np.radians(15 * (hour - 6))) if 6 <= hour <= 18 else 293

def RH_ambient(hour):
    return 0.6 + 0.1 * np.sin(np.radians(15 * (hour - 12)))

# 不同季节的仿真
seasons = [
    (80, 30, 'Spring'),   # 春分
    (172, 35, 'Summer'),  # 夏至
    (266, 25, 'Autumn'),  # 秋分
    (355, 15, 'Winter')   # 冬至
]

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

colors = ['green', 'red', 'orange', 'blue']

for idx, (day, G_max, season) in enumerate(seasons):
    result = still.daily_simulation(day, latitude=35, G_max=G_max*1.2,
                                   T_ambient_func=T_ambient, RH_func=RH_ambient)
    
    axes[0, 0].plot(result['hours'], result['G_solar'], 
                    color=colors[idx], label=season, linewidth=2)
    axes[0, 1].plot(result['hours'], result['T_water']-273.15, 
                    color=colors[idx], label=season, linewidth=2)
    axes[1, 0].plot(result['hours'], result['evaporation_rate']*1000, 
                    color=colors[idx], label=season, linewidth=2)
    axes[1, 1].plot(result['hours'], result['cumulative_yield'], 
                    color=colors[idx], label=season, linewidth=2)

axes[0, 0].set_xlabel('Hour of Day', fontsize=11)
axes[0, 0].set_ylabel('Solar Radiation (W/m²)', fontsize=11)
axes[0, 0].set_title('Solar Radiation Profile', fontsize=12, fontweight='bold')
axes[0, 0].legend(fontsize=9)
axes[0, 0].grid(True, alpha=0.3)
axes[0, 0].set_xlim(0, 24)

axes[0, 1].set_xlabel('Hour of Day', fontsize=11)
axes[0, 1].set_ylabel('Water Temperature (°C)', fontsize=11)
axes[0, 1].set_title('Water Temperature Evolution', fontsize=12, fontweight='bold')
axes[0, 1].legend(fontsize=9)
axes[0, 1].grid(True, alpha=0.3)
axes[0, 1].set_xlim(0, 24)

axes[1, 0].set_xlabel('Hour of Day', fontsize=11)
axes[1, 0].set_ylabel('Evaporation Rate (g/h)', fontsize=11)
axes[1, 0].set_title('Evaporation Rate', fontsize=12, fontweight='bold')
axes[1, 0].legend(fontsize=9)
axes[1, 0].grid(True, alpha=0.3)
axes[1, 0].set_xlim(0, 24)

axes[1, 1].set_xlabel('Hour of Day', fontsize=11)
axes[1, 1].set_ylabel('Cumulative Yield (g)', fontsize=11)
axes[1, 1].set_title('Cumulative Fresh Water Yield', fontsize=12, fontweight='bold')
axes[1, 1].legend(fontsize=9)
axes[1, 1].grid(True, alpha=0.3)
axes[1, 1].set_xlim(0, 24)

plt.tight_layout()
plt.savefig(f'{output_dir}/case5_solar_still.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 案例5图像已保存: case5_solar_still.png")


# ==================== 案例6：辐射-化学反应耦合 ====================
print("\n" + "=" * 70)
print("案例6：光催化反应器中的辐射-反应耦合")
print("=" * 70)

class PhotocatalyticReactor:
    """光催化反应器模型"""
    
    def __init__(self, length, width, height, n_points=50):
        """
        length: 反应器长度 [m]
        width: 反应器宽度 [m]
        height: 反应器高度 [m]
        """
        self.length = length
        self.width = width
        self.height = height
        self.nx = n_points
        self.x = np.linspace(0, length, n_points)
        self.dx = length / (n_points - 1)
        
    def radiation_distribution(self, I0, kappa, scattering=False):
        """
        计算反应器内辐射分布
        I0: 入射辐射强度 [W/m²]
        kappa: 吸收系数 [1/m]
        """
        if scattering:
            # 考虑散射（简化扩散近似）
            D_diffusion = 1 / (3 * kappa)
            # 解析解
            I = I0 * np.exp(-np.sqrt(kappa / D_diffusion) * self.x)
        else:
            # 纯吸收（Beer-Lambert定律）
            I = I0 * np.exp(-kappa * self.x)
        
        return I
    
    def reaction_rate(self, I, C, T, quantum_yield=0.1, k_thermal=0.01):
        """
        计算反应速率
        I: 辐射强度 [W/m²]
        C: 反应物浓度 [mol/m³]
        T: 温度 [K]
        """
        # 光化学反应速率
        # 假设：反应速率正比于吸收的光子数和反应物浓度
        r_photo = quantum_yield * kappa * I * C / (100 + I)  # 饱和模型
        
        # 热化学反应速率（Arrhenius）
        Ea = 50000  # 活化能 [J/mol]
        A = 1e5  # 指前因子
        r_thermal = A * np.exp(-Ea / (R_universal * T)) * C
        
        return r_photo, r_thermal
    
    def simulate(self, I0, C0, T0, flow_velocity, t_max, dt, kappa=10):
        """
        模拟反应器性能
        """
        n_steps = int(t_max / dt) + 1
        
        # 初始化
        C = np.ones(self.nx) * C0
        T = np.ones(self.nx) * T0
        
        C_history = [C.copy()]
        T_history = [T.copy()]
        conversion_history = []
        
        for step in range(n_steps - 1):
            # 1. 计算辐射分布
            I = self.radiation_distribution(I0, kappa)
            
            # 2. 计算反应速率
            r_photo = np.zeros(self.nx)
            r_thermal = np.zeros(self.nx)
            
            for i in range(self.nx):
                rp, rt = self.reaction_rate(I[i], C[i], T[i])
                r_photo[i] = rp
                r_thermal[i] = rt
            
            r_total = r_photo + r_thermal
            
            # 3. 浓度方程（对流+反应）
            C_old = C.copy()
            for i in range(1, self.nx):
                # 对流项（迎风格式）
                adv = -flow_velocity * (C_old[i] - C_old[i-1]) / self.dx
                # 反应项
                react = -r_total[i]
                
                C[i] = C_old[i] + (adv + react) * dt
            
            C[0] = C0  # 入口边界
            C = np.maximum(C, 0)
            
            # 4. 能量方程（对流+反应热+辐射加热）
            T_old = T.copy()
            rho_mix = 1.2  # 混合物密度 [kg/m³]
            cp_mix = 1200  # 混合物比热 [J/kgK]
            DeltaH = -50000  # 反应热 [J/mol]
            
            for i in range(1, self.nx):
                adv_T = -flow_velocity * (T_old[i] - T_old[i-1]) / self.dx
                react_heat = -r_total[i] * DeltaH / (rho_mix * cp_mix)
                rad_heat = kappa * I[i] / (rho_mix * cp_mix)
                
                T[i] = T_old[i] + (adv_T + react_heat + rad_heat) * dt
            
            T[0] = T0
            
            # 5. 记录
            if step % 100 == 0:
                C_history.append(C.copy())
                T_history.append(T.copy())
                conversion = 1 - C[-1] / C0
                conversion_history.append(conversion)
        
        return {
            'x': self.x,
            'I': I,
            'C_final': C,
            'T_final': T,
            'C_history': np.array(C_history),
            'T_history': np.array(T_history),
            'conversion': np.array(conversion_history),
            'r_photo': r_photo,
            'r_thermal': r_thermal
        }


# 案例6仿真
print("\n正在进行案例6仿真：光催化反应器...")

reactor = PhotocatalyticReactor(length=0.5, width=0.1, height=0.05, n_points=50)

# 基础仿真
result = reactor.simulate(
    I0=1000,  # W/m²
    C0=10,    # mol/m³
    T0=298,   # K
    flow_velocity=0.01,  # m/s
    t_max=100,  # s
    dt=0.01,
    kappa=5
)

# 不同入射辐射强度的影响
I0_values = [100, 500, 1000, 2000]

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

colors = ['blue', 'green', 'orange', 'red']

# 1. 辐射分布
for idx, I0 in enumerate(I0_values):
    I = reactor.radiation_distribution(I0, kappa=5)
    axes[0, 0].plot(reactor.x*100, I, color=colors[idx], 
                    label=f'I₀={I0} W/m²', linewidth=2)

axes[0, 0].set_xlabel('Position (cm)', fontsize=11)
axes[0, 0].set_ylabel('Radiation Intensity (W/m²)', fontsize=11)
axes[0, 0].set_title('Radiation Distribution in Reactor', fontsize=12, fontweight='bold')
axes[0, 0].legend(fontsize=9)
axes[0, 0].grid(True, alpha=0.3)

# 2. 浓度分布演变
time_indices = [0, len(result['C_history'])//3, 2*len(result['C_history'])//3, -1]
for idx, t_idx in enumerate(time_indices):
    axes[0, 1].plot(reactor.x*100, result['C_history'][t_idx], 
                    color=colors[idx], label=f't={t_idx*0.01*100:.0f}s', linewidth=2)

axes[0, 1].set_xlabel('Position (cm)', fontsize=11)
axes[0, 1].set_ylabel('Concentration (mol/m³)', fontsize=11)
axes[0, 1].set_title('Concentration Evolution', fontsize=12, fontweight='bold')
axes[0, 1].legend(fontsize=9)
axes[0, 1].grid(True, alpha=0.3)

# 3. 反应速率分布
axes[1, 0].plot(reactor.x*100, result['r_photo']*1000, 'b-', 
                linewidth=2, label='Photochemical')
axes[1, 0].plot(reactor.x*100, result['r_thermal']*1000, 'r--', 
                linewidth=2, label='Thermal')
axes[1, 0].set_xlabel('Position (cm)', fontsize=11)
axes[1, 0].set_ylabel('Reaction Rate (mmol/m³s)', fontsize=11)
axes[1, 0].set_title('Reaction Rate Distribution', fontsize=12, fontweight='bold')
axes[1, 0].legend(fontsize=9)
axes[1, 0].grid(True, alpha=0.3)

# 4. 转化率随辐射强度变化
conversions = []
for I0 in I0_values:
    res = reactor.simulate(I0=I0, C0=10, T0=298, flow_velocity=0.01,
                          t_max=100, dt=0.01, kappa=5)
    conversions.append(res['C_final'][-1] / 10 * 100)

axes[1, 1].bar(range(len(I0_values)), conversions, color=colors, alpha=0.7)
axes[1, 1].set_xticks(range(len(I0_values)))
axes[1, 1].set_xticklabels([f'{I0}' for I0 in I0_values])
axes[1, 1].set_xlabel('Incident Radiation (W/m²)', fontsize=11)
axes[1, 1].set_ylabel('Remaining Concentration (%)', fontsize=11)
axes[1, 1].set_title('Effect of Radiation Intensity', fontsize=12, fontweight='bold')
axes[1, 1].grid(True, alpha=0.3, axis='y')

plt.tight_layout()
plt.savefig(f'{output_dir}/case6_photocatalytic_reactor.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 案例6图像已保存: case6_photocatalytic_reactor.png")


# ==================== 总结 ====================
print("\n" + "=" * 70)
print("所有仿真案例完成！")
print("=" * 70)
print("\n生成的图像文件：")
print("  1. case1_droplet_evaporation.png - 单液滴蒸发")
print("  2. case2_humid_air_radiation.png - 湿空气辐射传输")
print("  3. case3_porous_drying.png - 多孔介质干燥")
print("  4. case4_spray_cooling.png - 喷雾冷却")
print("  5. case5_solar_still.png - 太阳能蒸馏")
print("  6. case6_photocatalytic_reactor.png - 光催化反应器")
print("\n所有文件保存在：", output_dir)
print("=" * 70)