1
This commit is contained in:
95
org/q2.py
Normal file
95
org/q2.py
Normal file
@@ -0,0 +1,95 @@
|
||||
import numpy as np
|
||||
from scipy.constants import h, c, k
|
||||
from scipy.integrate import simps
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
from org.q1 import emissivity
|
||||
|
||||
|
||||
# ------------------------------------------------
|
||||
# 1. 引用第一题的 emissivity() 函数
|
||||
# ------------------------------------------------
|
||||
# (假设已运行第一题代码)
|
||||
# emissivity(lambda_um, thickness_um)
|
||||
|
||||
|
||||
# ------------------------------------------------
|
||||
# 2. 黑体辐射谱
|
||||
# ------------------------------------------------
|
||||
def planck_lambda(lambda_um, T):
|
||||
lambda_m = lambda_um * 1e-6
|
||||
return (2*h*c**2)/(lambda_m**5) / (np.exp(h*c/(lambda_m*k*T)) - 1)
|
||||
|
||||
# ------------------------------------------------
|
||||
# 3. 大气透过率(示例,可替换 MODTRAN 数据)
|
||||
# ------------------------------------------------
|
||||
def atmosphere_tau(lambda_um):
|
||||
tau = np.ones_like(lambda_um)
|
||||
window = (lambda_um >= 8) & (lambda_um <= 13)
|
||||
tau[window] = 0.8
|
||||
tau[~window] = 0.1
|
||||
return tau
|
||||
|
||||
# ------------------------------------------------
|
||||
# 4. 太阳辐照谱(简化 AM1.5)
|
||||
# ------------------------------------------------
|
||||
def solar_spectrum(lambda_um):
|
||||
I0 = 1.5e3 # simplified scaling
|
||||
return I0 * np.exp(-(lambda_um - 0.5)**2 / 0.4)
|
||||
|
||||
# ------------------------------------------------
|
||||
# 5. 能量项计算
|
||||
# ------------------------------------------------
|
||||
def P_rad(Ts, lambda_um, eps):
|
||||
return simps(eps * planck_lambda(lambda_um, Ts), lambda_um)
|
||||
|
||||
def P_atm(Ta, lambda_um, eps):
|
||||
return simps(eps * planck_lambda(lambda_um, Ta) * (1 - atmosphere_tau(lambda_um)), lambda_um)
|
||||
|
||||
def P_solar(lambda_um, eps):
|
||||
A = eps # approximate absorption
|
||||
return simps(A * solar_spectrum(lambda_um), lambda_um)
|
||||
|
||||
def P_conv(Ts, Ta, h=5):
|
||||
return h * (Ts - Ta)
|
||||
|
||||
# ------------------------------------------------
|
||||
# 6. 净冷却功率
|
||||
# ------------------------------------------------
|
||||
def P_net(Ts, Ta, lambda_um, eps, h=5):
|
||||
return P_rad(Ts, lambda_um, eps) - P_atm(Ta, lambda_um, eps) - P_solar(lambda_um, eps) - P_conv(Ts, Ta, h)
|
||||
|
||||
# ------------------------------------------------
|
||||
# 7. 求稳态温度(解 Pnet=0)
|
||||
# ------------------------------------------------
|
||||
def solve_temperature(Ta, lambda_um, eps):
|
||||
Ts = Ta
|
||||
for _ in range(1000):
|
||||
f = P_net(Ts, Ta, lambda_um, eps)
|
||||
Ts -= 0.1 * f # simple iteration
|
||||
return Ts
|
||||
|
||||
# ------------------------------------------------
|
||||
# 8. 运行示例
|
||||
# ------------------------------------------------
|
||||
lambda_range = np.linspace(0.35, 20, 2000)
|
||||
d_pdms = 10
|
||||
eps = emissivity(lambda_range, d_pdms)
|
||||
|
||||
Ta = 300 # 环境温度
|
||||
Ts = solve_temperature(Ta, lambda_range, eps)
|
||||
|
||||
print("Steady-state temperature:", Ts)
|
||||
|
||||
# 曲线可视化
|
||||
T_list = np.linspace(250, 330, 200)
|
||||
P_list = [P_net(T, Ta, lambda_range, eps) for T in T_list]
|
||||
|
||||
plt.figure(figsize=(10,5))
|
||||
plt.plot(T_list, P_list)
|
||||
plt.axhline(0, color='r')
|
||||
plt.xlabel("Temperature (K)")
|
||||
plt.ylabel("Net Cooling Power")
|
||||
plt.title("Cooling Power vs Temperature")
|
||||
plt.grid(True)
|
||||
plt.show()
|
||||
Reference in New Issue
Block a user