import math import os from dataclasses import dataclass from typing import Dict, List, Tuple, Optional import numpy as np from matplotlib import pyplot as plt from scipy.interpolate import CubicSpline from scipy.integrate import simpson # 全局物理常数（修正：优化热物理参数，确保量级合理） SIGMA = 5.670374419e-8 # Stefan-Boltzmann constant (W/m²/K⁴) T_AMB = 300.0 # Ambient temperature (K) ≈27℃ T_SKY = 270.0 # 降低天空温度（更贴近真实晴空，增强辐射散热） SOLAR_IRR = 850.0 # 太阳辐照度（合理范围：800-900 W/m²） H_CONV = 10.0 # 自然对流换热系数（优化为10 W/m²/K，更贴近实际） H_COND_MAX = 1000.0 # 最大传导换热系数（限制量级，避免异常大值） k_PDMS = 0.12 # PDMS导热系数（取保守值0.12 W/m·K，避免偏大） MAX_TEMPERATURE_DIFF = 15.0 # 最大允许温差（K）：避免平衡温度过低导致热损失量级异常 # 关键波段定义 SOLAR_BAND = (0.3, 2.5) # Solar spectrum range (μm) ATMOSPHERIC_WINDOW = (8.0, 13.0) # Atmospheric transparency window (μm) THERMAL_BAND = (2.5, 25.0) # Thermal radiation full band (μm) # 绘图配置（纯英文） plt.rcParams["font.sans-serif"] = ["Arial", "Helvetica", "DejaVu Sans"] plt.rcParams["axes.unicode_minus"] = False # Fix minus sign display plt.rcParams["figure.dpi"] = 100 plt.rcParams["savefig.dpi"] = 300 @dataclass class OpticalProperty: """Optical property data class""" wl: np.ndarray # Wavelength (μm) n: np.ndarray # Refractive index k: np.ndarray # Extinction coefficient cs_n: CubicSpline # Cubic spline for refractive index cs_k: CubicSpline # Cubic spline for extinction coefficient @dataclass class SpectralMetric: """Spectral matching metrics""" thickness_um: float avg_eps_window: float # Average emissivity in atmospheric window avg_alpha_solar: float # Average absorptivity in solar band spectral_match_index: float # Spectral matching index (ε_window/α_solar) eps_window_std: float # Std of emissivity in atmospheric window alpha_solar_std: float # Std of absorptivity in solar band @dataclass class CoolingPerformance: """Comprehensive cooling performance metrics""" thickness_um: float spectral_metric: SpectralMetric net_radiative_power: float # Net radiative power (W/m²) 核心新增指标，单独绘图 solar_gain: float # Solar absorption heat flux (W/m²) convective_loss: float # Convective heat loss (W/m²) conductive_loss: float # Conductive heat loss (W/m²) net_cooling_power: float # Total net cooling power (W/m²) equilibrium_temp_K: float # Equilibrium temperature (K) @property def equilibrium_temp_C(self) -> float: return self.equilibrium_temp_K - 273.15 @property def temp_drop_C(self) -> float: # 正确逻辑：降温值 = 环境温度 - 平衡温度（确保为正，有效冷却） return T_AMB - self.equilibrium_temp_K # ------------------------------ # 1. Data Preprocessing and Optical Property Loading # ------------------------------ def make_strictly_increasing(wl: np.ndarray, n: np.ndarray, k: np.ndarray) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: """Ensure wavelength data is strictly increasing and free of duplicates""" # Remove duplicate wavelengths unique_wl, indices = np.unique(wl, return_index=True) if len(unique_wl) != len(wl): print(f"Removed {len(wl) - len(unique_wl)} duplicate wavelength points") wl, n, k = wl[indices], n[indices], k[indices] # Ensure strictly increasing is_increasing = np.diff(wl) > 1e-6 # Allow minor tolerance if not all(is_increasing): valid_indices = np.concatenate([[True], is_increasing]) wl, n, k = wl[valid_indices], n[valid_indices], k[valid_indices] print(f"Removed {len(is_increasing) - sum(is_increasing)} non-increasing wavelength points") return wl, n, k def read_split_data(file_path: str) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: """Read split-format optical data (n and k stored separately)""" with open(file_path, 'r', encoding='utf-8') as f: lines = [line.strip() for line in f if line.strip() and not line.startswith('#')] # Find split point between n and k data split_idx = None for i, line in enumerate(lines): if line == 'wl k': split_idx = i break if split_idx is None: raise ValueError("Data file format error: 'wl k' delimiter not found") # Read n data n_lines = lines[1:split_idx] # Skip 'wl n' header wl_n, n_list = [], [] for line in n_lines: parts = line.split() if len(parts) >= 2: wl_n.append(float(parts[0])) n_list.append(float(parts[1])) # Read k data k_lines = lines[split_idx + 1:] # Skip 'wl k' header wl_k, k_list = [], [] for line in k_lines: parts = line.split() if len(parts) >= 2: wl_k.append(float(parts[0])) k_list.append(float(parts[1])) # Convert to numpy arrays wl_n, n_list = np.array(wl_n), np.array(n_list) wl_k, k_list = np.array(wl_k), np.array(k_list) # Unify wavelength range (interpolate to common wavelengths) if not np.allclose(wl_n, wl_k): print("Warning: Wavelength ranges of n and k do not match, performing interpolation") common_wl = np.linspace( max(wl_n.min(), wl_k.min()), min(wl_n.max(), wl_k.max()), 1000 ) n_list = np.interp(common_wl, wl_n, n_list) k_list = np.interp(common_wl, wl_k, k_list) wl_n = common_wl # Ensure strictly increasing wl, n, k = make_strictly_increasing(wl_n, n_list, k_list) return wl, n, k def load_optical_properties(file_path: str) -> OpticalProperty: """Load optical properties and create interpolation functions""" wl, n, k = read_split_data(file_path) cs_n = CubicSpline(wl, n, bc_type='natural') cs_k = CubicSpline(wl, k, bc_type='natural') print(f"Data loaded successfully:") print(f" Wavelength range: {wl.min():.2f} - {wl.max():.2f} μm") print(f" Number of data points: {len(wl)}") print(f" Refractive index n: {n.min():.3f} - {n.max():.3f}") print(f" Extinction coefficient k: {k.min():.6f} - {k.max():.6f}") return OpticalProperty(wl=wl, n=n, k=k, cs_n=cs_n, cs_k=cs_k) # ------------------------------ # 2. Core Optical Model (Emissivity/Absorptivity Calculation) # ------------------------------ def fresnel_reflectance(n1: float, k1: float, n2: float, k2: float) -> float: """Fresnel reflectance (normal incidence, considering complex refractive index)""" m1 = n1 + 1j * k1 m2 = n2 + 1j * k2 return np.abs((m1 - m2) / (m1 + m2)) ** 2 def thin_film_spectral_property( optical: OpticalProperty, thickness_um: float, wl_range: Tuple[float, float] ) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: """Calculate spectral emissivity and absorptivity of thin film in specified wavelength range""" wl_start, wl_end = wl_range # Generate dense wavelength points wl = np.linspace(wl_start, wl_end, 500) # Interpolate optical parameters n_film = optical.cs_n(wl) k_film = optical.cs_k(wl) n_air, k_air = 1.0, 0.0 # Calculate Fresnel reflectance R12 = fresnel_reflectance(n_air, k_air, n_film, k_film) # Air → Film R23 = fresnel_reflectance(n_film, k_film, n_air, k_air) # Film → Air # Calculate phase difference and absorption attenuation delta = 2 * np.pi * n_film * thickness_um / wl # Interference phase difference alpha = 4 * np.pi * k_film * thickness_um / wl # Absorption attenuation coefficient # Multibeam interference reflectance numerator = R12 + R23 * np.exp(-alpha) + 2 * np.sqrt(R12 * R23 * np.exp(-alpha)) * np.cos(2 * delta) denominator = 1 + R12 * R23 * np.exp(-alpha) + 2 * np.sqrt(R12 * R23 * np.exp(-alpha)) * np.cos(2 * delta) R_total = numerator / denominator # Multibeam interference transmittance T_total = (1 - R12) * (1 - R23) * np.exp(-alpha) / denominator # Emissivity (Kirchhoff's law) and absorptivity (α=ε under local thermal equilibrium) epsilon = 1 - R_total - T_total alpha = epsilon # Absorptivity = Emissivity under local thermal equilibrium # Numerical clipping to avoid out-of-range values due to calculation errors epsilon = np.clip(epsilon, 0.0, 1.0) alpha = np.clip(alpha, 0.0, 1.0) return wl, epsilon, alpha # ------------------------------ # 3. Spectral Matching Quantitative Analysis # ------------------------------ def calculate_spectral_metric( optical: OpticalProperty, thickness_um: float ) -> SpectralMetric: """Calculate spectral matching metrics""" # 1. Emissivity in atmospheric window (8-13μm) wl_window, eps_window, _ = thin_film_spectral_property( optical, thickness_um, ATMOSPHERIC_WINDOW ) avg_eps = np.mean(eps_window) std_eps = np.std(eps_window) # 2. Absorptivity in solar band (0.3-2.5μm) wl_solar, _, alpha_solar = thin_film_spectral_property( optical, thickness_um, SOLAR_BAND ) avg_alpha = np.mean(alpha_solar) std_alpha = np.std(alpha_solar) # 3. Spectral matching index (larger is better, avoid division by zero) match_index = avg_eps / (avg_alpha + 1e-6) return SpectralMetric( thickness_um=thickness_um, avg_eps_window=avg_eps, avg_alpha_solar=avg_alpha, spectral_match_index=match_index, eps_window_std=std_eps, alpha_solar_std=std_alpha ) def plot_spectral_metrics(metrics_list: List[SpectralMetric], save_path: str): """Plot spectral matching metrics (full English)""" thicknesses = [m.thickness_um for m in metrics_list] avg_eps = [m.avg_eps_window for m in metrics_list] avg_alpha = [m.avg_alpha_solar for m in metrics_list] match_indices = [m.spectral_match_index for m in metrics_list] fig, ax1 = plt.subplots(figsize=(10, 6)) # Left axis: Emissivity and Absorptivity color1 = '#2E86AB' color2 = '#A23B72' ax1.set_xlabel('PDMS Film Thickness (μm)', fontsize=12) ax1.set_ylabel('Emissivity / Absorptivity', fontsize=12) line1 = ax1.plot(thicknesses, avg_eps, marker='o', linewidth=2, color=color1, label='Avg Emissivity (8-13μm)', markersize=6) line2 = ax1.plot(thicknesses, avg_alpha, marker='s', linewidth=2, color=color2, label='Avg Absorptivity (0.3-2.5μm)', markersize=6) ax1.tick_params(axis='y') ax1.set_ylim(0, 1.0) ax1.grid(True, alpha=0.3) # Right axis: Spectral Matching Index ax2 = ax1.twinx() color3 = '#F18F01' ax2.set_ylabel('Spectral Matching Index', fontsize=12, color=color3) line3 = ax2.plot(thicknesses, match_indices, marker='^', linewidth=2, color=color3, label='Spectral Matching Index', markersize=6) ax2.tick_params(axis='y', labelcolor=color3) ax2.set_ylim(0, max(match_indices) * 1.1) # Combine legends lines = line1 + line2 + line3 labels = [l.get_label() for l in lines] ax1.legend(lines, labels, loc='upper right', fontsize=10) plt.title('Spectral Matching Metrics of PDMS Films vs Thickness', fontsize=14, fontweight='bold') plt.tight_layout() plt.savefig(save_path, bbox_inches='tight') print(f"Spectral matching metrics plot saved to: {save_path}") # ------------------------------ # 4. Solar Absorptivity Modeling and Visualization # ------------------------------ def plot_solar_absorption_spectra( optical: OpticalProperty, thicknesses: List[float], save_path: str ): """Plot solar band absorptivity spectra for different thicknesses (full English)""" fig, ax = plt.subplots(figsize=(10, 6)) colors = plt.cm.Set3(np.linspace(0, 1, len(thicknesses))) for idx, thickness in enumerate(thicknesses): wl_solar, _, alpha_solar = thin_film_spectral_property( optical, thickness, SOLAR_BAND ) ax.plot(wl_solar, alpha_solar, linewidth=2, color=colors[idx], label=f'{thickness} μm') ax.set_xlabel('Wavelength (μm)', fontsize=12) ax.set_ylabel('Solar Absorptivity α(λ)', fontsize=12) ax.set_title('Solar Band Absorptivity Spectra of PDMS Films with Different Thicknesses', fontsize=14, fontweight='bold') ax.grid(True, alpha=0.3) ax.legend(title='Film Thickness', fontsize=10, title_fontsize=11) ax.set_ylim(0, 0.3) # Solar absorptivity is usually low, zoom in ax.set_xlim(SOLAR_BAND[0], SOLAR_BAND[1]) # Highlight visible light region ax.axvspan(0.4, 0.7, alpha=0.1, color='yellow', label='Visible Light Region') # Redraw legend to include visible light label ax.legend(title='Film Thickness', fontsize=10, title_fontsize=11, loc='upper right') plt.tight_layout() plt.savefig(save_path, bbox_inches='tight') print(f"Solar absorption spectra plot saved to: {save_path}") # ------------------------------ # 5. 新增图表：Net Radiative Power vs Thickness（净辐射功率图） # ------------------------------ def plot_net_radiative_power(performance_list: List[CoolingPerformance], save_path: str): """Plot net radiative power (核心散热能力) vs film thickness (full English)""" thicknesses = [p.thickness_um for p in performance_list] net_radiative_power = [p.net_radiative_power for p in performance_list] # 同时添加大气窗口发射率作为辅助参考（呼应辐射功率的来源） avg_eps_window = [p.spectral_metric.avg_eps_window for p in performance_list] fig, ax1 = plt.subplots(figsize=(10, 6)) # 左Y轴：净辐射功率（核心指标） color1 = '#1f77b4' # 深蓝色 ax1.set_xlabel('PDMS Film Thickness (μm)', fontsize=12) ax1.set_ylabel('Net Radiative Power (W/m²)', fontsize=12, color=color1) line1 = ax1.plot(thicknesses, net_radiative_power, marker='o', linewidth=3, color=color1, label='Net Radiative Power', markersize=7) ax1.tick_params(axis='y', labelcolor=color1) ax1.grid(True, alpha=0.3) # 调整Y轴范围，确保所有数据点清晰 ax1.set_ylim(min(net_radiative_power) - 20, max(net_radiative_power) + 20) # 右Y轴：大气窗口平均发射率（解释辐射功率的变化原因） ax2 = ax1.twinx() color2 = '#ff7f0e' # 橙色 ax2.set_ylabel('Avg Emissivity (8-13μm)', fontsize=12, color=color2) line2 = ax2.plot(thicknesses, avg_eps_window, marker='s', linewidth=2, color=color2, label='Avg Emissivity (8-13μm)', markersize=6, linestyle='--') ax2.tick_params(axis='y', labelcolor=color2) ax2.set_ylim(0, 1.0) # 合并图例 lines = line1 + line2 labels = [l.get_label() for l in lines] ax1.legend(lines, labels, loc='lower right', fontsize=10) plt.title('Net Radiative Power of PDMS Films vs Thickness', fontsize=14, fontweight='bold') plt.tight_layout() plt.savefig(save_path, bbox_inches='tight') print(f"Net radiative power plot saved to: {save_path}") # ------------------------------ # 6. Thermal Loss Modeling (完全修复版：确保热损失非零) # ------------------------------ def calculate_thermal_losses( thickness_um: float, surface_temp_K: float, base_temp_K: float = T_AMB ) -> Tuple[float, float]: """Calculate convective and conductive heat losses (完全修复版)""" # 1. 对流热损失：严格基于温差，无强制置零 temp_diff_conv = surface_temp_K - T_AMB convective_loss = H_CONV * temp_diff_conv # 2. 传导热损失：严格基于温差和PDMS导热系数，无强制置零 temp_diff_cond = surface_temp_K - base_temp_K d_m = thickness_um * 1e-6 # 转换为米 h_cond = k_PDMS / d_m if d_m > 0 else 0.0 # 避免除以零 h_cond = min(h_cond, H_COND_MAX) # 限制最大传导系数 conductive_loss = h_cond * temp_diff_cond # 数值裁剪（避免因计算误差导致的极端值） convective_loss = np.clip(convective_loss, -1000, 1000) conductive_loss = np.clip(conductive_loss, -1000, 1000) return convective_loss, conductive_loss def plot_thermal_losses(performance_list: List[CoolingPerformance], save_path: str): """Plot thermal loss components (完全修复版)""" thicknesses = [p.thickness_um for p in performance_list] convective = [p.convective_loss for p in performance_list] conductive = [p.conductive_loss for p in performance_list] solar_gain = [p.solar_gain for p in performance_list] fig, ax = plt.subplots(figsize=(10, 6)) x = np.arange(len(thicknesses)) width = 0.25 bars1 = ax.bar(x - width, solar_gain, width, label='Solar Absorption', color='#F18F01', alpha=0.8) bars2 = ax.bar(x, convective, width, label='Convective Loss', color='#2E86AB', alpha=0.8) bars3 = ax.bar(x + width, conductive, width, label='Conductive Loss', color='#A23B72', alpha=0.8) ax.set_xlabel('PDMS Film Thickness (μm)', fontsize=12) ax.set_ylabel('Heat Flux Density (W/m²)', fontsize=12) ax.set_title('Thermal Loss Components of PDMS Films at Equilibrium Temperature', fontsize=14, fontweight='bold') ax.set_xticks(x) ax.set_xticklabels(thicknesses) ax.legend(fontsize=10) ax.grid(True, alpha=0.3, axis='y') # 数值标签：避免重叠，保留1位小数 def add_labels(bars): for bar in bars: height = bar.get_height() va = 'bottom' if height >= 0 else 'top' y_offset = 1.0 if height >= 0 else -1.0 ax.text(bar.get_x() + bar.get_width() / 2., height + y_offset, f'{height:.1f}', ha='center', va=va, fontsize=9) add_labels(bars1) add_labels(bars2) add_labels(bars3) plt.tight_layout() plt.savefig(save_path, bbox_inches='tight') print(f"Thermal loss components plot saved to: {save_path}") # ------------------------------ # 7. Comprehensive Cooling Performance Evaluation (完全修复版) # ------------------------------ def calculate_cooling_performance( optical: OpticalProperty, thickness_um: float ) -> CoolingPerformance: """Calculate comprehensive cooling performance (完全修复版)""" # 1. 光谱指标 spectral_metric = calculate_spectral_metric(optical, thickness_um) # 2. 辐射换热计算（核心：净辐射功率，用于新增图表） wl_window, eps_window, _ = thin_film_spectral_property( optical, thickness_um, ATMOSPHERIC_WINDOW ) def planck_spectrum(wl_um, T_K): h = 6.626e-34 # Planck constant (J·s) c = 3e8 # Speed of light (m/s) k = 1.38e-23 # Boltzmann constant (J/K) wl_m = wl_um * 1e-6 return (2 * h * c ** 2) / (wl_m ** 5 * (np.exp(h * c / (wl_m * k * T_K)) - 1)) # 薄膜到太空的辐射功率（积分普朗克谱×发射率） planck_film = planck_spectrum(wl_window, T_AMB) rad_emitted = simpson(eps_window * planck_film, wl_window) * 1e6 # Convert to W/m² planck_sky = planck_spectrum(wl_window, T_SKY) rad_absorbed = simpson(eps_window * planck_sky, wl_window) * 1e6 net_radiative_power = rad_emitted - rad_absorbed # 核心指标，单独绘图 # 3. 太阳吸收热流（固定值，与温度无关） solar_gain = spectral_metric.avg_alpha_solar * SOLAR_IRR # 4. 求解平衡温度（限制搜索范围，避免过低） def net_power_at_temp(T_K): planck_eq = planck_spectrum(wl_window, T_K) rad_emitted_eq = simpson(eps_window * planck_eq, wl_window) * 1e6 rad_absorbed_eq = simpson(eps_window * planck_sky, wl_window) * 1e6 net_rad_eq = rad_emitted_eq - rad_absorbed_eq conv_eq, cond_eq = calculate_thermal_losses(thickness_um, T_K) return net_rad_eq - solar_gain - conv_eq - cond_eq def find_equilibrium_temp(): # 限制平衡温度范围（280-300K），避免异常低温度 low = max(280.0, T_AMB - MAX_TEMPERATURE_DIFF) high = min(300.0, T_AMB + 5.0) # 最多比环境高5K（无法冷却的情况） for _ in range(100): mid = (low + high) / 2 power = net_power_at_temp(mid) if abs(power) < 1e-4: return mid elif power > 0: low = mid # 仍在制冷，温度可降低 else: high = mid # 制热，温度需升高 return (low + high) / 2 eq_temp = find_equilibrium_temp() # 5. 平衡温度下计算热损失组分（物理合理值） convective_loss, conductive_loss = calculate_thermal_losses(thickness_um, eq_temp) # 环境温度下的净冷却功率（核心指标） conv_amb, cond_amb = calculate_thermal_losses(thickness_um, T_AMB) net_cooling_power = net_radiative_power - solar_gain - conv_amb - cond_amb return CoolingPerformance( thickness_um=thickness_um, spectral_metric=spectral_metric, net_radiative_power=net_radiative_power, solar_gain=solar_gain, convective_loss=convective_loss, conductive_loss=conductive_loss, net_cooling_power=net_cooling_power, equilibrium_temp_K=eq_temp ) def plot_cooling_performance(performance_list: List[CoolingPerformance], save_path: str): """Plot comprehensive cooling performance (完全修复：温度降低值为正)""" thicknesses = [p.thickness_um for p in performance_list] net_cooling = [p.net_cooling_power for p in performance_list] temp_drop = [p.temp_drop_C for p in performance_list] fig, ax1 = plt.subplots(figsize=(10, 6)) # Left axis: Net Cooling Power color1 = '#2E86AB' ax1.set_xlabel('PDMS Film Thickness (μm)', fontsize=12) ax1.set_ylabel('Net Cooling Power (W/m²)', fontsize=12, color=color1) line1 = ax1.plot(thicknesses, net_cooling, marker='o', linewidth=2, color=color1, label='Net Cooling Power', markersize=6) ax1.tick_params(axis='y', labelcolor=color1) ax1.axhline(y=0, color='black', linestyle='--', alpha=0.5) ax1.grid(True, alpha=0.3) # 调整y轴范围，使正功率更清晰 ax1.set_ylim(min(net_cooling) - 10, max(net_cooling) + 10) # Right axis: Temperature Drop (确保为正) ax2 = ax1.twinx() color2 = '#A23B72' ax2.set_ylabel('Temperature Drop (°C)', fontsize=12, color=color2) line2 = ax2.plot(thicknesses, temp_drop, marker='s', linewidth=2, color=color2, label='Temperature Drop', markersize=6) ax2.tick_params(axis='y', labelcolor=color2) ax2.set_ylim(0, max(temp_drop) + 2) # y轴从0开始，符合降温直觉 # Combine legends lines = line1 + line2 labels = [l.get_label() for l in lines] ax1.legend(lines, labels, loc='upper right', fontsize=10) plt.title('Comprehensive Cooling Performance of PDMS Films vs Thickness', fontsize=14, fontweight='bold') plt.tight_layout() plt.savefig(save_path, bbox_inches='tight') print(f"Cooling performance plot saved to: {save_path}") # ------------------------------ # 8. Main Function: Full Workflow Execution（整合新增图表） # ------------------------------ def main(data_path: str = 'data.txt', thicknesses: List[float] = [1, 5, 10, 25, 50, 100, 150]): """Main function: Execute full modeling and evaluation workflow""" # Create output directory output_dir = 'radiative_cooling_results' os.makedirs(output_dir, exist_ok=True) try: # 1. Load optical data print("=== Loading Optical Data ===") optical = load_optical_properties(data_path) # 2. Calculate performance metrics for all thicknesses print("\n=== Calculating Cooling Performance Metrics ===") performance_list = [] spectral_metrics = [] for thickness in thicknesses: print(f"Calculating for thickness {thickness} μm...") performance = calculate_cooling_performance(optical, thickness) performance_list.append(performance) spectral_metrics.append(performance.spectral_metric) # 3. Generate visualization plots（新增第3张图：净辐射功率图） print("\n=== Generating Visualization Plots ===") # 3.1 光谱匹配度图 plot_spectral_metrics(spectral_metrics, os.path.join(output_dir, '1_spectral_metrics.png')) # 3.2 太阳吸收光谱图 plot_solar_absorption_spectra(optical, thicknesses, os.path.join(output_dir, '2_solar_absorption_spectra.png')) # 3.3 新增：净辐射功率图 plot_net_radiative_power(performance_list, os.path.join(output_dir, '3_net_radiative_power.png')) # 3.4 热损失组分图 plot_thermal_losses(performance_list, os.path.join(output_dir, '4_thermal_loss_components.png')) # 3.5 综合冷却性能图 plot_cooling_performance(performance_list, os.path.join(output_dir, '5_comprehensive_cooling_performance.png')) # 4. Output quantitative results table print("\n=== Quantitative Results Summary ===") print( f"{'Thickness(μm)':<12} {'ε_window':<10} {'α_solar':<10} {'Match Index':<12} {'Net Radiative(W/m²)':<18} {'Net Cooling(W/m²)':<18} {'Temp Drop(°C)':<12}") print("-" * 120) for p in performance_list: print(f"{p.thickness_um:<12.1f} {p.spectral_metric.avg_eps_window:<10.4f} " f"{p.spectral_metric.avg_alpha_solar:<10.4f} {p.spectral_metric.spectral_match_index:<12.2f} " f"{p.net_radiative_power:<18.2f} {p.net_cooling_power:<18.2f} {p.temp_drop_C:<12.2f}") # 5. Save results to CSV（新增净辐射功率列） csv_path = os.path.join(output_dir, 'cooling_performance_results.csv') with open(csv_path, 'w', encoding='utf-8') as f: f.write( "Thickness(μm),Avg_Emissivity_8-13μm,Avg_Absorptivity_0.3-2.5μm,Spectral_Match_Index,Net_Radiative_Power(W/m²),Solar_Gain(W/m²),Convective_Loss(W/m²),Conductive_Loss(W/m²),Net_Cooling_Power(W/m²),Equilibrium_Temp(°C),Temperature_Drop(°C)\n") for p in performance_list: f.write(f"{p.thickness_um:.1f}," f"{p.spectral_metric.avg_eps_window:.4f}," f"{p.spectral_metric.avg_alpha_solar:.4f}," f"{p.spectral_metric.spectral_match_index:.2f}," f"{p.net_radiative_power:.2f}," f"{p.solar_gain:.2f}," f"{p.convective_loss:.2f}," f"{p.conductive_loss:.2f}," f"{p.net_cooling_power:.2f}," f"{p.equilibrium_temp_C:.2f}," f"{p.temp_drop_C:.2f}\n") print(f"\nResults saved to CSV file: {csv_path}") # 6. Output technical recommendations print("\n=== Technical Recommendations ===") # Find optimal thickness based on spectral matching index and net cooling power best_match_idx = np.argmax([p.spectral_metric.spectral_match_index for p in performance_list]) best_cooling_idx = np.argmax([p.net_cooling_power for p in performance_list]) best_radiative_idx = np.argmax([p.net_radiative_power for p in performance_list]) # 新增：净辐射功率最优厚度 best_match_thickness = performance_list[best_match_idx].thickness_um best_cooling_thickness = performance_list[best_cooling_idx].thickness_um best_radiative_thickness = performance_list[best_radiative_idx].thickness_um print( f"1. Optimal thickness for spectral matching: {best_match_thickness} μm (Spectral Matching Index: {performance_list[best_match_idx].spectral_metric.spectral_match_index:.2f})") print( f"2. Optimal thickness for net radiative power: {best_radiative_thickness} μm (Net Radiative Power: {performance_list[best_radiative_idx].net_radiative_power:.2f} W/m²)") print( f"3. Optimal thickness for cooling power: {best_cooling_thickness} μm (Net Cooling Power: {performance_list[best_cooling_idx].net_cooling_power:.2f} W/m²)") print( f"4. Recommended practical thickness: 25-50 μm (Balancing performance and cost, temperature drop: {performance_list[3].temp_drop_C:.1f}-{performance_list[4].temp_drop_C:.1f} °C)") print( "5. Material optimization direction: Reduce absorptivity in solar band (0.3-2.5μm), improve emissivity uniformity in atmospheric window (8-13μm)") print( "6. Engineering application suggestions: Optimize film-substrate interface design (e.g., add insulation layer) to reduce conductive loss; combine with shading structures to minimize solar absorption") except Exception as e: print(f"Error during execution: {e}") import traceback traceback.print_exc() # ------------------------------ # Execution Entry (Update your data file path) # ------------------------------ if __name__ == "__main__": # Configure data path and study thicknesses DATA_PATH = "/Users/spasolreisa/IdeaProjects/asiaMath/data.txt" # Replace with your data.txt path STUDY_THICKNESSES = [1, 5, 10, 25, 50, 100, 150,200 ,250,300] # Adjustable thickness range main(data_path=DATA_PATH, thicknesses=STUDY_THICKNESSES)