结构热应力仿真-主题060_地热能源系统-主题060_地热能源系统教程
主题060: 地热能源系统仿真
目录
- 引言
- 地源热泵系统原理
- 热响应测试理论
- 数学模型与数值方法
- Python代码实现
- 案例实战
- 结果分析与讨论
- 进阶应用
- 总结与习题
#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
主题060: 地热能源系统仿真
实例一: 地源热泵热响应测试与长期性能分析
本代码模拟地源热泵系统的热响应测试(TRT)和长期运行性能，
包括钻孔换热器、土壤热传导、热响应测试分析等内容。
"""
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from scipy.special import exp1
from scipy.integrate import odeint
from scipy.optimize import curve_fit
import warnings
warnings.filterwarnings('ignore')
# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False
class GeothermalSystem:
"""地热能源系统模拟器"""
def __init__(self, soil_type='clay'):
"""
初始化地热系统模拟器
参数:
soil_type: 土壤类型 ('clay'-粘土, 'sand'-砂土, 'rock'-岩石, 'gravel'-砾石)
"""
self.soil_type = soil_type
# 根据土壤类型设置参数
if soil_type == 'clay':
# 粘土
self.rho_soil = 1800 # 土壤密度 (kg/m³)
self.cp_soil = 1500 # 土壤比热容 (J/kg·K)
self.k_soil = 1.5 # 土壤导热系数 (W/m·K)
self.alpha_soil = self.k_soil / (self.rho_soil * self.cp_soil) # 热扩散系数
self.T_g = 15.0 # 地温 (°C)
self.porosity = 0.4 # 孔隙率
elif soil_type == 'sand':
# 砂土
self.rho_soil = 1600
self.cp_soil = 1200
self.k_soil = 2.0
self.alpha_soil = self.k_soil / (self.rho_soil * self.cp_soil)
self.T_g = 14.0
self.porosity = 0.35
elif soil_type == 'rock':
# 岩石
self.rho_soil = 2500
self.cp_soil = 1000
self.k_soil = 3.0
self.alpha_soil = self.k_soil / (self.rho_soil * self.cp_soil)
self.T_g = 12.0
self.porosity = 0.05
else: # gravel
# 砾石
self.rho_soil = 2000
self.cp_soil = 1400
self.k_soil = 2.5
self.alpha_soil = self.k_soil / (self.rho_soil * self.cp_soil)
self.T_g = 13.0
self.porosity = 0.3
# 钻孔换热器参数 (默认值)
self.H = 100 # 钻孔深度 (m)
self.r_b = 0.075 # 钻孔半径 (m)
self.r_p = 0.02 # 管道外半径 (m)
self.d = 0.05 # 管间距 (m)
self.k_grout = 1.5 # 回填材料导热系数 (W/m·K)
self.k_pipe = 0.4 # 管道导热系数 (W/m·K)
# 流体参数
self.rho_fluid = 1000 # 流体密度 (kg/m³)
self.cp_fluid = 4180 # 流体比热容 (J/kg·K)
self.k_fluid = 0.6 # 流体导热系数 (W/m·K)
self.mu_fluid = 0.001 # 流体动力粘度 (Pa·s)
def g_function_ils(self, Fo):
"""
无限长线源(ILS) g-函数
参数:
Fo: 傅里叶数 = alpha * t / r_b^2
返回:
g: 无量纲温度响应
"""
if Fo <= 0:
return 0
# ILS解: g = 0.5 * E1(1/(4*Fo))
# 对于大Fo，使用近似: g ≈ 0.5 * (ln(4*Fo) - gamma)
# 其中 gamma ≈ 0.5772 (欧拉常数)
gamma = 0.5772156649
if Fo > 5:
# 大时间近似
g = 0.5 * (np.log(4 * Fo) - gamma)
else:
# 精确解
g = 0.5 * exp1(1 / (4 * Fo))
return g
def g_function_ics(self, Fo):
"""
无限长圆柱源(ICS) g-函数
更精确的模型，考虑钻孔半径
"""
if Fo <= 0:
return 0
gamma = 0.5772156649
if Fo > 20:
# 大时间近似
g = np.log(4 * Fo) - gamma - 1 / (4 * Fo)
else:
# 级数解
g = 0.5 * (exp1(1 / (4 * Fo)) + 0.5 / Fo * np.exp(-1 / (4 * Fo)))
return g
def thermal_response_test(self, Q_heat, t_test, m_flow, T_inlet, model='ILS'):
"""
热响应测试(TRT)仿真
参数:
Q_heat: 加热功率 (W)
t_test: 测试时间数组 (s)
m_flow: 流体质量流量 (kg/s)
T_inlet: 进口流体温度 (°C)
model: 模型类型 ('ILS'-线源, 'ICS'-圆柱源)
返回:
dict: 包含温度响应结果
"""
# 计算傅里叶数
Fo = self.alpha_soil * t_test / self.r_b**2
# 选择g-函数
if model == 'ILS':
g_func = np.array([self.g_function_ils(f) for f in Fo])
else:
g_func = np.array([self.g_function_ics(f) for f in Fo])
# 计算钻孔壁温度 (线源理论)
# T_b = T_g + Q/(2*pi*k_soil*H) * g(Fo)
T_b = self.T_g + Q_heat / (2 * np.pi * self.k_soil * self.H) * g_func
# 计算总热阻
# R_total = R_grout + R_pipe + R_convection
# 对流换热系数 (Dittus-Boelter公式)
Re = 4 * m_flow / (np.pi * (2*self.r_p)**2 * self.mu_fluid) # 雷诺数
Pr = self.cp_fluid * self.mu_fluid / self.k_fluid # 普朗特数
if Re > 2300:
# 湍流
Nu = 0.023 * Re**0.8 * Pr**0.4
else:
# 层流
Nu = 3.66
h_conv = Nu * self.k_fluid / (2 * self.r_p) # 对流换热系数
# 单位长度热阻
R_grout = 1 / (2 * np.pi * self.k_grout) * np.log(self.r_b / (self.r_b - 0.01))
R_pipe = 1 / (2 * np.pi * self.k_pipe) * np.log(self.r_p / (self.r_p - 0.002))
R_conv = 1 / (2 * np.pi * self.r_p * h_conv)
R_total = R_grout + R_pipe + R_conv
# 计算流体平均温度和出口温度
# T_f_avg = T_b + Q * R_total / H
T_f_avg = T_b + Q_heat * R_total / self.H
# 假设对称双U型管，进出口温差
delta_T = Q_heat / (m_flow * self.cp_fluid)
T_outlet = T_f_avg + delta_T / 2
T_inlet_calc = T_f_avg - delta_T / 2
return {
't': t_test,
'T_b': T_b,
'T_f_avg': T_f_avg,
'T_outlet': T_outlet,
'T_inlet': T_inlet_calc,
'Fo': Fo,
'g_func': g_func,
'R_total': R_total
}
def estimate_soil_conductivity(self, Q_heat, t_test, T_f_avg, t_start=None, t_end=None):
"""
从TRT数据估算土壤导热系数
使用线源模型反演:
T_f = T_g + Q/(4*pi*k*H) * (ln(t) + ln(4*alpha/r_b^2) - gamma) + Q*R_b/H
参数:
Q_heat: 加热功率 (W)
t_test: 时间数组 (s)
T_f_avg: 平均流体温度 (°C)
t_start, t_end: 用于拟合的时间范围
返回:
dict: 估算的土壤参数
"""
if t_start is None:
t_start = 3600 # 从1小时开始
if t_end is None:
t_end = t_test[-1]
# 选择拟合数据
mask = (t_test >= t_start) & (t_test <= t_end)
t_fit = t_test[mask]
T_fit = T_f_avg[mask]
# 线源模型: T = m * ln(t) + c
# 其中 m = Q/(4*pi*k*H)
ln_t = np.log(t_fit)
# 线性拟合
coeffs = np.polyfit(ln_t, T_fit, 1)
m = coeffs[0]
c = coeffs[1]
# 估算导热系数
k_est = Q_heat / (4 * np.pi * self.H * m)
# 估算热阻
R_b_est = (c - self.T_g) * self.H / Q_heat - 1 / (4 * np.pi * k_est) *
(np.log(4 * self.alpha_soil / self.r_b**2) - 0.5772)
# 计算拟合质量
T_pred = m * ln_t + c
r_squared = 1 - np.sum((T_fit - T_pred)**2) / np.sum((T_fit - np.mean(T_fit))**2)
return {
'k_soil_estimated': k_est,
'R_b_estimated': R_b_est,
'slope': m,
'intercept': c,
'R_squared': r_squared,
't_fit': t_fit,
'T_fit': T_fit,
'T_pred': T_pred
}
def long_term_performance(self, years, Q_annual, config='balanced'):
"""
长期性能分析
参数:
years: 模拟年数
Q_annual: 年热负荷 (kWh/year)
config: 配置 ('balanced'-平衡, 'heating'-供热为主, 'cooling'-制冷为主)
返回:
dict: 长期性能结果
"""
# 时间数组 (秒)
t_years = np.linspace(0, years, years * 12) # 月度数据
t_seconds = t_years * 365 * 24 * 3600
# 根据配置设置季节性负荷
if config == 'balanced':
# 平衡模式：冬夏负荷相等
Q_seasonal = Q_annual / 8760 * (1 + 0.5 * np.sin(2 * np.pi * t_years))
elif config == 'heating':
# 供热为主：冬季负荷大
Q_seasonal = Q_annual / 8760 * (1 + 0.8 * np.cos(2 * np.pi * t_years))
else: # cooling
# 制冷为主：夏季负荷大
Q_seasonal = Q_annual / 8760 * (1 - 0.8 * np.cos(2 * np.pi * t_years))
# 计算g-函数 (脉冲叠加法简化)
# 使用 Eskilson 的 g-函数近似
g_cumulative = np.zeros_like(t_seconds)
for i, t in enumerate(t_seconds):
if t > 0:
Fo = self.alpha_soil * t / self.r_b**2
g_cumulative[i] = self.g_function_ils(Fo)
# 计算温度变化 (简化模型)
# 考虑热负荷累积效应
delta_T_soil = np.zeros_like(t_seconds)
for i in range(1, len(t_seconds)):
dt = t_seconds[i] - t_seconds[i-1]
Q_avg = (Q_seasonal[i] + Q_seasonal[i-1]) / 2 * 1000 # 转换为W
# 温度变化率
dT_dt = Q_avg / (self.rho_soil * self.cp_soil * np.pi * self.H * (10**2))
delta_T_soil[i] = delta_T_soil[i-1] + dT_dt * dt * 0.01 # 简化系数
# 计算流体温度
T_f = self.T_g + delta_T_soil + Q_seasonal * 1000 / (2 * np.pi * self.k_soil * self.H) * g_cumulative
# 计算性能系数 (COP)
# 简化模型: COP随温差降低
COP_heating = 4.0 - 0.05 * (T_f - self.T_g)
COP_cooling = 5.0 - 0.05 * (self.T_g - T_f)
return {
'years': t_years,
'Q_seasonal': Q_seasonal,
'T_f': T_f,
'delta_T_soil': delta_T_soil,
'COP_heating': COP_heating,
'COP_cooling': COP_cooling,
'g_cumulative': g_cumulative
}
def multiple_boreholes(self, n_boreholes, spacing, Q_total, t_test):
"""
多钻孔系统热干扰分析
参数:
n_boreholes: 钻孔数量
spacing: 钻孔间距 (m)
Q_total: 总加热功率 (W)
t_test: 测试时间 (s)
返回:
dict: 多钻孔系统结果
"""
# 单孔功率
Q_single = Q_total / n_boreholes
# 计算各钻孔位置 (假设线性排列)
positions = np.arange(n_boreholes) * spacing
# 计算中心钻孔的温度 (受所有钻孔影响)
center_idx = n_boreholes // 2
# 自响应
Fo_self = self.alpha_soil * t_test / self.r_b**2
g_self = self.g_function_ils(Fo_self)
# 相邻钻孔干扰
g_neighbors = 0
for i in range(n_boreholes):
if i != center_idx:
r_ij = abs(positions[i] - positions[center_idx])
Fo_ij = self.alpha_soil * t_test / r_ij**2
g_ij = self.g_function_ils(Fo_ij)
g_neighbors += g_ij
# 总g-函数
g_total = g_self + g_neighbors
# 中心钻孔温度
T_center = self.T_g + Q_single / (2 * np.pi * self.k_soil * self.H) * g_total
# 单孔温度 (无干扰)
T_single = self.T_g + Q_single / (2 * np.pi * self.k_soil * self.H) * g_self
# 温度升高比
temp_ratio = (T_center - self.T_g) / (T_single - self.T_g)
return {
'T_center': T_center,
'T_single': T_single,
'temp_ratio': temp_ratio,
'g_total': g_total,
'g_self': g_self,
'positions': positions
}
def visualize_trt_results():
"""可视化热响应测试结果"""
print("="*60)
print("案例1: 热响应测试(TRT)仿真")
print("="*60)
# 创建地热系统
geo = GeothermalSystem(soil_type='clay')
# TRT参数
Q_heat = 5000 # 加热功率 5 kW
t_test = np.logspace(2, 6, 200) # 时间从100秒到10万秒
m_flow = 0.3 # 流量 0.3 kg/s
T_inlet = 15.0 # 进口温度
# ILS模型
result_ils = geo.thermal_response_test(Q_heat, t_test, m_flow, T_inlet, model='ILS')
# ICS模型
result_ics = geo.thermal_response_test(Q_heat, t_test, m_flow, T_inlet, model='ICS')
# 绘制结果
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# 图1: 温度响应
ax1 = axes[0, 0]
ax1.semilogx(result_ils['t']/3600, result_ils['T_f_avg'], 'b-', linewidth=2, label='ILS模型')
ax1.semilogx(result_ics['t']/3600, result_ics['T_f_avg'], 'r--', linewidth=2, label='ICS模型')
ax1.axhline(y=geo.T_g, color='g', linestyle=':', label=f'初始地温 ({geo.T_g}°C)')
ax1.set_xlabel('时间 (小时)', fontsize=11)
ax1.set_ylabel('平均流体温度 (°C)', fontsize=11)
ax1.set_title('热响应测试: 流体温度变化', fontsize=12, fontweight='bold')
ax1.legend()
ax1.grid(True, alpha=0.3)
# 图2: g-函数
ax2 = axes[0, 1]
ax2.semilogx(result_ils['Fo'], result_ils['g_func'], 'b-', linewidth=2, label='ILS')
ax2.semilogx(result_ics['Fo'], result_ics['g_func'], 'r--', linewidth=2, label='ICS')
ax2.set_xlabel('傅里叶数 Fo', fontsize=11)
ax2.set_ylabel('g-函数', fontsize=11)
ax2.set_title('无量纲g-函数对比', fontsize=12, fontweight='bold')
ax2.legend()
ax2.grid(True, alpha=0.3)
# 图3: 进出口温度
ax3 = axes[1, 0]
ax3.semilogx(result_ils['t']/3600, result_ils['T_inlet'], 'b-', linewidth=2, label='进口')
ax3.semilogx(result_ils['t']/3600, result_ils['T_outlet'], 'r-', linewidth=2, label='出口')
ax3.fill_between(result_ils['t']/3600, result_ils['T_inlet'], result_ils['T_outlet'],
alpha=0.3, color='yellow', label='温差')
ax3.set_xlabel('时间 (小时)', fontsize=11)
ax3.set_ylabel('温度 (°C)', fontsize=11)
ax3.set_title('进出口温度与温差', fontsize=12, fontweight='bold')
ax3.legend()
ax3.grid(True, alpha=0.3)
# 图4: 温度-对数时间曲线 (用于估算导热系数)
ax4 = axes[1, 1]
ln_t = np.log(result_ils['t'])
ax4.plot(ln_t, result_ils['T_f_avg'], 'b-', linewidth=2, label='TRT数据')
# 线性拟合 (从1小时后)
mask = result_ils['t'] >= 3600
if np.any(mask):
fit_result = geo.estimate_soil_conductivity(Q_heat, result_ils['t'],
result_ils['T_f_avg'])
ax4.plot(np.log(fit_result['t_fit']), fit_result['T_pred'], 'r--',
linewidth=2, label=f'线性拟合 (R²={fit_result["R_squared"]:.4f})')
ax4.text(0.05, 0.95, f'估算导热系数: {fit_result["k_soil_estimated"]:.2f} W/m·Kn'
f'实际导热系数: {geo.k_soil:.2f} W/m·Kn'
f'估算钻孔热阻: {fit_result["R_b_estimated"]:.4f} m·K/W',
transform=ax4.transAxes, fontsize=10, verticalalignment='top',
bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))
ax4.set_xlabel('ln(t) [t单位为秒]', fontsize=11)
ax4.set_ylabel('平均流体温度 (°C)', fontsize=11)
ax4.set_title('温度-对数时间曲线 (导热系数估算)', fontsize=12, fontweight='bold')
ax4.legend()
ax4.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('热响应测试仿真.png', dpi=150, bbox_inches='tight')
print(" 图表已保存: 热响应测试仿真.png")
plt.close()
print(f"nTRT结果:")
print(f" 50小时后流体温度: {result_ils['T_f_avg'][-1]:.2f}°C")
print(f" 温度升高: {result_ils['T_f_avg'][-1] - geo.T_g:.2f}°C")
print(f" 估算导热系数: {fit_result['k_soil_estimated']:.2f} W/m·K")
def visualize_soil_comparison():
"""不同土壤类型对比"""
print("n" + "="*60)
print("案例2: 不同土壤类型热响应对比")
print("="*60)
soil_types = ['clay', 'sand', 'rock', 'gravel']
soil_names = ['粘土', '砂土', '岩石', '砾石']
colors = ['brown', 'orange', 'gray', 'yellowgreen']
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# TRT参数
Q_heat = 5000
t_test = np.logspace(2, 6, 200)
m_flow = 0.3
T_inlet = 15.0
# 图1: 温度响应对比
ax1 = axes[0]
for soil_type, name, color in zip(soil_types, soil_names, colors):
geo = GeothermalSystem(soil_type=soil_type)
result = geo.thermal_response_test(Q_heat, t_test, m_flow, T_inlet)
ax1.semilogx(result['t']/3600, result['T_f_avg'], color=color,
linewidth=2, label=f'{name} (k={geo.k_soil} W/m·K)')
ax1.set_xlabel('时间 (小时)', fontsize=11)
ax1.set_ylabel('平均流体温度 (°C)', fontsize=11)
ax1.set_title('不同土壤类型的TRT响应', fontsize=12, fontweight='bold')
ax1.legend()
ax1.grid(True, alpha=0.3)
# 图2: 物性参数对比
ax2 = axes[1]
k_values = []
alpha_values = []
T_g_values = []
for soil_type in soil_types:
geo = GeothermalSystem(soil_type=soil_type)
k_values.append(geo.k_soil)
alpha_values.append(geo.alpha_soil * 1e6) # 转换为10^-6 m²/s
T_g_values.append(geo.T_g)
x = np.arange(len(soil_types))
width = 0.25
ax2_twin = ax2.twinx()
ax2_twin2 = ax2.twinx()
ax2_twin2.spines['right'].set_position(('outward', 60))
bars1 = ax2.bar(x - width, k_values, width, label='导热系数', color='steelblue')
bars2 = ax2_twin.bar(x, alpha_values, width, label='热扩散系数', color='coral')
bars3 = ax2_twin2.bar(x + width, T_g_values, width, label='地温', color='green')
ax2.set_xlabel('土壤类型', fontsize=11)
ax2.set_ylabel('导热系数 (W/m·K)', fontsize=11, color='steelblue')
ax2_twin.set_ylabel('热扩散系数 (10⁻⁶ m²/s)', fontsize=11, color='coral')
ax2_twin2.set_ylabel('地温 (°C)', fontsize=11, color='green')
ax2.set_xticks(x)
ax2.set_xticklabels(soil_names)
ax2.set_title('土壤热物性参数对比', fontsize=12, fontweight='bold')
# 添加图例
lines1, labels1 = ax2.get_legend_handles_labels()
lines2, labels2 = ax2_twin.get_legend_handles_labels()
lines3, labels3 = ax2_twin2.get_legend_handles_labels()
ax2.legend(lines1 + lines2 + lines3, labels1 + labels2 + labels3, loc='upper right')
plt.tight_layout()
plt.savefig('土壤类型对比.png', dpi=150, bbox_inches='tight')
print(" 图表已保存: 土壤类型对比.png")
plt.close()
def visualize_long_term():
"""长期性能分析"""
print("n" + "="*60)
print("案例3: 长期运行性能分析")
print("="*60)
geo = GeothermalSystem(soil_type='clay')
# 模拟10年
years = 10
Q_annual = 15000 # 年热负荷 15 MWh
# 不同配置
configs = ['balanced', 'heating', 'cooling']
config_names = ['平衡模式', '供热模式', '制冷模式']
colors = ['green', 'red', 'blue']
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
for config, name, color in zip(configs, config_names, colors):
result = geo.long_term_performance(years, Q_annual, config)
# 图1: 流体温度
ax1 = axes[0, 0]
ax1.plot(result['years'], result['T_f'], color=color, linewidth=2, label=name)
ax1.axhline(y=geo.T_g, color='k', linestyle='--', alpha=0.5)
# 图2: 土壤温度变化
ax2 = axes[0, 1]
ax2.plot(result['years'], result['delta_T_soil'], color=color, linewidth=2, label=name)
# 图3: 季节性负荷
ax3 = axes[1, 0]
ax3.plot(result['years'], result['Q_seasonal']/1000, color=color, linewidth=2, label=name)
# 图4: COP
ax4 = axes[1, 1]
if config == 'cooling':
ax4.plot(result['years'], result['COP_cooling'], color=color, linewidth=2, label=name)
else:
ax4.plot(result['years'], result['COP_heating'], color=color, linewidth=2, label=name)
# 设置标签
axes[0, 0].set_xlabel('时间 (年)', fontsize=11)
axes[0, 0].set_ylabel('流体温度 (°C)', fontsize=11)
axes[0, 0].set_title('长期流体温度变化', fontsize=12, fontweight='bold')
axes[0, 0].legend()
axes[0, 0].grid(True, alpha=0.3)
axes[0, 1].set_xlabel('时间 (年)', fontsize=11)
axes[0, 1].set_ylabel('土壤温升 (°C)', fontsize=11)
axes[0, 1].set_title('土壤热累积效应', fontsize=12, fontweight='bold')
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3)
axes[1, 0].set_xlabel('时间 (年)', fontsize=11)
axes[1, 0].set_ylabel('热负荷 (kW)', fontsize=11)
axes[1, 0].set_title('季节性热负荷', fontsize=12, fontweight='bold')
axes[1, 0].legend()
axes[1, 0].grid(True, alpha=0.3)
axes[1, 1].set_xlabel('时间 (年)', fontsize=11)
axes[1, 1].set_ylabel('性能系数 COP', fontsize=11)
axes[1, 1].set_title('系统COP变化', fontsize=12, fontweight='bold')
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('长期性能分析.png', dpi=150, bbox_inches='tight')
print(" 图表已保存: 长期性能分析.png")
plt.close()
def visualize_multiple_boreholes():
"""多钻孔热干扰分析"""
print("n" + "="*60)
print("案例4: 多钻孔系统热干扰分析")
print("="*60)
geo = GeothermalSystem(soil_type='clay')
# 分析不同钻孔数量和间距
n_boreholes_list = [1, 3, 5, 7, 9]
spacings = [3, 5, 7, 10] # 米
t_test = 3600 * 24 * 30 # 30天
Q_total = 50000 # 总功率 50 kW
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# 图1: 不同钻孔数量的温度升高比
ax1 = axes[0]
for spacing in spacings:
temp_ratios = []
for n in n_boreholes_list:
result = geo.multiple_boreholes(n, spacing, Q_total, t_test)
temp_ratios.append(result['temp_ratio'])
ax1.plot(n_boreholes_list, temp_ratios, 'o-', linewidth=2,
markersize=8, label=f'间距 {spacing}m')
ax1.axhline(y=1.0, color='k', linestyle='--', alpha=0.5, label='无干扰')
ax1.set_xlabel('钻孔数量', fontsize=11)
ax1.set_ylabel('温度升高比', fontsize=11)
ax1.set_title('热干扰效应 (运行30天)', fontsize=12, fontweight='bold')
ax1.legend()
ax1.grid(True, alpha=0.3)
# 图2: 不同时间的热干扰
ax2 = axes[1]
n_boreholes = 5
spacing = 5
times = np.logspace(3, 7, 50) # 从1000秒到约100天
temp_ratios_time = []
for t in times:
result = geo.multiple_boreholes(n_boreholes, spacing, Q_total, t)
temp_ratios_time.append(result['temp_ratio'])
ax2.semilogx(times/3600/24, temp_ratios_time, 'b-', linewidth=2)
ax2.axhline(y=1.0, color='k', linestyle='--', alpha=0.5)
ax2.set_xlabel('时间 (天)', fontsize=11)
ax2.set_ylabel('温度升高比', fontsize=11)
ax2.set_title(f'{n_boreholes}个钻孔的热干扰随时间变化 (间距{spacing}m)',
fontsize=12, fontweight='bold')
ax2.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('多钻孔热干扰分析.png', dpi=150, bbox_inches='tight')
print(" 图表已保存: 多钻孔热干扰分析.png")
plt.close()
def visualize_borehole_design():
"""钻孔设计参数分析"""
print("n" + "="*60)
print("案例5: 钻孔设计参数优化")
print("="*60)
geo = GeothermalSystem(soil_type='clay')
# 参数范围
H_range = np.linspace(50, 200, 20) # 深度 50-200m
r_b_range = np.linspace(0.05, 0.15, 20) # 半径 5-15cm
k_grout_range = np.linspace(0.5, 3.0, 20) # 回填导热系数
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# 固定参数
Q_heat = 5000
t_test = 3600 * 50 # 50小时
m_flow = 0.3
# 图1: 深度影响
ax1 = axes[0, 0]
T_f_values = []
for H in H_range:
geo.H = H
result = geo.thermal_response_test(Q_heat, np.array([t_test]), m_flow, 15.0)
T_f_values.append(result['T_f_avg'][0])
ax1.plot(H_range, T_f_values, 'b-', linewidth=2)
ax1.set_xlabel('钻孔深度 (m)', fontsize=11)
ax1.set_ylabel('流体温度 (°C)', fontsize=11)
ax1.set_title('钻孔深度对温度的影响', fontsize=12, fontweight='bold')
ax1.grid(True, alpha=0.3)
# 图2: 半径影响
ax2 = axes[0, 1]
T_f_values = []
geo.H = 100 # 恢复默认
for r_b in r_b_range:
geo.r_b = r_b
result = geo.thermal_response_test(Q_heat, np.array([t_test]), m_flow, 15.0)
T_f_values.append(result['T_f_avg'][0])
ax2.plot(r_b_range*100, T_f_values, 'r-', linewidth=2)
ax2.set_xlabel('钻孔半径 (cm)', fontsize=11)
ax2.set_ylabel('流体温度 (°C)', fontsize=11)
ax2.set_title('钻孔半径对温度的影响', fontsize=12, fontweight='bold')
ax2.grid(True, alpha=0.3)
geo.r_b = 0.075 # 恢复默认
# 图3: 回填材料影响
ax3 = axes[1, 0]
T_f_values = []
R_b_values = []
for k_grout in k_grout_range:
geo.k_grout = k_grout
result = geo.thermal_response_test(Q_heat, np.array([t_test]), m_flow, 15.0)
T_f_values.append(result['T_f_avg'][0])
R_b_values.append(result['R_total'])
ax3.plot(k_grout_range, T_f_values, 'g-', linewidth=2)
ax3.set_xlabel('回填材料导热系数 (W/m·K)', fontsize=11)
ax3.set_ylabel('流体温度 (°C)', fontsize=11)
ax3.set_title('回填材料对温度的影响', fontsize=12, fontweight='bold')
ax3.grid(True, alpha=0.3)
# 图4: 热阻分析
ax4 = axes[1, 1]
ax4.plot(k_grout_range, R_b_values, 'm-', linewidth=2)
ax4.set_xlabel('回填材料导热系数 (W/m·K)', fontsize=11)
ax4.set_ylabel('总热阻 (m·K/W)', fontsize=11)
ax4.set_title('回填材料导热系数与热阻关系', fontsize=12, fontweight='bold')
ax4.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('钻孔设计参数分析.png', dpi=150, bbox_inches='tight')
print(" 图表已保存: 钻孔设计参数分析.png")
plt.close()
if __name__ == "__main__":
print("n" + "="*60)
print("主题060: 地热能源系统仿真")
print("地源热泵热响应测试与长期性能分析")
print("="*60 + "n")
# 运行所有案例
visualize_trt_results()
visualize_soil_comparison()
visualize_long_term()
visualize_multiple_boreholes()
visualize_borehole_design()
print("n" + "="*60)
print("仿真完成!")
print("="*60)
print("n生成的文件:")
print(" 1. 热响应测试仿真.png")
print(" 2. 土壤类型对比.png")
print(" 3. 长期性能分析.png")
print(" 4. 多钻孔热干扰分析.png")
print(" 5. 钻孔设计参数分析.png")
© 版权声明
文章版权归作者所有，未经允许请勿转载。
