diff --git a/org/other/outputs/question1_emissivity.png b/org/other/outputs/question1_emissivity.png new file mode 100644 index 0000000..ceedf5b Binary files /dev/null and b/org/other/outputs/question1_emissivity.png differ diff --git a/org/other/outputs/question1_emissivity_summary.csv b/org/other/outputs/question1_emissivity_summary.csv new file mode 100644 index 0000000..a29527a --- /dev/null +++ b/org/other/outputs/question1_emissivity_summary.csv @@ -0,0 +1,7 @@ +thickness_um,avg_emissivity_8_13um +1,0.3238 +5,0.7365 +10,0.8286 +25,0.8987 +50,0.9282 +100,0.9420 \ No newline at end of file diff --git a/org/other/outputs/question2_cooling_results.png b/org/other/outputs/question2_cooling_results.png new file mode 100644 index 0000000..47971bb Binary files /dev/null and b/org/other/outputs/question2_cooling_results.png differ diff --git a/org/other/outputs/question2_cooling_summary.csv b/org/other/outputs/question2_cooling_summary.csv new file mode 100644 index 0000000..24f3708 --- /dev/null +++ b/org/other/outputs/question2_cooling_summary.csv @@ -0,0 +1,7 @@ +thickness_um,eps_8_13,alpha_solar,net_power_amb_Wm2,eq_temp_C,delta_T_C +1,0.3238,0.051,-10.03,25.19,-1.66 +5,0.7365,0.055,32.08,36.68,9.83 +10,0.8286,0.060,37.78,41.76,14.91 +25,0.8987,0.075,32.05,42.47,15.62 +50,0.9282,0.100,12.81,32.83,5.98 +100,0.9420,0.150,-30.65,14.91,-11.94 diff --git a/org/other/outputs/question3_multilayer_summary.csv b/org/other/outputs/question3_multilayer_summary.csv new file mode 100644 index 0000000..aac0210 --- /dev/null +++ b/org/other/outputs/question3_multilayer_summary.csv @@ -0,0 +1,16 @@ +rank,score,epsilon,alpha,layers +1,0.8599,0.9182,0.1944,PDMS@12.857um;HfO2@0.411um;SiO2@1.089um +2,0.8585,0.9185,0.1999,PDMS@14.393um;SiO2@0.470um;HfO2@0.302um +3,0.8506,0.9325,0.2729,PDMS@20.398um;Al2O3@0.510um;SiO2@1.906um +4,0.8499,0.9223,0.2411,PDMS@17.833um;HfO2@0.316um;SiO2@1.127um +5,0.8492,0.9183,0.2304,PDMS@13.284um;HfO2@0.081um;Al2O3@0.371um;Si3N4@0.723um +6,0.8482,0.9355,0.2911,PDMS@21.867um;SiO2@0.309um;Si3N4@0.707um;TiO2@0.057um +7,0.8475,0.9223,0.2495,PDMS@17.611um;SiO2@1.097um;HfO2@0.613um +8,0.8470,0.9269,0.2663,PDMS@19.765um;SiO2@1.968um;Al2O3@0.609um +9,0.8468,0.9221,0.2510,PDMS@15.789um;Si3N4@1.199um;SiO2@1.265um +10,0.8464,0.9264,0.2667,PDMS@18.329um;Si3N4@0.322um;Al2O3@0.532um +11,0.8463,0.9145,0.2273,PDMS@15.894um;Al2O3@1.159um;SiO2@0.497um +12,0.8461,0.9240,0.2595,PDMS@18.179um;Al2O3@0.324um;HfO2@0.734um +13,0.8459,0.9053,0.1980,PDMS@12.520um;HfO2@0.132um;Si3N4@0.431um;SiO2@0.435um +14,0.8446,0.9091,0.2150,PDMS@12.045um;Al2O3@0.692um;HfO2@0.660um;SiO2@0.207um +15,0.8445,0.9401,0.3186,PDMS@26.732um;SiO2@0.202um;HfO2@0.346um diff --git a/org/other/outputs/question3_pareto.png b/org/other/outputs/question3_pareto.png new file mode 100644 index 0000000..dd95a08 Binary files /dev/null and b/org/other/outputs/question3_pareto.png differ diff --git a/org/other/q2_1.py b/org/other/q2_1.py new file mode 100644 index 0000000..541c857 --- /dev/null +++ b/org/other/q2_1.py @@ -0,0 +1,667 @@ +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) \ No newline at end of file diff --git a/org/other/radiative_cooling_results/1_spectral_metrics.png b/org/other/radiative_cooling_results/1_spectral_metrics.png new file mode 100644 index 0000000..6058ee1 Binary files /dev/null and b/org/other/radiative_cooling_results/1_spectral_metrics.png differ diff --git a/org/other/radiative_cooling_results/2_solar_absorption_spectra.png b/org/other/radiative_cooling_results/2_solar_absorption_spectra.png new file mode 100644 index 0000000..070f3e5 Binary files /dev/null and b/org/other/radiative_cooling_results/2_solar_absorption_spectra.png differ diff --git a/org/other/radiative_cooling_results/3_net_radiative_power.png b/org/other/radiative_cooling_results/3_net_radiative_power.png new file mode 100644 index 0000000..4430b07 Binary files /dev/null and b/org/other/radiative_cooling_results/3_net_radiative_power.png differ diff --git a/org/other/radiative_cooling_results/4_thermal_loss_components.png b/org/other/radiative_cooling_results/4_thermal_loss_components.png new file mode 100644 index 0000000..f5a38c1 Binary files /dev/null and b/org/other/radiative_cooling_results/4_thermal_loss_components.png differ diff --git a/org/other/radiative_cooling_results/5_comprehensive_cooling_performance.png b/org/other/radiative_cooling_results/5_comprehensive_cooling_performance.png new file mode 100644 index 0000000..cd3012c Binary files /dev/null and b/org/other/radiative_cooling_results/5_comprehensive_cooling_performance.png differ diff --git a/org/other/radiative_cooling_results/cooling_performance_results.csv b/org/other/radiative_cooling_results/cooling_performance_results.csv new file mode 100644 index 0000000..2ff3047 --- /dev/null +++ b/org/other/radiative_cooling_results/cooling_performance_results.csv @@ -0,0 +1,11 @@ +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) +1.0,0.2198,0.0001,1832.21,4213386159920.65,0.10,0.00,0.00,4213386159920.55,26.85,0.00 +5.0,0.5527,0.0006,938.24,10471946524462.82,0.50,0.00,0.00,10471946524462.32,26.85,0.00 +10.0,0.6845,0.0012,573.18,13000419684891.89,1.01,0.00,0.00,13000419684890.88,26.85,0.00 +25.0,0.8256,0.0030,278.83,15732438394151.59,2.52,0.00,0.00,15732438394149.08,26.85,0.00 +50.0,0.8964,0.0059,151.79,17117210018907.72,5.02,0.00,0.00,17117210018902.70,26.85,0.00 +100.0,0.9338,0.0117,79.78,17842980949637.75,9.95,0.00,0.00,17842980949627.80,26.85,0.00 +150.0,0.9423,0.0174,54.22,18006953306435.60,14.77,0.00,0.00,18006953306420.83,26.85,0.00 +200.0,0.9446,0.0229,41.22,18051728737920.30,19.48,0.00,0.00,18051728737900.82,26.85,0.00 +250.0,0.9453,0.0286,33.09,18064943852359.92,24.28,0.00,0.00,18064943852335.64,26.85,0.00 +300.0,0.9455,0.0339,27.91,18069013293873.77,28.79,0.00,0.00,18069013293844.97,26.85,0.00 diff --git a/org/use/PDMS_emissivity_spectrum.png b/org/use/PDMS_emissivity_spectrum.png new file mode 100644 index 0000000..2f7d7c0 Binary files /dev/null and b/org/use/PDMS_emissivity_spectrum.png differ diff --git a/org/use/outputs/advanced_emissivity_spectrum.png b/org/use/outputs/advanced_emissivity_spectrum.png new file mode 100644 index 0000000..381eb53 Binary files /dev/null and b/org/use/outputs/advanced_emissivity_spectrum.png differ diff --git a/org/use/outputs/question2_cooling_results_advanced.png b/org/use/outputs/question2_cooling_results_advanced.png new file mode 100644 index 0000000..d90b090 Binary files /dev/null and b/org/use/outputs/question2_cooling_results_advanced.png differ diff --git a/org/use/outputs/question2_cooling_summary_advanced.csv b/org/use/outputs/question2_cooling_summary_advanced.csv new file mode 100644 index 0000000..a0a5d5f --- /dev/null +++ b/org/use/outputs/question2_cooling_summary_advanced.csv @@ -0,0 +1,7 @@ +thickness_um,eps_8_13,alpha_solar,net_power_amb_Wm2,eq_temp_C,delta_T_C +1,0.2197,0.000,24.23,30.50,3.65 +5,0.5527,0.001,60.67,40.72,13.87 +10,0.6846,0.001,74.71,49.42,22.57 +25,0.8257,0.003,88.68,56.85,30.00 +50,0.8965,0.006,93.75,56.85,30.00 +100,0.9338,0.012,92.50,56.85,30.00 diff --git a/org/use/outputs/question3_multilayer_summary.csv b/org/use/outputs/question3_multilayer_summary.csv new file mode 100644 index 0000000..09ecbbe --- /dev/null +++ b/org/use/outputs/question3_multilayer_summary.csv @@ -0,0 +1,16 @@ +rank,score,epsilon,alpha,layers +1,0.9268,0.9334,0.0219,PDMS@49.640um;TiO2@0.118um;SiO2@0.411um +2,0.9245,0.9327,0.0275,PDMS@48.004um;Si3N4@0.391um;SiO2@0.066um +3,0.9222,0.9321,0.0332,PDMS@47.171um;SiO2@0.265um;HfO2@0.473um +4,0.9218,0.9324,0.0353,PDMS@48.703um;HfO2@0.469um;SiO2@1.280um +5,0.9200,0.9334,0.0448,PDMS@49.757um;TiO2@0.260um;SiO2@0.636um +6,0.9190,0.9307,0.0390,PDMS@46.213um;SiO2@1.752um;HfO2@0.496um +7,0.9188,0.9266,0.0261,PDMS@40.876um;HfO2@0.378um;SiO2@0.157um +8,0.9176,0.9241,0.0219,PDMS@39.132um;SiO2@1.429um;Al2O3@0.277um +9,0.9174,0.9268,0.0316,PDMS@41.884um;TiO2@0.146um;SiO2@1.911um +10,0.9172,0.9334,0.0538,PDMS@49.777um;Al2O3@0.683um;Si3N4@0.289um +11,0.9171,0.9258,0.0290,PDMS@41.247um;Al2O3@0.434um;SiO2@1.647um +12,0.9163,0.9293,0.0433,PDMS@43.750um;SiO2@0.225um;Al2O3@0.822um +13,0.9158,0.9264,0.0356,PDMS@41.092um;SiO2@0.313um;Al2O3@0.628um +14,0.9150,0.9289,0.0462,PDMS@44.669um;SiO2@1.676um;Al2O3@0.840um +15,0.9145,0.9331,0.0620,PDMS@48.926um;Si3N4@0.865um;SiO2@1.509um diff --git a/org/use/outputs/question3_pareto.png b/org/use/outputs/question3_pareto.png new file mode 100644 index 0000000..54be2f3 Binary files /dev/null and b/org/use/outputs/question3_pareto.png differ diff --git a/org/use/q2.py b/org/use/q2.py deleted file mode 100644 index 8553eb8..0000000 --- a/org/use/q2.py +++ /dev/null @@ -1,385 +0,0 @@ -import importlib.util -import math -import os -from dataclasses import dataclass -from typing import Dict, List, Tuple - -import numpy as np -from matplotlib import pyplot as plt -from scipy.interpolate import CubicSpline - -plt.rcParams["font.sans-serif"] = ["DejaVu Sans", "Arial", "Helvetica"] -plt.rcParams["axes.unicode_minus"] = False - -SIGMA = 5.670374419e-8 -T_AMB = 300.0 # K (≈27 ℃) -T_SKY = 280.0 # Clear dry sky equivalent radiation temperature -SOLAR_IRR = 900.0 # W/m^2, clear sky noon -H_CONV = 8.0 # W/m^2/K, natural convection + light wind - - -def make_strictly_increasing(wl, n, k): - """确保波长数据严格递增""" - 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] - - is_increasing = np.diff(wl) > 0 - if not all(is_increasing): - valid_indices = np.concatenate([[True], is_increasing]) - wl, n, k = wl[valid_indices], n[valid_indices], k[valid_indices] - - return wl, n, k - - -def read_split_data(file_path): - """读取分块格式的光学数据""" - with open(file_path, 'r', encoding='utf-8') as f: - lines = [line.strip() for line in f if line.strip() and not line.startswith('#')] - - split_idx = None - for i, line in enumerate(lines): - if line == 'wl k': - split_idx = i - break - - n_lines = lines[1:split_idx] - 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])) - - k_lines = lines[split_idx + 1:] - 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])) - - wl_n = np.array(wl_n) - n_list = np.array(n_list) - wl_k = np.array(wl_k) - k_list = np.array(k_list) - - # 检查波长范围是否一致 - if not np.array_equal(wl_n, wl_k): - print(f"警告: n和k的波长范围不完全一致") - print(f"n波长范围: {wl_n.min():.2f} - {wl_n.max():.2f} μm") - print(f"k波长范围: {wl_k.min():.2f} - {wl_k.max():.2f} μm") - - # 找到共同的波长范围 - common_wl = np.intersect1d(wl_n, wl_k) - if len(common_wl) == 0: - raise ValueError("n和k没有共同的波长点！") - - # 插值到共同的波长点 - n_common = np.interp(common_wl, wl_n, n_list) - k_common = np.interp(common_wl, wl_k, k_list) - - wl_n, n_list, wl_k, k_list = common_wl, n_common, common_wl, k_common - - sorted_idx = np.argsort(wl_n) - sorted_wl = wl_n[sorted_idx] - sorted_n = n_list[sorted_idx] - sorted_k = k_list[sorted_idx] - - return sorted_wl, sorted_n, sorted_k - - -def fresnel_reflectance(n1, k1, n2, k2): - """菲涅尔反射率计算""" - m1 = n1 + 1j * k1 - m2 = n2 + 1j * k2 - return np.abs((m1 - m2) / (m1 + m2)) ** 2 - - -def thin_film_emissivity(n_film, k_film, d, wl, n_air=1.0, k_air=0.0): - """精确的薄膜发射率计算""" - R12 = fresnel_reflectance(n_air, k_air, n_film, k_film) - R23 = fresnel_reflectance(n_film, k_film, n_air, k_air) - - delta = 2 * np.pi * n_film * d / wl - alpha = 4 * np.pi * k_film * d / wl - - 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 - - T_total = (1 - R12) * (1 - R23) * np.exp(-alpha) / denominator - - epsilon = 1 - R_total - T_total - return np.clip(epsilon, 0, 1) - - -def solar_absorptance(thickness_um: float, cs_n, cs_k, wl_range) -> float: - """ - 基于精确光学模型计算太阳吸收率 - """ - # 太阳光谱范围 - 确保在数据范围内 - wl_solar_min = max(0.3, wl_range.min()) - wl_solar_max = min(2.5, wl_range.max()) - - if wl_solar_min >= wl_solar_max: - print(f"警告: 数据波长范围 {wl_range.min():.2f}-{wl_range.max():.2f}μm 不覆盖太阳光谱") - return min(0.05 + 0.001 * thickness_um, 0.15) - - wl_solar = np.linspace(wl_solar_min, wl_solar_max, 200) - - try: - n_solar = cs_n(wl_solar) - k_solar = cs_k(wl_solar) - alpha_spectral = thin_film_emissivity(n_solar, k_solar, thickness_um, wl_solar) - avg_alpha = np.mean(alpha_spectral) - return float(np.clip(avg_alpha, 0, 1)) - - except Exception as e: - print(f"太阳吸收率计算失败 {thickness_um}μm: {e}") - return min(0.05 + 0.001 * thickness_um, 0.15) - - -@dataclass -class CoolingResult: - thickness_um: float - eps_window: float - alpha_solar: float - net_power_at_amb: float - eq_temp_K: float - - @property - def eq_temp_C(self) -> float: - return self.eq_temp_K - 273.15 - - @property - def delta_T(self) -> float: - return self.eq_temp_K - T_AMB - - -def net_cooling_power(temp_K: float, emissivity: float, alpha_s: float) -> float: - """计算净冷却功率""" - radiative = emissivity * SIGMA * (temp_K ** 4 - T_SKY ** 4) - solar_gain = alpha_s * SOLAR_IRR - convective = H_CONV * (temp_K - T_AMB) - return -solar_gain - convective + radiative - - -def solve_equilibrium(emissivity: float, alpha_s: float) -> float: - """求解平衡温度""" - low, high = 250.0, 330.0 - f_low = net_cooling_power(low, emissivity, alpha_s) - f_high = net_cooling_power(high, emissivity, alpha_s) - - if f_low * f_high > 0: - return low if abs(f_low) < abs(f_high) else high - - for _ in range(80): - mid = 0.5 * (low + high) - f_mid = net_cooling_power(mid, emissivity, alpha_s) - if abs(f_mid) < 1e-4: - return mid - if f_low * f_mid <= 0: - high, f_high = mid, f_mid - else: - low, f_low = mid, f_mid - return 0.5 * (low + high) - - -def load_optical_data(): - """加载光学数据并创建插值函数""" - data_path = '/Users/spasolreisa/IdeaProjects/asiaMath/data.txt' - print(f"正在加载光学数据: {data_path}") - wl_all, n_all, k_all = read_split_data(data_path) - wl_all, n_all, k_all = make_strictly_increasing(wl_all, n_all, k_all) - - print(f"数据波长范围: {wl_all.min():.2f} - {wl_all.max():.2f} μm") - print(f"数据点数: {len(wl_all)}") - print(f"折射率n范围: {n_all.min():.3f} - {n_all.max():.3f}") - print(f"消光系数k范围: {k_all.min():.3f} - {k_all.max():.3f}") - - cs_n = CubicSpline(wl_all, n_all) - cs_k = CubicSpline(wl_all, k_all) - - return wl_all, cs_n, cs_k - - -def compute_emissivity_for_thickness(thicknesses_um, cs_n, cs_k, wl_range): - """计算各厚度在大气窗口的平均发射率""" - # 大气窗口波长范围 - 确保在数据范围内 - wl_window_min = max(8.0, wl_range.min()) - wl_window_max = min(13.0, wl_range.max()) - - if wl_window_min >= wl_window_max: - print(f"警告: 数据波长范围 {wl_range.min():.2f}-{wl_range.max():.2f}μm 不覆盖大气窗口") - # 返回默认值 - return {}, {d: 0.8 for d in thicknesses_um} - - wl_window = np.linspace(wl_window_min, wl_window_max, 300) - - emissivity_map = {} - window_avg = {} - - for d in thicknesses_um: - try: - n_window = cs_n(wl_window) - k_window = cs_k(wl_window) - epsilon_window = thin_film_emissivity(n_window, k_window, d, wl_window) - avg_epsilon = np.mean(epsilon_window) - - emissivity_map[d] = epsilon_window - window_avg[d] = float(avg_epsilon) - print(f"厚度 {d}μm: 大气窗口发射率 = {avg_epsilon:.4f}") - except Exception as e: - print(f"发射率计算失败 {d}μm: {e}") - emissivity_map[d] = np.full_like(wl_window, 0.8) - window_avg[d] = 0.8 - - return emissivity_map, window_avg - - -def evaluate() -> Tuple[List[CoolingResult], str]: - """主评估函数""" - # 加载光学数据 - try: - wl_all, cs_n, cs_k = load_optical_data() - except Exception as e: - print(f"加载光学数据失败: {e}") - return [], "" - - # 定义研究的厚度 - 使用与原始代码相同的厚度 - THICKNESSES = [1, 5, 10, 25, 50, 100] - - # 计算发射率 - print("\n计算大气窗口发射率...") - emissivity_map, window_avg = compute_emissivity_for_thickness(THICKNESSES, cs_n, cs_k, wl_all) - - results: List[CoolingResult] = [] - for d in THICKNESSES: - eps_window = window_avg[d] - alpha_s = solar_absorptance(d, cs_n, cs_k, wl_all) - net_amb = net_cooling_power(T_AMB, eps_window, alpha_s) - eq_temp = solve_equilibrium(eps_window, alpha_s) - - results.append( - CoolingResult( - thickness_um=d, - eps_window=eps_window, - alpha_solar=alpha_s, - net_power_at_amb=net_amb, - eq_temp_K=eq_temp, - ) - ) - - # 保存结果 - outdir = os.path.join(os.path.dirname(__file__), "outputs") - os.makedirs(outdir, exist_ok=True) - - csv_path = os.path.join(outdir, "question2_cooling_summary_advanced.csv") - with open(csv_path, "w", encoding="utf-8") as f: - f.write("thickness_um,eps_8_13,alpha_solar,net_power_amb_Wm2,eq_temp_C,delta_T_C\n") - for res in results: - f.write( - f"{res.thickness_um},{res.eps_window:.4f},{res.alpha_solar:.3f}," - f"{res.net_power_at_amb:.2f},{res.eq_temp_C:.2f},{res.delta_T:.2f}\n" - ) - - fig_path = os.path.join(outdir, "question2_cooling_results_advanced.png") - plot_results(results, fig_path) - - # 绘制发射率光谱 - plot_emissivity_spectrum(THICKNESSES, cs_n, cs_k, wl_all, outdir) - - return results, csv_path - - -def plot_results(results: List[CoolingResult], fig_path: str) -> None: - """绘制冷却性能结果 - 使用与原始代码相同的格式""" - thickness = [r.thickness_um for r in results] - net_power = [r.net_power_at_amb for r in results] - delta_T = [r.delta_T for r in results] - - fig, ax1 = plt.subplots(figsize=(9, 5)) - - # 使用与原始代码相同的条形图格式 - ax1.bar(thickness, net_power, width=4, alpha=0.6, label="Net Cooling Power @T_amb") - ax1.set_xlabel("PDMS Film Thickness (µm)") - ax1.set_ylabel("Net Cooling Power (W/m²)") - ax1.axhline(0, color="black", linewidth=0.8) - - ax2 = ax1.twinx() - ax2.plot(thickness, delta_T, color="tab:red", marker="o", - linewidth=2, markersize=6, label="Equilibrium Temperature Difference") - ax2.set_ylabel("Equilibrium Temperature Difference (K)") - - lines, labels = ax1.get_legend_handles_labels() - lines2, labels2 = ax2.get_legend_handles_labels() - ax1.legend(lines + lines2, labels + labels2, loc="upper right") - ax1.grid(True, alpha=0.3) - fig.tight_layout() - plt.savefig(fig_path, dpi=300) - plt.close(fig) - print(f"冷却性能图保存至: {fig_path}") - - -def plot_emissivity_spectrum(thicknesses, cs_n, cs_k, wl_range, outdir): - """绘制发射率光谱""" - # 使用原始数据的完整波长范围 - wl_min = wl_range.min() - wl_max = wl_range.max() - wl_fine = np.linspace(wl_min, wl_max, 1000) - - plt.figure(figsize=(12, 7)) - - for d in thicknesses: - try: - n_film = cs_n(wl_fine) - k_film = cs_k(wl_fine) - epsilon = thin_film_emissivity(n_film, k_film, d, wl_fine) - plt.plot(wl_fine, epsilon, linewidth=2, label=f'{d} μm') - except Exception as e: - print(f"绘图失败 {d}μm: {e}") - continue - - # 标注大气窗口（如果数据覆盖） - if wl_min <= 13 and wl_max >= 8: - plt.axvspan(8, 13, alpha=0.15, color='red', label='Atmospheric Window (8-13 μm)') - - plt.xlabel('Wavelength (μm)', fontsize=14) - plt.ylabel('Emissivity ε(λ)', fontsize=14) - plt.title('PDMS Thin Film Spectral Emissivity', fontsize=16, fontweight='bold') - plt.grid(True, alpha=0.3, linestyle='--') - plt.legend(fontsize=10, loc='best') - plt.ylim(0, 1.05) - plt.xlim(wl_min, wl_max) - - fig_path = os.path.join(outdir, "advanced_emissivity_spectrum.png") - plt.tight_layout() - plt.savefig(fig_path, dpi=300, bbox_inches='tight') - plt.close() - print(f"发射率光谱图保存至: {fig_path}") - - -def main(): - """主函数""" - print("开始辐射冷却性能分析...") - results, csv_path = evaluate() - - if results: - print(f"\n结果已保存至: {csv_path}") - print("\n=== 辐射冷却性能汇总 ===") - print("厚度(μm) | 大气窗口发射率 | 太阳吸收率 | 净冷却功率(W/m²) | 平衡温度(°C) | ΔT(K)") - print("-" * 85) - for r in results: - print( - f"{r.thickness_um:>7.1f} | {r.eps_window:>13.3f} | {r.alpha_solar:>10.3f} | " - f"{r.net_power_at_amb:>15.1f} | {r.eq_temp_C:>11.1f} | " - f"{r.delta_T:>5.1f}" - ) - else: - print("分析失败，无结果输出。") - - -if __name__ == "__main__": - main() \ No newline at end of file diff --git a/org/use/q2_1.py b/org/use/q2_1.py new file mode 100644 index 0000000..541c857 --- /dev/null +++ b/org/use/q2_1.py @@ -0,0 +1,667 @@ +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) \ No newline at end of file