传热学仿真-主题088-多孔介质相变传热

kkchenjj

157人浏览 · 2026-03-21 04:55:07

kkchenjj · 2026-03-21 04:55:07 发布

第八十八篇：多孔介质相变传热

摘要

多孔介质相变传热是地热工程、蓄能技术和石油开采等领域的重要研究课题，涉及复杂的传热传质耦合过程。本教程系统介绍多孔介质中相变的物理机制，讨论局部热平衡与非平衡模型、有效介质理论以及数值求解方法。通过Python实现多孔介质中相变的数值模拟，包括温度场分布、相变界面移动和传热性能分析。本教程涵盖达西定律、渗透率、孔隙率等核心内容。通过本教程的学习，读者将掌握多孔介质相变传热的建模方法，能够进行地热系统和蓄能装置的优化设计。

关键词

多孔介质，相变传热，局部热平衡，达西定律，有效介质理论，渗透率，地热系统

在这里插入图片描述

1. 多孔介质基础

1.1 基本参数

孔隙率：
$ϕ=VvoidVtotal\phi = \frac{V_{void}}{V_{total}}$

渗透率（达西定律）：
$u=−Kμ∇p\mathbf{u} = -\frac{K}{\mu} \nabla p$

1.2 有效介质理论

有效导热系数：
$keff=ϕkf+(1−ϕ)ksk_{eff} = \phi k_f + (1-\phi) k_s$

有效热容：
$(ρcp)eff=ϕ(ρcp)f+(1−ϕ)(ρcp)s(\rho c_p)_{eff} = \phi (\rho c_p)_f + (1-\phi) (\rho c_p)_s$

1.3 局部热平衡

局部热平衡假设：
$T_s = T_f = T$

能量方程：
$(ρcp)eff∂T∂t=∇⋅(keff∇T)+ϕρL∂fl∂t(\rho c_p)_{eff} \frac{\partial T}{\partial t} = \nabla \cdot (k_{eff} \nabla T) + \phi \rho L \frac{\partial f_l}{\partial t}$

2. 数值模型

2.1 相变模型

液相分数：
$fl={0T<TsT−TsTl−TsTs≤T≤Tl1T>Tlf_l = \begin{cases} 0 & T < T_s \\ \frac{T-T_s}{T_l-T_s} & T_s \leq T \leq T_l \\ 1 & T > T_l \end{cases}$

2.2 有效热容法

等效热容：
$cp,eff={csT<Tscm+LTl−TsTs≤T≤TlclT>Tlc_{p,eff} = \begin{cases} c_s & T < T_s \\ c_m + \frac{L}{T_l-T_s} & T_s \leq T \leq T_l \\ c_l & T > T_l \end{cases}$

2.3 边界条件

恒定温度边界：
$T_w \text{ at } r = r_w$

绝热边界：
$r=re\frac{\partial T}{\partial n} = 0 \text{ at } r = r_e$

3. Python仿真实现

3.1 多孔介质相变模拟

"""
主题088：多孔介质相变传热
多孔介质中相变过程仿真
"""
import numpy as np
import matplotlib.pyplot as plt
import os

# 创建输出目录
output_dir = r'd:\文档\500仿真领域\工程仿真\传热学仿真\主题088'
os.makedirs(output_dir, exist_ok=True)

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

# 多孔介质参数（岩石-水系统）
phi = 0.3  # 孔隙率
k_rock = 2.5  # W/(m·K)
k_water = 0.6  # W/(m·K)
rho_rock = 2500  # kg/m³
rho_water = 1000  # kg/m³
c_rock = 800  # J/(kg·K)
c_water = 4186  # J/(kg·K)

# 有效参数
k_eff = phi * k_water + (1-phi) * k_rock
rho_eff = phi * rho_water + (1-phi) * rho_rock
c_eff = (phi * rho_water * c_water + (1-phi) * rho_rock * c_rock) / rho_eff
alpha_eff = k_eff / (rho_eff * c_eff)

# 相变参数
T_melt = 273.15  # K
L_fusion = 334e3  # J/kg
delta_T = 1.0  # K

def effective_properties(T):
    """计算等效物性参数"""
    if T < (T_melt - delta_T):
        c_p = c_eff
        k = k_eff
        rho = rho_eff
    elif T > (T_melt + delta_T):
        c_p = c_eff
        k = k_eff
        rho = rho_eff
    else:
        # 相变区域
        c_p = c_eff + L_fusion / (2 * delta_T)
        k = k_eff
        rho = rho_eff
    
    return c_p, k, rho

def solve_porous_media_phase_change(Nr, Nt, dt, t_end):
    """多孔介质相变数值求解"""
    r_w = 0.1  # m (井筒半径)
    r_e = 10.0  # m (外边界半径)
    
    r = np.linspace(r_w, r_e, Nr)
    dr = r[1] - r[0]
    
    # 初始温度场
    T = np.ones(Nr) * 278.15  # K (5°C)
    T[0] = 263.15  # K (-10°C，井筒温度)
    
    T_history = [T.copy()]
    t_history = [0]
    
    t = 0
    while t < t_end:
        T_new = T.copy()
        
        # 内部节点（柱坐标）
        for i in range(1, Nr-1):
            c_p, k, rho = effective_properties(T[i])
            alpha = k / (rho * c_p)
            
            # 柱坐标拉普拉斯算子
            dT_dr = (T[i+1] - T[i-1]) / (2*dr)
            d2T_dr2 = (T[i+1] - 2*T[i] + T[i-1]) / dr**2
            
            laplacian = d2T_dr2 + (1/r[i]) * dT_dr
            
            T_new[i] = T[i] + alpha * dt * laplacian
        
        # 边界条件
        T_new[0] = 263.15  # 恒定温度
        T_new[-1] = T[-2]  # 绝热边界
        
        T = T_new
        t += dt
        
        if len(T_history) < 50 and t % (t_end/20) < dt:
            T_history.append(T.copy())
            t_history.append(t)
    
    return r, T_history, t_history

# 1. 多孔介质相变模拟
fig, axes = plt.subplots(2, 2, figsize=(14, 12))

Nr, Nt = 100, 1000
dt = 3600  # 1小时
t_end = 86400 * 30  # 30天

r, T_history, t_history = solve_porous_media_phase_change(Nr, Nt, dt, t_end)

# 温度分布
for i, (T, t) in enumerate(zip(T_history, t_history)):
    if i % 5 == 0:
        axes[0, 0].plot(r, T - 273.15, linewidth=2, label=f't={t/86400:.0f}天' if i < 15 else '')

axes[0, 0].axhline(y=0, color='r', linestyle='--', label='冰点')
axes[0, 0].set_xlabel('半径 (m)', fontsize=11)
axes[0, 0].set_ylabel('温度 (°C)', fontsize=11)
axes[0, 0].set_title('多孔介质温度分布', fontsize=12)
axes[0, 0].legend()
axes[0, 0].grid(True, alpha=0.3)

# 冻结前沿位置
frozen_front = []
for T in T_history:
    frozen_idx = np.where(T <= T_melt)[0]
    if len(frozen_idx) > 0:
        frozen_front.append(r[frozen_idx[-1]])
    else:
        frozen_front.append(r[0])

axes[0, 1].plot(np.array(t_history)/86400, np.array(frozen_front), 'b-o', linewidth=2, markersize=4)
axes[0, 1].set_xlabel('时间 (天)', fontsize=11)
axes[0, 1].set_ylabel('冻结半径 (m)', fontsize=11)
axes[0, 1].set_title('冻结前沿演化', fontsize=12)
axes[0, 1].grid(True, alpha=0.3)

# 不同孔隙率对比
porosities = [0.2, 0.3, 0.4]
colors = ['blue', 'green', 'red']

for phi_val, color in zip(porosities, colors):
    k_eff_val = phi_val * k_water + (1-phi_val) * k_rock
    rho_eff_val = phi_val * rho_water + (1-phi_val) * rho_rock
    c_eff_val = (phi_val * rho_water * c_water + (1-phi_val) * rho_rock * c_rock) / rho_eff_val
    alpha_val = k_eff_val / (rho_eff_val * c_eff_val)
    
    # 简化的温度衰减
    t_sim = np.linspace(0, t_end, 100)
    T_center = 278.15 - (278.15 - 263.15) * (1 - np.exp(-alpha_val * t_sim / r_w**2))
    
    axes[1, 0].plot(t_sim/86400, T_center - 273.15, color=color, linewidth=2, 
                    label=f'φ={phi_val}')

axes[1, 0].set_xlabel('时间 (天)', fontsize=11)
axes[1, 0].set_ylabel('中心温度 (°C)', fontsize=11)
axes[1, 0].set_title('不同孔隙率的影响', fontsize=12)
axes[1, 0].legend()
axes[1, 0].grid(True, alpha=0.3)

# 温度梯度分布
T_final = T_history[-1]
dT_dr = np.gradient(T_final, r)

axes[1, 1].plot(r[1:-1], dT_dr[1:-1], 'g-', linewidth=2)
axes[1, 1].set_xlabel('半径 (m)', fontsize=11)
axes[1, 1].set_ylabel('温度梯度 (K/m)', fontsize=11)
axes[1, 1].set_title('温度梯度分布', fontsize=12)
axes[1, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/porous_media_phase_change.png', dpi=150, bbox_inches='tight')
plt.close()
print("多孔介质相变模拟图已保存")

print("\n" + "="*60)
print("仿真3：多孔介质中冰透镜体形成")
print("="*60)

# 模拟冻土中冰透镜体的形成过程

# 土壤参数
porosity_soil = 0.4  # 孔隙率
k_soil_frozen = 2.5  # W/(m·K)，冻结土壤导热系数
k_soil_thawed = 1.5  # W/(m·K)，融化土壤导热系数
rho_soil = 1600  # kg/m³，土壤密度
c_soil = 800  # J/(kg·K)，土壤比热容

# 水分迁移参数
k_hydraulic = 1e-8  # m²，水力传导系数
dP_suction = 1000  # Pa，吸力梯度

# 空间离散
z_soil = np.linspace(0, 2, 100)  # m，深度
dz_soil = z_soil[1] - z_soil[0]

# 初始条件
T_initial_soil = 5 + 273.15  # K，初始温度
T_surface_cooling = -10 + 273.15  # K，地表冷却温度

# 时间参数
t_days = 30
t_soil = np.linspace(0, t_days*86400, 500)  # s
dt_soil = t_soil[1] - t_soil[0]

# 温度场求解（简化的一维模型）
T_soil = np.ones((len(t_soil), len(z_soil))) * T_initial_soil
T_soil[:, 0] = T_surface_cooling  # 地表边界条件

# 冻结锋面位置
freezing_front_position = []

for n in range(1, len(t_soil)):
    for i in range(1, len(z_soil)-1):
        # 根据温度选择导热系数
        if T_soil[n-1, i] < 273.15:
            k_eff = k_soil_frozen
        else:
            k_eff = k_soil_thawed
        
        # 显式差分
        alpha_soil = k_eff / (rho_soil * c_soil)
        T_soil[n, i] = T_soil[n-1, i] + alpha_soil * dt_soil / dz_soil**2 * \
                       (T_soil[n-1, i+1] - 2*T_soil[n-1, i] + T_soil[n-1, i-1])
    
    # 找到冻结锋面位置
    frozen_indices = np.where(T_soil[n] < 273.15)[0]
    if len(frozen_indices) > 0:
        freezing_front_position.append(z_soil[frozen_indices[-1]])
    else:
        freezing_front_position.append(0)

# 冰透镜体厚度估算（简化模型）
# 假设水分向冻结锋面迁移并冻结
ice_lens_thickness = []
for i, z_front in enumerate(freezing_front_position):
    if i > 0 and z_front > freezing_front_position[i-1]:
        # 冰透镜体增长
        v_migration = k_hydraulic * dP_suction / (z_front + 0.01)  # 水分迁移速度
        d_ice = v_migration * dt_soil * porosity_soil  # 冰透镜体厚度增量
        ice_lens_thickness.append(d_ice * 1000)  # mm
    else:
        ice_lens_thickness.append(0)

# 累积冰透镜体厚度
cumulative_ice = np.cumsum(ice_lens_thickness)

print("多孔介质中冰透镜体形成:")
print(f"土壤孔隙率: {porosity_soil}")
print(f"冻结导热系数: {k_soil_frozen} W/(m·K)")
print(f"融化导热系数: {k_soil_thawed} W/(m·K)")
print(f"模拟时间: {t_days} 天")

print(f"\n冻结特征:")
print(f"  最终冻深: {freezing_front_position[-1]:.2f} m")
print(f"  最大冰透镜体厚度: {np.max(cumulative_ice):.2f} mm")
print(f"  平均冻胀速率: {freezing_front_position[-1]/t_days*1000:.2f} mm/天")

# 可视化
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# 温度场演化
time_plot_days = [0, 5, 10, 20, 30]
for t_day in time_plot_days:
    idx = int(t_day * 86400 / dt_soil)
    if idx < len(T_soil):
        axes[0].plot(T_soil[idx] - 273.15, z_soil, linewidth=2, label=f'{t_day} days')

axes[0].axvline(x=0, color='k', linestyle='--', alpha=0.5, label='Freezing point')
axes[0].set_xlabel('Temperature (°C)', fontsize=11)
axes[0].set_ylabel('Depth (m)', fontsize=11)
axes[0].set_title('Temperature Profile Evolution', fontsize=12, fontweight='bold')
axes[0].invert_yaxis()
axes[0].legend(fontsize=9)
axes[0].grid(True, alpha=0.3)

# 冻深和冰透镜体演化
t_days_plot = np.linspace(0, t_days, len(freezing_front_position))
axes[1].plot(t_days_plot, freezing_front_position, 'b-', linewidth=2, label='Freezing Front')
ax2 = axes[1].twinx()
ax2.plot(t_days_plot, cumulative_ice, 'r--', linewidth=2, label='Cumulative Ice')
axes[1].set_xlabel('Time (days)', fontsize=11)
axes[1].set_ylabel('Freezing Depth (m)', fontsize=11, color='b')
ax2.set_ylabel('Cumulative Ice Thickness (mm)', fontsize=11, color='r')
axes[1].set_title('Freezing Front and Ice Lens Formation', fontsize=12, fontweight='bold')
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/ice_lens_formation.png', dpi=150, bbox_inches='tight')
plt.close()

print("\n图3：多孔介质中冰透镜体形成已保存")

print("\n" + "="*60)
print("仿真4：相变储能填充床性能分析")
print("="*60)

# 分析相变材料填充床的储热性能

# 填充床参数
D_bed = 0.5  # m，床层直径
H_bed = 1.0  # m，床层高度
V_bed = np.pi * (D_bed/2)**2 * H_bed  # m³

# 相变材料球参数（石蜡）
D_pcm = 0.02  # m，PCM球直径
V_pcm_single = 4/3 * np.pi * (D_pcm/2)**3  # m³

# 填充参数
void_fraction = 0.4  # 空隙率
N_pcm = int(V_bed * (1 - void_fraction) / V_pcm_single)  # PCM球数量
V_pcm_total = N_pcm * V_pcm_single

# PCM物性
T_m_pcm = 58 + 273.15  # K，熔点
L_pcm = 200e3  # J/kg，潜热
c_pcm = 2500  # J/(kg·K)，比热容
rho_pcm = 900  # kg/m³，密度
k_pcm = 0.2  # W/(m·K)，导热系数

# 传热流体参数（水）
T_inlet_pcm = 70 + 273.15  # K，进口温度
m_dot_pcm = 0.1  # kg/s，质量流量
cp_water_pcm = 4200  # J/(kg·K)

# 储能计算
m_pcm = rho_pcm * V_pcm_total
Q_sensible_solid = m_pcm * c_pcm * (T_m_pcm - (20 + 273.15))  # 固态显热
Q_latent = m_pcm * L_pcm  # 潜热
Q_sensible_liquid = m_pcm * c_pcm * (T_inlet_pcm - T_m_pcm)  # 液态显热
Q_total_pcm = Q_sensible_solid + Q_latent + Q_sensible_liquid

# 传热单元数
A_surface_pcm = N_pcm * np.pi * D_pcm**2  # m²，总表面积
h_conv_pcm = 100  # W/(m²·K)，对流换热系数
NTU_pcm = h_conv_pcm * A_surface_pcm / (m_dot_pcm * cp_water_pcm)
effectiveness_pcm = 1 - np.exp(-NTU_pcm)

# 充电时间估算
Q_max_pcm = m_dot_pcm * cp_water_pcm * (T_inlet_pcm - (20 + 273.15))
t_charge_pcm = Q_total_pcm / (effectiveness_pcm * Q_max_pcm) / 3600  # 小时

print("相变储能填充床性能分析:")
print(f"床层尺寸: D={D_bed*1000:.0f}mm, H={H_bed*1000:.0f}mm")
print(f"PCM球数量: {N_pcm}")
print(f"PCM总体积: {V_pcm_total*1e6:.0f} cm³")
print(f"空隙率: {void_fraction}")

print(f"\n储能容量:")
print(f"  固态显热: {Q_sensible_solid/3600/1000:.2f} kWh")
print(f"  潜热: {Q_latent/3600/1000:.2f} kWh")
print(f"  液态显热: {Q_sensible_liquid/3600/1000:.2f} kWh")
print(f"  总储能: {Q_total_pcm/3600/1000:.2f} kWh")

print(f"\n传热性能:")
print(f"  NTU: {NTU_pcm:.2f}")
print(f"  传热效率: {effectiveness_pcm:.3f}")
print(f"  估算充电时间: {t_charge_pcm:.1f} 小时")

# 不同进口温度的影响
T_inlet_range = np.linspace(60, 80, 10) + 273.15
charge_times = []
for T_in in T_inlet_range:
    Q_sens_liq = m_pcm * c_pcm * (T_in - T_m_pcm)
    Q_tot = Q_sensible_solid + Q_latent + Q_sens_liq
    Q_max_t = m_dot_pcm * cp_water_pcm * (T_in - (20 + 273.15))
    t_ch = Q_tot / (effectiveness_pcm * Q_max_t) / 3600
    charge_times.append(t_ch)

# 可视化
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# 储能容量组成
categories = ['Sensible\n(Solid)', 'Latent', 'Sensible\n(Liquid)']
values = [Q_sensible_solid/3600/1000, Q_latent/3600/1000, Q_sensible_liquid/3600/1000]
colors = ['lightblue', 'orange', 'lightcoral']

axes[0].bar(categories, values, color=colors, alpha=0.7)
axes[0].set_ylabel('Energy Capacity (kWh)', fontsize=11)
axes[0].set_title('PCM Bed Energy Storage Capacity', fontsize=12, fontweight='bold')
axes[0].grid(True, alpha=0.3, axis='y')

# 充电时间vs进口温度
axes[1].plot(T_inlet_range - 273.15, charge_times, 'g-o', linewidth=2, markersize=6)
axes[1].fill_between(T_inlet_range - 273.15, charge_times, alpha=0.3, color='green')
axes[1].axvline(x=70, color='r', linestyle='--', linewidth=2, label='Design: 70°C')
axes[1].set_xlabel('Inlet Temperature (°C)', fontsize=11)
axes[1].set_ylabel('Charging Time (h)', fontsize=11)
axes[1].set_title('Charging Time vs Inlet Temperature', fontsize=12, fontweight='bold')
axes[1].legend(fontsize=10)
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/pcm_packed_bed.png', dpi=150, bbox_inches='tight')
plt.close()

print("\n图4：相变储能填充床性能分析已保存")
print("\n所有仿真完成！")