主题017：多孔介质中的对流换热

多孔介质中的对流换热在能源、环境、化工和生物医学等领域具有广泛的应用。理解多孔介质中的流动与传热机理对于优化相关工程系统至关重要。本教程系统介绍多孔介质的结构特征与输运机理，深入分析达西定律、Brinkman方程和Forchheimer修正等控制方程。详细阐述局部热平衡与非热平衡模型的适用条件，建立多孔介质对流换热的无量纲关联式。通过填充床、热管、地热系统等典型案例，展示多孔介质强化换热的机理与设计方法。教程包含完整的Python仿真代码，实现多孔介质流动与换热的数值模拟，为相关工程应用提供理论基础和计算工具。

摘要

多孔介质中的对流换热在能源、环境、化工和生物医学等领域具有广泛的应用。本教程系统介绍多孔介质的结构特征与输运机理，深入分析达西定律、Brinkman方程和Forchheimer修正等控制方程。详细阐述局部热平衡与非热平衡模型的适用条件，建立多孔介质对流换热的无量纲关联式。通过填充床、热管、地热系统等典型案例，展示多孔介质强化换热的机理与设计方法。教程包含完整的Python仿真代码，实现多孔介质流动与换热的数值模拟，为相关工程应用提供理论基础和计算工具，帮助读者掌握多孔介质传热理论。

关键词

多孔介质，达西定律，Brinkman方程，Forchheimer修正，局部热平衡，非热平衡，对流换热，有效导热系数，渗透率，孔隙率

---

1. 引言

1.1 多孔介质的定义与分类

多孔介质是指由固体骨架和相互连通的孔隙组成的材料。孔隙中充满流体（气体或液体），形成复杂的输运网络。

按孔隙结构分类：

• 颗粒型多孔介质：如砂层、填充床、催化剂颗粒
• 纤维型多孔介质：如滤纸、隔热材料、纺织物
• 泡沫型多孔介质：如金属泡沫、陶瓷泡沫
• 裂隙型多孔介质：如岩石、混凝土、土壤

按孔隙尺度分类：

• 宏观多孔介质：孔隙尺度远大于分子平均自由程
• 介观多孔介质：孔隙尺度与分子平均自由程相当
• 微纳多孔介质：孔隙尺度在微米或纳米量级

按应用领域分类：

• 天然多孔介质：土壤、岩石、生物组织
• 人造多孔介质：催化剂、过滤器、隔热材料
• 功能多孔介质：金属泡沫、相变储能材料

1.2 多孔介质对流换热的重要性

多孔介质对流换热在众多工程领域具有关键作用：

1. 能源工程：地热开发、油藏开采、燃料电池
2. 化工过程：催化反应器、吸附分离、干燥过程
3. 环境科学：地下水污染、土壤修复、大气扩散
4. 生物医学：组织工程、药物输送、血液灌注
5. 电子冷却：热管、均热板、相变储能
6. 建筑节能：多孔保温材料的热性能优化

1.3 本教程内容安排

本教程将从以下几个方面系统阐述多孔介质对流换热：

• 第2节：多孔介质的基本参数与表征
• 第3节：多孔介质流动控制方程
• 第4节：多孔介质传热模型
• 第5节：强迫对流换热分析
• 第6节：自然对流与混合对流
• 第7节：Python仿真实现
• 第8节：工程应用与总结

---

2. 多孔介质的基本参数与表征

2.1 孔隙率与渗透率

2.1.1 孔隙率（Porosity）

孔隙率是多孔介质最基本的结构参数：

$$\epsilon =\frac{V_{void}}{V_{total}}=\frac{V_{void}}{V_{solid}+V_{void}}$$

其中：

• $\epsilon$：孔隙率（无量纲，0 < ε < 1）
• $V_{void}$：孔隙体积
• $V_{solid}$：固体骨架体积
• $V_{total}$：总体积

典型多孔介质的孔隙率：

材料	孔隙率范围
砂层	0.30-0.50
土壤	0.40-0.60
泡沫金属	0.80-0.95
催化剂颗粒	0.40-0.60
纤维材料	0.80-0.99
岩石	0.01-0.30
活性炭	0.50-0.80
陶瓷泡沫	0.70-0.90
金属海绵	0.60-0.90

孔隙率的工程意义：

• 孔隙率越高，流体通过能力越强，但固体骨架强度越低
• 孔隙率影响有效导热系数和比表面积
• 孔隙率分布的均匀性对流动和传热有重要影响

2.1.2 渗透率（Permeability）

渗透率表征多孔介质允许流体通过的能力：

$$K=\frac{\mu QL}{A\Delta P}$$

其中：

• $K$：渗透率（m²）
• $\mu$：流体动力粘度（Pa·s）
• $Q$：体积流量（m³/s）
• $L$：多孔介质厚度（m）
• $A$：横截面积（m²）
• $\Delta P$：压差（Pa）

典型渗透率范围：

材料	渗透率K (m²)
清洁砂	10⁻¹² - 10⁻¹⁰
砂岩	10⁻¹⁵ - 10⁻¹³
石灰岩	10⁻¹⁶ - 10⁻¹⁴
混凝土	10⁻¹⁸ - 10⁻¹⁶
泡沫金属	10⁻⁸ - 10⁻⁶
活性炭	10⁻¹⁹ - 10⁻¹⁷
玻璃纤维	10⁻¹³ - 10⁻¹¹

渗透率的物理意义：

• 渗透率是介质的固有属性，与流体性质无关
• 渗透率越高，流体通过越容易，压降越小
• 渗透率与孔隙率、孔隙连通性密切相关

达西定律的适用条件：

• 流动为层流（低雷诺数）
• 流体为牛顿流体
• 流体与固体不发生化学反应
• 孔隙结构稳定

2.2 比表面积与曲折因子

2.2.1 比表面积（Specific Surface Area）

比表面积定义为单位体积多孔介质的表面积：

$$a_s=\frac{A_{surface}}{V_{total}}\text{(m²/m³)}$$

对于球形颗粒填充床：
$$a_s=\frac{6(1-\epsilon )}{d_p}$$

其中$d_p$为颗粒直径。

2.2.2 曲折因子（Tortuosity）

曲折因子表征孔隙通道的弯曲程度：

$$\tau ={\left(\frac{L_{actual}}{L_{straight}}\right)}^2$$

其中：

• $L_{actual}$：实际流动路径长度
• $L_{straight}$：直线距离

曲折因子通常大于1，典型值为1.5-3.0。

2.3 孔径分布与连通性

2.3.1 孔径分布

多孔介质的孔径分布影响流动和传热特性：

• 平均孔径：$d_{avg}=\frac{4\epsilon }{a_s}$
• 水力直径：$D_h=\frac{4\epsilon }{a_s(1-\epsilon )}$

2.3.2 连通性

孔隙连通性决定流体的可及性：

• 配位数：每个孔隙连接的喉道数
• 连通概率：孔隙相互连通的概率

孔径分布的测量方法：

• 压汞法：适用于大孔径范围
• 气体吸附法：适用于微孔和介孔
• 图像分析法：基于显微图像分析

连通性的重要性：

• 连通性差的介质有效渗透率远低于理论值
• 孤立孔隙不参与流动和传热
• 优先流动通道导致非均匀的温度分布

---

3. 多孔介质流动控制方程

3.1 达西定律（Darcy’s Law）

3.1.1 基本形式

达西定律描述低速流动时压力梯度与速度的关系：

$$\vec{u}=-\frac{K}{\mu }\nabla P$$

其中：

• $\vec{u}$：达西速度（表观速度，m/s）
• $K$：渗透率（m²）
• $\mu$：流体动力粘度（Pa·s）
• $\nabla P$：压力梯度（Pa/m）

物理意义：达西速度是体积平均速度，实际孔隙中的速度（间隙速度）为：
$$\vec{v}=\frac{\vec{u}}{\epsilon }$$

3.1.2 适用范围

达西定律适用于：

• 层流流动（Re < 1-10）
• 牛顿流体
• 不可压缩流动
• 各向同性多孔介质

雷诺数定义：
$$Re=\frac{\rho u{d}_p}{\mu }$$

达西定律的局限性：

• 不适用于高流速（惯性效应显著）
• 不适用于边界附近（粘性效应重要）
• 不适用于各向异性介质
• 不适用于非牛顿流体

修正方法：

• 高流速：Forchheimer修正
• 边界效应：Brinkman修正
• 各向异性：张量形式的渗透率

3.2 Brinkman方程

当需要考虑边界效应或高孔隙率多孔介质时，使用Brinkman修正：

$$\nabla P=-\frac{\mu }{K}\vec{u}+\tilde{\mu }{\nabla }^2\vec{u}$$

其中：

• $\tilde{\mu }$：有效粘度（通常$\tilde{\mu }\approx \mu$）

Brinkman方程在壁面附近自动过渡到粘性边界层。

3.3 Forchheimer修正

对于高速流动，需要考虑惯性效应：

$$\nabla P=-\frac{\mu }{K}\vec{u}-\frac{\rho F}{\sqrt{K}}|\vec{u}|\vec{u}$$

其中：

• $F$：Forchheimer系数（无量纲，通常0.1-1.0）

Forchheimer系数经验关联：
$$F=\frac{1.75}{\sqrt{150{\epsilon }^3}}$$

3.4 通用动量方程

综合考虑粘性、惯性和边界效应：

$$\nabla P=-\frac{\mu }{K}\vec{u}-\frac{\rho F}{\sqrt{K}}|\vec{u}|\vec{u}+\tilde{\mu }{\nabla }^2\vec{u}$$

方程选择的指导原则：

• 低流速、远离边界：达西定律
• 低流速、边界附近：Brinkman方程
• 高流速、远离边界：Forchheimer修正
• 高流速、边界附近：通用动量方程

无量纲数：

• 达西数：$Da=K/L^2$，表征渗透率与特征尺度的比值
• 雷诺数：$Re=\rho u{d}_p/\mu$，表征惯性力与粘性力的比值
• Forchheimer数：$Fo=FReD{a}^{1/2}$，表征惯性效应的相对重要性

---

4. 多孔介质传热模型

4.1 有效导热系数

4.1.1 定义

有效导热系数表征多孔介质的综合导热能力：

$$k_{eff}=\epsilon k_f+(1-\epsilon )k_s$$

其中：

• $k_f$：流体导热系数
• $k_s$：固体导热系数

4.1.2 常用模型

并联模型（上限）：
$$k_{eff,parallel}=\epsilon k_f+(1-\epsilon )k_s$$

串联模型（下限）：
$$\frac{1}{k_{eff,series}}=\frac{\epsilon }{k_f}+\frac{1-\epsilon }{k_s}$$

几何平均模型：
$$k_{eff}=k_f^{\epsilon }{k}_s^{1-\epsilon }$$

经验关联式（适用于填充床）：
$$\frac{k_{eff}}{k_f}=\frac{k_s}{k_f}\left[1-\frac{1-\epsilon }{(1-\epsilon +\epsilon k_f/k_s{)}^{0.5}}\right]$$

影响有效导热系数的因素：

• 孔隙率：孔隙率增加通常降低有效导热系数
• 固液导热系数比：比值越大，有效导热系数越高
• 孔隙结构：连通性好的固体骨架导热性能更好
• 温度：高温下辐射传热可能变得重要

工程应用中的考虑：

• 对于金属泡沫，有效导热系数接近固体值
• 对于绝热材料，有效导热系数接近气体值
• 相变材料填充的多孔介质需考虑相变潜热

4.2 局部热平衡（LTE）模型

假设流体和固体骨架温度相同：

$$(\rho c_p{)}_{eff}\frac{\partial T}{\partial t}+(\rho c_p{)}_f\vec{u}\cdot \nabla T=\nabla \cdot (k_{eff}\nabla T)+q^{\prime \prime \prime }$$

其中有效热容：
$$(\rho c_p{)}_{eff}=\epsilon (\rho c_p{)}_f+(1-\epsilon )(\rho c_p{)}_s$$

适用条件：

• 流体与固体热容相近
• 流速较低
• 孔隙尺度较小

4.3 局部非热平衡（LTNE）模型

分别求解流体和固体的能量方程：

流体能量方程：
$$\epsilon (\rho c_p{)}_f\frac{\partial T_f}{\partial t}+(\rho c_p{)}_f\vec{u}\cdot \nabla T_f=\epsilon k_f{\nabla }^2{T}_f+h_{sf}{a}_s(T_s-T_f)$$

固体能量方程：
$$(1-\epsilon )(\rho c_p{)}_s\frac{\partial T_s}{\partial t}=(1-\epsilon )k_s{\nabla }^2{T}_s-h_{sf}{a}_s(T_s-T_f)$$

其中：

• $h_{sf}$：流体-固体对流换热系数（W/m²·K）
• $a_s$：比表面积（m²/m³）

体积换热系数：
$$h_v=h_{sf}{a}_s\text{(W/m³⋅K)}$$

LTNE模型的适用条件：

• 流体与固体热容差异大
• 流速较高
• 孔隙尺度较大
• 快速瞬态过程

模型选择的判据：

• 毕奥数 $Bi=h_{sf}{a}_s{L}^2/k_{eff}<0.1$：可使用LTE模型
• 达西数 $Da<10^{-6}$：通常满足LTE条件
• 时间尺度分析：热平衡建立时间远小于过程时间

工程应用建议：

• 大多数填充床反应器可使用LTE模型
• 金属泡沫换热器通常需要LTNE模型
• 快速瞬态过程建议使用LTNE模型

4.4 流体-固体换热系数

4.4.1 填充床关联式

Wakao关联式：
$$N{u}_{sf}=\frac{h_{sf}{d}_p}{k_f}=2+1.1R{e}^{0.6}P{r}^{1/3}$$

适用范围：$15<Re<8500$，$0.6<Pr<400$

Gnielinski修正：
$$N{u}_{sf}=2+\sqrt{N{u}_{lam}^2+N{u}_{turb}^2}$$

其中：
$$N{u}_{lam}=0.664R{e}^{0.5}P{r}^{1/3}$$
$$N{u}_{turb}=\frac{0.037R{e}^{0.8}Pr}{1+2.443R{e}^{-0.1}(P{r}^{2/3}-1)}$$

4.4.2 泡沫金属关联式

对于高孔隙率金属泡沫：
$$N{u}_{sf}=C\cdot R{e}^m\cdot P{r}^{1/3}$$

其中C和m取决于泡沫结构（孔径、孔隙率）。

流体-固体换热系数的影响因素：

• 流速：流速增加，换热系数增大
• 孔径：孔径减小，比表面积增大，换热增强
• 流体物性：普朗特数影响温度边界层发展
• 孔隙结构：连通性好的结构有利于换热

工程应用中的考虑：

• 填充床反应器：使用Wakao关联式
• 金属泡沫换热器：使用实验测定值或专门关联式
• 纤维材料：考虑纤维直径和排列方式

---

5. 强迫对流换热分析

5.1 填充床换热

5.1.1 努塞尔数关联式

对于填充床管道流动：

低雷诺数（Re < 100）：
$$N{u}_D=\frac{hD}{k_f}=C_1+C_{2}R{e}_D^{0.5}P{r}^{1/3}$$

高雷诺数（Re > 100）：
$$N{u}_D=C_{3}R{e}_D^{0.8}P{r}^{1/3}$$

其中经验常数：

• $C_1=0.2-0.5$（取决于颗粒形状）
• $C_2=0.1-0.3$
• $C_3=0.02-0.05$

5.1.2 入口段效应

热入口段长度：
$$\frac{L_{t,fd}}{D}=0.05R{e}_{D}Pr\cdot f(\epsilon )$$

其中$f(\epsilon )$为孔隙率修正函数。

5.2 泡沫金属换热

5.2.1 结构参数

金属泡沫的关键参数：

• 孔隙率：通常0.85-0.95
• 孔密度（PPI）：每英寸孔数，10-100 PPI
• 孔径：$d_p=\frac{25.4}{PPI}$（mm）

5.2.2 换热特性

金属泡沫的换热增强机理：

1. 高比表面积：提供大量换热面积
2. 流体混合：复杂孔隙结构增强湍流
3. 导热增强：金属骨架提供高效导热路径

努塞尔数关联式：
$$Nu=C\cdot R{e}^m\cdot P{r}^{1/3}\cdot (1-\epsilon {)}^n$$

典型值：

• $C=0.1-0.5$
• $m=0.5-0.8$
• $n=0.5-1.0$

5.3 压降特性

5.3.1 Ergun方程

Ergun方程描述填充床的压降：

$$\frac{\Delta P}{L}=\frac{150(1-\epsilon {)}^2}{{\epsilon }^3}\frac{\mu u}{d_p^2}+\frac{1.75(1-\epsilon )}{{\epsilon }^3}\frac{\rho u^2}{d_p}$$

第一项：粘性损失（达西项）
第二项：惯性损失（Forchheimer项）

5.3.2 渗透率与惯性系数

从Ergun方程可得：

$$K=\frac{{\epsilon }^3{d}_p^2}{150(1-\epsilon {)}^2}$$

$$F=\frac{1.75}{\sqrt{150{\epsilon }^3}}$$

---

6. 自然对流与混合对流

6.1 自然对流控制方程

6.1.1 达西-瑞利模型

对于低速自然对流，使用达西定律结合Boussinesq近似：

$$\vec{u}=-\frac{K}{\mu }(\nabla P-{\rho }_{0}g\beta (T-T_0))$$

能量方程：
$$(\rho c_p{)}_{eff}\frac{\partial T}{\partial t}+(\rho c_p{)}_f\vec{u}\cdot \nabla T=k_{eff}{\nabla }^{2}T$$

6.1.2 瑞利数

多孔介质自然对流的瑞利数：

$$R{a}^{\ast }=\frac{g\beta K\Delta TL}{\nu {\alpha }_{eff}}$$

其中：

• ${\alpha }_{eff}=\frac{k_{eff}}{(\rho c_p{)}_f}$：有效热扩散系数

6.2 自然对流换热关联式

6.2.1 垂直平板

层流（Ra < 100）*：
$$N{u}_L=0.5R{a}^{\ast 0.5}$$

湍流（Ra > 100）*：
$$N{u}_L=0.2R{a}^{\ast 0.33}$$

6.2.2 水平层

底部加热：

• $R{a}^{\ast }<4{\pi }^2$：纯导热
• $R{a}^{\ast }>4{\pi }^2$：对流开始

临界瑞利数：$R{a}_c^{\ast }=4{\pi }^2\approx 39.5$

6.3 混合对流

6.3.1 混合对流参数

定义混合对流参数：

$$\frac{G{r}^{\ast }}{R{e}^2}=\frac{g\beta K\Delta T}{\nu u}$$

流动状态判据：

• $\frac{G{r}^{\ast }}{R{e}^2}\ll 1$：强迫对流主导
• $\frac{G{r}^{\ast }}{R{e}^2}\approx 1$：混合对流
• $\frac{G{r}^{\ast }}{R{e}^2}\gg 1$：自然对流主导

6.3.2 换热关联式

混合对流努塞尔数：
$$N{u}^n=N{u}_{forced}^n+N{u}_{natural}^n$$

其中n通常取3或4。

---

7. Python仿真实现

7.1 案例一：达西流动与压降分析

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
import matplotlib
matplotlib.use('Agg')
import warnings
warnings.filterwarnings('ignore')

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

print("=" * 60)
print("主题017：多孔介质中的对流换热")
print("=" * 60)

# ============================================
# 案例一：达西流动与压降分析
# ============================================
print("\n案例一：达西流动与压降分析")
print("-" * 50)

def darcy_flow(K, mu, dP_dx):
    """
    计算达西速度
    """
    u = -K / mu * dP_dx
    return u

def ergun_equation(epsilon, dp, mu, rho, u):
    """
    计算Ergun方程压降
    
    参数:
        epsilon: 孔隙率
        dp: 颗粒直径
        mu: 动力粘度
        rho: 密度
        u: 达西速度
    """
    # 粘性项
    viscous_term = 150 * (1 - epsilon)**2 / epsilon**3 * mu * u / dp**2
    
    # 惯性项
    inertial_term = 1.75 * (1 - epsilon) / epsilon**3 * rho * u**2 / dp
    
    # 总压降梯度
    dP_dx = viscous_term + inertial_term
    
    return dP_dx, viscous_term, inertial_term

def permeability_from_ergun(epsilon, dp):
    """
    从Ergun方程计算渗透率
    """
    K = epsilon**3 * dp**2 / (150 * (1 - epsilon)**2)
    return K

# 参数设置
epsilon_values = np.linspace(0.3, 0.95, 20)
dp = 0.001  # 颗粒直径 1mm
mu = 1e-3   # 水的粘度
rho = 1000  # 水的密度
u = 0.01    # 达西速度 0.01 m/s

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

# 渗透率随孔隙率变化
ax1 = axes[0, 0]
K_values = [permeability_from_ergun(eps, dp) for eps in epsilon_values]
ax1.semilogy(epsilon_values, K_values, 'b-', linewidth=2)
ax1.set_xlabel('Porosity', fontsize=11)
ax1.set_ylabel('Permeability K (m²)', fontsize=11)
ax1.set_title('Permeability vs Porosity', fontsize=12)
ax1.grid(True, alpha=0.3, which='both')

# Ergun方程压降分解
ax2 = axes[0, 1]
viscous_contrib = []
inertial_contrib = []
total_contrib = []

for eps in epsilon_values:
    dP, visc, inert = ergun_equation(eps, dp, mu, rho, u)
    viscous_contrib.append(visc)
    inertial_contrib.append(inert)
    total_contrib.append(dP)

ax2.semilogy(epsilon_values, viscous_contrib, 'b-', linewidth=2, label='Viscous term')
ax2.semilogy(epsilon_values, inertial_contrib, 'r--', linewidth=2, label='Inertial term')
ax2.semilogy(epsilon_values, total_contrib, 'k-', linewidth=2, label='Total')
ax2.set_xlabel('Porosity', fontsize=11)
ax2.set_ylabel('Pressure Gradient (Pa/m)', fontsize=11)
ax2.set_title('Ergun Equation Terms', fontsize=12)
ax2.legend(fontsize=9)
ax2.grid(True, alpha=0.3, which='both')

# 不同流速下的压降
ax3 = axes[1, 0]
u_range = np.linspace(0.001, 0.1, 50)
epsilon_cases = [0.4, 0.6, 0.8]
colors = ['blue', 'green', 'red']

for i, eps in enumerate(epsilon_cases):
    dP_values = []
    for u_val in u_range:
        dP, _, _ = ergun_equation(eps, dp, mu, rho, u_val)
        dP_values.append(dP)
    ax3.plot(u_range, dP_values, color=colors[i], linewidth=2, label=f'ε = {eps}')

ax3.set_xlabel('Velocity (m/s)', fontsize=11)
ax3.set_ylabel('Pressure Gradient (Pa/m)', fontsize=11)
ax3.set_title('Pressure Drop vs Velocity', fontsize=12)
ax3.legend(fontsize=9)
ax3.grid(True, alpha=0.3)

# 达西定律验证
ax4 = axes[1, 1]
K_test = 1e-10  # 测试渗透率
mu_test = 1e-3
dP_range = np.linspace(-1000, 1000, 100)
u_darcy = [-K_test / mu_test * dP for dP in dP_range]

ax4.plot(dP_range, u_darcy, 'b-', linewidth=2)
ax4.axhline(y=0, color='k', linestyle='-', alpha=0.3)
ax4.axvline(x=0, color='k', linestyle='-', alpha=0.3)
ax4.set_xlabel('Pressure Gradient (Pa/m)', fontsize=11)
ax4.set_ylabel('Darcy Velocity (m/s)', fontsize=11)
ax4.set_title("Darcy's Law", fontsize=12)
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('case1_darcy_flow.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 达西流动与压降分析图已保存")

7.2 案例二：有效导热系数计算

# ============================================
# 案例二：有效导热系数计算
# ============================================
print("\n案例二：有效导热系数计算")
print("-" * 50)

def effective_conductivity_parallel(epsilon, kf, ks):
    """并联模型（上限）"""
    return epsilon * kf + (1 - epsilon) * ks

def effective_conductivity_series(epsilon, kf, ks):
    """串联模型（下限）"""
    return 1 / (epsilon / kf + (1 - epsilon) / ks)

def effective_conductivity_geometric(epsilon, kf, ks):
    """几何平均模型"""
    return kf**epsilon * ks**(1-epsilon)

def effective_conductivity_kunii(epsilon, kf, ks):
    """Kunii-Smith模型（填充床）"""
    k_ratio = ks / kf
    k_eff = kf * (k_ratio * (1 - (1 - epsilon) / (1 - epsilon + epsilon / k_ratio)**0.5))
    return k_eff

# 参数设置
kf = 0.6  # 水的导热系数 W/m·K
ks_values = [1, 10, 50, 100, 400]  # 不同固体导热系数
materials = ['Soil', 'Sand', 'Ceramic', 'Steel', 'Copper']
epsilon_range = np.linspace(0.2, 0.9, 50)

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

# 不同模型的对比
ax1 = axes[0, 0]
ks = 50  # 陶瓷
epsilon_plot = np.linspace(0.2, 0.9, 100)

k_parallel = [effective_conductivity_parallel(eps, kf, ks) for eps in epsilon_plot]
k_series = [effective_conductivity_series(eps, kf, ks) for eps in epsilon_plot]
k_geometric = [effective_conductivity_geometric(eps, kf, ks) for eps in epsilon_plot]

ax1.plot(epsilon_plot, k_parallel, 'b-', linewidth=2, label='Parallel (Upper bound)')
ax1.plot(epsilon_plot, k_series, 'r--', linewidth=2, label='Series (Lower bound)')
ax1.plot(epsilon_plot, k_geometric, 'g:', linewidth=2, label='Geometric mean')
ax1.set_xlabel('Porosity', fontsize=11)
ax1.set_ylabel('k_eff (W/m·K)', fontsize=11)
ax1.set_title(f'Effective Conductivity Models (k_s = {ks})', fontsize=12)
ax1.legend(fontsize=9)
ax1.grid(True, alpha=0.3)

# 不同固体材料的影响
ax2 = axes[0, 1]
epsilon_fixed = 0.4
k_eff_materials = []

for i, (ks, mat) in enumerate(zip(ks_values, materials)):
    k_eff = effective_conductivity_parallel(epsilon_fixed, kf, ks)
    k_eff_materials.append(k_eff)
    print(f"  {mat:8s} (k_s = {ks:3d}): k_eff = {k_eff:.2f} W/m·K")

colors = plt.cm.viridis(np.linspace(0, 1, len(materials)))
bars = ax2.bar(range(len(materials)), k_eff_materials, color=colors, alpha=0.7, edgecolor='black')
ax2.set_xticks(range(len(materials)))
ax2.set_xticklabels(materials, fontsize=9, rotation=45)
ax2.set_ylabel('k_eff (W/m·K)', fontsize=11)
ax2.set_title(f'Effective Conductivity (ε = {epsilon_fixed})', fontsize=12)
ax2.grid(True, alpha=0.3, axis='y')

# 孔隙率对k_eff的影响（不同材料）
ax3 = axes[1, 0]
for i, (ks, mat) in enumerate(zip(ks_values, materials)):
    k_eff_range = [effective_conductivity_parallel(eps, kf, ks) for eps in epsilon_range]
    ax3.plot(epsilon_range, k_eff_range, color=colors[i], linewidth=2, label=mat)

ax3.set_xlabel('Porosity', fontsize=11)
ax3.set_ylabel('k_eff (W/m·K)', fontsize=11)
ax3.set_title('k_eff vs Porosity for Different Materials', fontsize=12)
ax3.legend(fontsize=8)
ax3.grid(True, alpha=0.3)

# 导热系数比的影响
ax4 = axes[1, 1]
k_ratio_range = np.logspace(0, 3, 50)
epsilon_cases = [0.3, 0.5, 0.7, 0.9]

for i, eps in enumerate(epsilon_cases):
    k_eff_ratio = []
    for k_ratio in k_ratio_range:
        k_eff = effective_conductivity_parallel(eps, 1.0, k_ratio)
        k_eff_ratio.append(k_eff)
    ax4.semilogx(k_ratio_range, k_eff_ratio, linewidth=2, label=f'ε = {eps}')

ax4.set_xlabel('k_s/k_f Ratio', fontsize=11)
ax4.set_ylabel('k_eff/k_f', fontsize=11)
ax4.set_title('Normalized Effective Conductivity', fontsize=12)
ax4.legend(fontsize=9)
ax4.grid(True, alpha=0.3, which='both')

plt.tight_layout()
plt.savefig('case2_effective_conductivity.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 有效导热系数计算图已保存")

7.3 案例三：填充床对流换热

# ============================================
# 案例三：填充床对流换热
# ============================================
print("\n案例三：填充床对流换热")
print("-" * 50)

def nusselt_packed_bed(Re, Pr, epsilon):
    """
    填充床努塞尔数关联式（Wakao）
    """
    Nu = 2 + 1.1 * Re**0.6 * Pr**(1/3)
    return Nu

def h_sf_from_nusselt(Nu, kf, dp):
    """
    从努塞尔数计算对流换热系数
    """
    h_sf = Nu * kf / dp
    return h_sf

def specific_surface_area(epsilon, dp):
    """
    计算比表面积（球形颗粒）
    """
    a_s = 6 * (1 - epsilon) / dp
    return a_s

def volume_heat_transfer_coeff(h_sf, a_s):
    """
    计算体积换热系数
    """
    h_v = h_sf * a_s
    return h_v

# 参数设置
dp = 0.005  # 颗粒直径 5mm
kf = 0.6    # 流体导热系数
Pr = 7.0    # 普朗特数（水）
epsilon_values = [0.35, 0.4, 0.45, 0.5]
Re_range = np.logspace(0, 3, 50)

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

# 努塞尔数随雷诺数变化
ax1 = axes[0, 0]
colors = plt.cm.viridis(np.linspace(0, 1, len(epsilon_values)))

for i, eps in enumerate(epsilon_values):
    Nu_values = [nusselt_packed_bed(Re, Pr, eps) for Re in Re_range]
    ax1.loglog(Re_range, Nu_values, color=colors[i], linewidth=2, label=f'ε = {eps}')

ax1.axhline(y=2, color='r', linestyle='--', alpha=0.5, label='Conduction limit')
ax1.set_xlabel('Reynolds Number', fontsize=11)
ax1.set_ylabel('Nu_sf', fontsize=11)
ax1.set_title('Nusselt Number vs Re (Wakao Correlation)', fontsize=12)
ax1.legend(fontsize=9)
ax1.grid(True, alpha=0.3, which='both')

# 对流换热系数
ax2 = axes[0, 1]
Re_fixed = 100
h_sf_values = []
a_s_values = []

for eps in epsilon_values:
    Nu = nusselt_packed_bed(Re_fixed, Pr, eps)
    h_sf = h_sf_from_nusselt(Nu, kf, dp)
    a_s = specific_surface_area(eps, dp)
    h_sf_values.append(h_sf)
    a_s_values.append(a_s)
    print(f"  ε = {eps:.2f}: h_sf = {h_sf:.1f} W/m²·K, a_s = {a_s:.0f} m²/m³")

bars = ax2.bar(range(len(epsilon_values)), h_sf_values, color=colors, alpha=0.7, edgecolor='black')
ax2.set_xticks(range(len(epsilon_values)))
ax2.set_xticklabels([f'ε={eps}' for eps in epsilon_values], fontsize=9)
ax2.set_ylabel('h_sf (W/m²·K)', fontsize=11)
ax2.set_title(f'Convection Coefficient (Re = {Re_fixed})', fontsize=12)
ax2.grid(True, alpha=0.3, axis='y')

# 体积换热系数
ax3 = axes[1, 0]
h_v_values = []
for i, eps in enumerate(epsilon_values):
    h_v = volume_heat_transfer_coeff(h_sf_values[i], a_s_values[i])
    h_v_values.append(h_v)

bars = ax3.bar(range(len(epsilon_values)), h_v_values, color=colors, alpha=0.7, edgecolor='black')
ax3.set_xticks(range(len(epsilon_values)))
ax3.set_xticklabels([f'ε={eps}' for eps in epsilon_values], fontsize=9)
ax3.set_ylabel('h_v (W/m³·K)', fontsize=11)
ax3.set_title('Volume Heat Transfer Coefficient', fontsize=12)
ax3.grid(True, alpha=0.3, axis='y')

# 比表面积随孔隙率变化
ax4 = axes[1, 1]
epsilon_range = np.linspace(0.3, 0.6, 50)
a_s_range = [specific_surface_area(eps, dp) for eps in epsilon_range]

ax4.plot(epsilon_range, a_s_range, 'b-', linewidth=2)
ax4.set_xlabel('Porosity', fontsize=11)
ax4.set_ylabel('Specific Surface Area (m²/m³)', fontsize=11)
ax4.set_title(f'Specific Surface Area (d_p = {dp*1000:.0f} mm)', fontsize=12)
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('case3_packed_bed.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 填充床对流换热图已保存")

7.4 案例四：金属泡沫换热

# ============================================
# 案例四：金属泡沫换热
# ============================================
print("\n案例四：金属泡沫换热")
print("-" * 50)

def metal_foam_nusselt(Re, Pr, epsilon, PPI):
    """
    金属泡沫努塞尔数关联式
    
    参数:
        Re: 雷诺数
        Pr: 普朗特数
        epsilon: 孔隙率
        PPI: 每英寸孔数
    """
    # 经验关联式
    C = 0.2 + 0.1 * (1 - epsilon)
    m = 0.6 + 0.1 * epsilon
    n_exp = 0.8 - 0.2 * epsilon
    
    Nu = C * Re**m * Pr**(1/3) * (1 - epsilon)**n_exp
    return Nu

def metal_foam_permeability(epsilon, PPI):
    """
    金属泡沫渗透率估算
    """
    # 孔径
    d_pore = 25.4e-3 / PPI  # m
    
    # 经验公式
    K = epsilon**3 * d_pore**2 / (150 * (1 - epsilon)**2)
    return K

def metal_foam_surface_area(epsilon, PPI):
    """
    金属泡沫比表面积估算
    """
    d_pore = 25.4e-3 / PPI
    a_s = 6 * (1 - epsilon) / d_pore
    return a_s

# 参数设置
PPI_values = [10, 20, 40, 60, 80]
epsilon_foam = 0.9
Pr = 0.7  # 空气
Re_range = np.logspace(1, 4, 50)

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

# 不同PPI的努塞尔数
ax1 = axes[0, 0]
colors = plt.cm.plasma(np.linspace(0, 1, len(PPI_values)))

for i, PPI in enumerate(PPI_values):
    Nu_values = [metal_foam_nusselt(Re, Pr, epsilon_foam, PPI) for Re in Re_range]
    ax1.loglog(Re_range, Nu_values, color=colors[i], linewidth=2, label=f'{PPI} PPI')

ax1.set_xlabel('Reynolds Number', fontsize=11)
ax1.set_ylabel('Nu', fontsize=11)
ax1.set_title(f'Metal Foam Heat Transfer (ε = {epsilon_foam})', fontsize=12)
ax1.legend(fontsize=9)
ax1.grid(True, alpha=0.3, which='both')

# 渗透率随PPI变化
ax2 = axes[0, 1]
K_foam_values = [metal_foam_permeability(epsilon_foam, PPI) for PPI in PPI_values]

bars = ax2.bar(range(len(PPI_values)), K_foam_values, color=colors, alpha=0.7, edgecolor='black')
ax2.set_xticks(range(len(PPI_values)))
ax2.set_xticklabels([f'{PPI}' for PPI in PPI_values], fontsize=9)
ax2.set_xlabel('PPI', fontsize=11)
ax2.set_ylabel('Permeability K (m²)', fontsize=11)
ax2.set_title(f'Permeability vs PPI (ε = {epsilon_foam})', fontsize=12)
ax2.grid(True, alpha=0.3, axis='y')

# 比表面积随PPI变化
ax3 = axes[1, 0]
a_s_foam_values = [metal_foam_surface_area(epsilon_foam, PPI) for PPI in PPI_values]

bars = ax3.bar(range(len(PPI_values)), a_s_foam_values, color=colors, alpha=0.7, edgecolor='black')
ax3.set_xticks(range(len(PPI_values)))
ax3.set_xticklabels([f'{PPI}' for PPI in PPI_values], fontsize=9)
ax3.set_xlabel('PPI', fontsize=11)
ax3.set_ylabel('Specific Surface Area (m²/m³)', fontsize=11)
ax3.set_title(f'Surface Area vs PPI (ε = {epsilon_foam})', fontsize=12)
ax3.grid(True, alpha=0.3, axis='y')

# 孔隙率对换热的影响
ax4 = axes[1, 1]
epsilon_range = np.linspace(0.85, 0.98, 20)
PPI_fixed = 40
Re_fixed = 1000

Nu_epsilon = [metal_foam_nusselt(Re_fixed, Pr, eps, PPI_fixed) for eps in epsilon_range]

ax4.plot(epsilon_range, Nu_epsilon, 'b-', linewidth=2)
ax4.set_xlabel('Porosity', fontsize=11)
ax4.set_ylabel('Nu', fontsize=11)
ax4.set_title(f'Nu vs Porosity ({PPI_fixed} PPI, Re = {Re_fixed})', fontsize=12)
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('case4_metal_foam.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 金属泡沫换热图已保存")

7.5 案例五：自然对流与瑞利数

# ============================================
# 案例五：自然对流与瑞利数
# ============================================
print("\n案例五：自然对流与瑞利数")
print("-" * 50)

def rayleigh_number_porous(g, beta, K, delta_T, L, nu, alpha_eff):
    """
    计算多孔介质瑞利数
    """
    Ra_star = g * beta * K * delta_T * L / (nu * alpha_eff)
    return Ra_star

def nusselt_natural_convection(Ra_star, geometry='vertical'):
    """
    自然对流努塞尔数关联式
    """
    if geometry == 'vertical':
        if Ra_star < 100:  # 层流
            Nu = 0.5 * Ra_star**0.5
        else:  # 湍流
            Nu = 0.2 * Ra_star**0.33
    elif geometry == 'horizontal':
        if Ra_star < 40:  # 纯导热
            Nu = 1.0
        else:
            Nu = 0.2 * Ra_star**0.33
    return Nu

# 参数设置
g = 9.81
beta = 2e-4  # 热膨胀系数
K = 1e-8     # 渗透率
delta_T = 50  # 温差
L = 0.1      # 特征长度
nu = 1e-6    # 运动粘度
alpha_eff = 1e-7  # 有效热扩散系数

Ra_range = np.logspace(0, 4, 50)

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

# 努塞尔数随瑞利数变化
ax1 = axes[0, 0]
Nu_vertical = [nusselt_natural_convection(Ra, 'vertical') for Ra in Ra_range]
Nu_horizontal = [nusselt_natural_convection(Ra, 'horizontal') for Ra in Ra_range]

ax1.loglog(Ra_range, Nu_vertical, 'b-', linewidth=2, label='Vertical plate')
ax1.loglog(Ra_range, Nu_horizontal, 'r--', linewidth=2, label='Horizontal layer')
ax1.axvline(x=40, color='gray', linestyle=':', alpha=0.5, label='Critical Ra*')
ax1.set_xlabel('Ra*', fontsize=11)
ax1.set_ylabel('Nu', fontsize=11)
ax1.set_title('Natural Convection Heat Transfer', fontsize=12)
ax1.legend(fontsize=9)
ax1.grid(True, alpha=0.3, which='both')

# 渗透率对自然对流的影响
ax2 = axes[0, 1]
K_range = np.logspace(-10, -6, 50)
Ra_K = [rayleigh_number_porous(g, beta, K_val, delta_T, L, nu, alpha_eff) for K_val in K_range]

ax2.semilogx(K_range, Ra_K, 'b-', linewidth=2)
ax2.axhline(y=40, color='r', linestyle='--', alpha=0.5, label='Critical Ra* = 40')
ax2.set_xlabel('Permeability K (m²)', fontsize=11)
ax2.set_ylabel('Ra*', fontsize=11)
ax2.set_title('Ra* vs Permeability', fontsize=12)
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3, which='both')

# 温差对自然对流的影响
ax3 = axes[1, 0]
delta_T_range = np.linspace(10, 100, 50)
Ra_dT = [rayleigh_number_porous(g, beta, K, dT, L, nu, alpha_eff) for dT in delta_T_range]
Nu_dT = [nusselt_natural_convection(Ra, 'vertical') for Ra in Ra_dT]

ax2_twin = ax3.twinx()
ax3.plot(delta_T_range, Ra_dT, 'b-', linewidth=2, label='Ra*')
ax2_twin.plot(delta_T_range, Nu_dT, 'r--', linewidth=2, label='Nu')
ax3.set_xlabel('Temperature Difference (K)', fontsize=11)
ax3.set_ylabel('Ra*', fontsize=11, color='b')
ax2_twin.set_ylabel('Nu', fontsize=11, color='r')
ax3.set_title('Natural Convection vs Temperature Difference', fontsize=12)
ax3.grid(True, alpha=0.3)

# 流动状态图
ax4 = axes[1, 1]
ax4.axis('off')

flow_regime_text = """
Natural Convection in Porous Media:

Flow Regimes:

1. Conduction Dominated (Ra* < 40)
   • Heat transfer by conduction
   • Nu ≈ 1
   • No fluid motion

2. Convection Dominated (Ra* > 40)
   • Cellular flow patterns
   • Enhanced heat transfer
   • Nu ∝ Ra*^0.5 (laminar)
   • Nu ∝ Ra*^0.33 (turbulent)

Critical Rayleigh Number:
• Vertical plate: Ra*_c ≈ 40
• Horizontal layer: Ra*_c = 4π² ≈ 39.5
• Confined cavity: Depends on aspect ratio

Engineering Applications:
• Geothermal systems
• Building insulation
• Electronic cooling
• Thermal energy storage
"""

ax4.text(0.05, 0.95, flow_regime_text, fontsize=9, verticalalignment='top',
        family='monospace', bbox=dict(boxstyle='round', facecolor='lightyellow', alpha=0.5))

plt.tight_layout()
plt.savefig('case5_natural_convection.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 自然对流与瑞利数图已保存")

7.6 综合对比图

# ============================================
# 综合对比图
# ============================================
print("\n生成综合对比图...")

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

# 1. 不同多孔介质的渗透率对比
ax1 = axes[0, 0]
media_types = ['Sand', 'Sandstone', 'Limestone', 'Concrete', 'Metal Foam']
K_typical = [1e-11, 1e-14, 1e-15, 1e-17, 1e-7]
colors_media = plt.cm.viridis(np.linspace(0, 1, len(media_types)))

bars = ax1.barh(media_types, K_typical, color=colors_media, alpha=0.7, edgecolor='black')
ax1.set_xlabel('Permeability K (m²)', fontsize=10)
ax1.set_title('Permeability of Different Media', fontsize=11)
ax1.set_xscale('log')
ax1.grid(True, alpha=0.3, axis='x')

# 2. 换热增强效果对比
ax2 = axes[0, 1]
enhancement_types = ['Packed Bed', 'Metal Foam', 'Micro-channel', 'Finned Tube']
enhancement_factor = [2.5, 8.0, 15.0, 5.0]
colors_enhance = ['blue', 'red', 'green', 'orange']

bars = ax2.bar(enhancement_types, enhancement_factor, color=colors_enhance, alpha=0.7, edgecolor='black')
ax2.set_ylabel('Enhancement Factor', fontsize=10)
ax2.set_title('Heat Transfer Enhancement', fontsize=11)
ax2.set_xticklabels(enhancement_types, fontsize=8, rotation=45)
ax2.grid(True, alpha=0.3, axis='y')

# 3. 不同模型的适用范围
ax3 = axes[0, 2]
models = ["Darcy's\nLaw", 'Brinkman\nModel', 'Forchheimer\nCorrection', 'Full\nEquation']
Re_range_model = ['Re < 1', 'Re < 10', 'Re > 1', 'All Re']
colors_model = ['green', 'yellow', 'orange', 'blue']

for i, (model, Re, color) in enumerate(zip(models, Re_range_model, colors_model)):
    ax3.barh(i, 1, color=color, alpha=0.7, edgecolor='black')
    ax3.text(0.5, i, f'{model}\n({Re})', ha='center', va='center', fontsize=8)

ax3.set_xlim(0, 1)
ax3.set_ylim(-0.5, len(models) - 0.5)
ax3.set_yticks([])
ax3.set_xticks([])
ax3.set_title('Model Applicability', fontsize=11)
ax3.axis('off')

# 4. 孔隙率对各项性能的影响
ax4 = axes[1, 0]
epsilon_analysis = np.linspace(0.3, 0.95, 50)

# 归一化性能指标
k_eff_norm = [(eps * 0.6 + (1-eps) * 50) / 50 for eps in epsilon_analysis]
K_norm = [eps**3 / (1-eps)**2 / 10 for eps in epsilon_analysis]
a_s_norm = [(1-eps) * 6 / 0.001 / 3000 for eps in epsilon_analysis]

ax4.plot(epsilon_analysis, k_eff_norm, 'b-', linewidth=2, label='k_eff (normalized)')
ax4.plot(epsilon_analysis, K_norm, 'r--', linewidth=2, label='K (normalized)')
ax4.plot(epsilon_analysis, a_s_norm, 'g:', linewidth=2, label='a_s (normalized)')
ax4.set_xlabel('Porosity', fontsize=10)
ax4.set_ylabel('Normalized Value', fontsize=10)
ax4.set_title('Porosity Effect on Properties', fontsize=11)
ax4.legend(fontsize=8)
ax4.grid(True, alpha=0.3)

# 5. 工程应用分布
ax5 = axes[1, 1]
applications = ['Energy', 'Chemical', 'Environmental', 'Biomedical', 'Electronics']
usage_frequency = [30, 25, 20, 15, 10]
colors_app = plt.cm.plasma(np.linspace(0, 1, len(applications)))

wedges, texts, autotexts = ax5.pie(usage_frequency, labels=applications, autopct='%1.0f%%',
                                     colors=colors_app, startangle=90)
ax5.set_title('Application Distribution', fontsize=11)

# 6. 关键公式总结
ax6 = axes[1, 2]
ax6.axis('off')

formula_text = """
Key Equations:

Darcy's Law:
u = -(K/μ)∇P

Ergun Equation:
ΔP/L = 150(1-ε)²μu/(ε³dp²)
       + 1.75(1-ε)ρu²/(ε³dp)

Effective Conductivity:
k_eff = εkf + (1-ε)ks

Nusselt Number (Packed Bed):
Nu_sf = 2 + 1.1Re^0.6Pr^(1/3)

Volume Heat Transfer:
h_v = h_sf × a_s

Rayleigh Number:
Ra* = gβKΔTL/(να_eff)

Natural Convection:
Nu = 0.5Ra*^0.5 (laminar)
Nu = 0.2Ra*^0.33 (turbulent)
"""

ax6.text(0.05, 0.95, formula_text, fontsize=8.5, verticalalignment='top',
        family='monospace', bbox=dict(boxstyle='round', facecolor='lightyellow', alpha=0.5))

plt.tight_layout()
plt.savefig('comprehensive_comparison.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 综合对比图已保存")

print("\n" + "=" * 60)
print("所有仿真案例已完成！")
print("=" * 60)
print("\n生成的文件：")
print("  1. case1_darcy_flow.png - 达西流动与压降分析")
print("  2. case2_effective_conductivity.png - 有效导热系数计算")
print("  3. case3_packed_bed.png - 填充床对流换热")
print("  4. case4_metal_foam.png - 金属泡沫换热")
print("  5. case5_natural_convection.png - 自然对流与瑞利数")
print("  6. comprehensive_comparison.png - 综合对比")
print("\n输出目录：d:\\文档\\500仿真领域\\工程仿真\\对流换热仿真\\主题017_多孔介质中的对流换热")
print("=" * 60)

---