231 lines
9.8 KiB
Python
231 lines
9.8 KiB
Python
|
import pandas as pd
|
|||
|
|
import math
|
|||
|
|
from scipy.optimize import fsolve
|
|||
|
|
import requests
|
|||
|
|
import json
|
|||
|
|
import os
|
|||
|
|
|
|||
|
|
def wind_farm_analysis(device_name, area_km2, electricity_price, file_path, latitude, longitude,
|
|||
|
|
lateral_spacing_factor=5, longitudinal_spacing_factor=10, q=0.02, altitude=11,
|
|||
|
|
hub_height=100, Cp=0.45, eta=0.8, cost_per_kw=5):
|
|||
|
|
"""
|
|||
|
|
封装函数：分析风电场的风机数量及各项经济和技术指标，使用 NASA POWER API 获取风速和温度数据
|
|||
|
|
|
|||
|
|
参数：
|
|||
|
|
device_name (str): 设备名称
|
|||
|
|
area_km2 (float): 风电场面积（平方公里）
|
|||
|
|
electricity_price (float): 电价（元/kWh）
|
|||
|
|
file_path (str): 风机参数 Excel 文件路径
|
|||
|
|
latitude (float): 目标地点的纬度（度）
|
|||
|
|
longitude (float): 目标地点的经度（度）
|
|||
|
|
lateral_spacing_factor (float): 横向间距因子（默认为 5D）
|
|||
|
|
longitudinal_spacing_factor (float): 纵向间距因子（默认为 10D）
|
|||
|
|
q (float): 运维成本占初始投资成本的比例（默认 0.02 表示 2%）
|
|||
|
|
altitude (float): 海拔高度（m），默认 11m
|
|||
|
|
hub_height (float): 轮毂高度（m），默认 100m
|
|||
|
|
Cp (float): 风能利用系数，默认 0.45
|
|||
|
|
eta (float): 总系统效率，默认 0.8
|
|||
|
|
cost_per_mw (float): 每 MW 投资成本（万元/MW），默认 5000 万元/MW
|
|||
|
|
|
|||
|
|
返回：
|
|||
|
|
dict: 包含风电场分析结果的字典
|
|||
|
|
"""
|
|||
|
|
def estimate_wind_turbine_count(area_km2, blade_diameter):
|
|||
|
|
area_m2 = area_km2 * 1_000_000
|
|||
|
|
lateral_spacing = lateral_spacing_factor * blade_diameter
|
|||
|
|
longitudinal_spacing = longitudinal_spacing_factor * blade_diameter
|
|||
|
|
turbine_area = lateral_spacing * longitudinal_spacing
|
|||
|
|
turbine_count = int(area_m2 / turbine_area)
|
|||
|
|
print(f"单台风机占地面积: {turbine_area:,} 平方米 "
|
|||
|
|
f"(横向间距: {lateral_spacing} 米, 纵向间距: {longitudinal_spacing} 米)")
|
|||
|
|
print(f"估算风机数量: {turbine_count} 台")
|
|||
|
|
return turbine_count
|
|||
|
|
|
|||
|
|
def get_wind_turbine_specs(device_name, file_path):
|
|||
|
|
try:
|
|||
|
|
df = pd.read_excel(file_path)
|
|||
|
|
match = df[df.iloc[:, 0] == device_name]
|
|||
|
|
if not match.empty:
|
|||
|
|
rated_power = match.iloc[0, 1]
|
|||
|
|
swept_area = match.iloc[0, 7] # 扫风面积
|
|||
|
|
blade_diameter = match.iloc[0, 6] # 叶片直径
|
|||
|
|
print(f"找到设备 '{device_name}'，额定功率: {rated_power} KW, "
|
|||
|
|
f"扫风面积: {swept_area} m², 叶片直径: {blade_diameter} 米")
|
|||
|
|
return rated_power, swept_area, blade_diameter
|
|||
|
|
else:
|
|||
|
|
raise ValueError(f"未找到设备名称: {device_name}")
|
|||
|
|
except FileNotFoundError:
|
|||
|
|
raise FileNotFoundError(f"文件未找到: {file_path}")
|
|||
|
|
except Exception as e:
|
|||
|
|
raise Exception(f"发生错误: {str(e)}")
|
|||
|
|
|
|||
|
|
def fetch_nasa_data(latitude, longitude, start_year="2023", end_year="2023"):
|
|||
|
|
"""
|
|||
|
|
从 NASA POWER API 获取 12 个月的平均气温和风速数据，不保存缓存
|
|||
|
|
"""
|
|||
|
|
try:
|
|||
|
|
url = (f"https://power.larc.nasa.gov/api/temporal/monthly/point?"
|
|||
|
|
f"parameters=T2M,WS10M&community=RE&longitude={longitude}&latitude={latitude}&"
|
|||
|
|
f"start={start_year}&end={end_year}&format=JSON")
|
|||
|
|
response = requests.get(url)
|
|||
|
|
response.raise_for_status() # 检查请求是否成功
|
|||
|
|
data = response.json()
|
|||
|
|
|
|||
|
|
# 提取气温和风速数据
|
|||
|
|
year = start_year
|
|||
|
|
temperatures = [data["properties"]["parameter"]["T2M"][f"{year}{str(i).zfill(2)}"] for i in range(1, 13)]
|
|||
|
|
wind_speeds = [data["properties"]["parameter"]["WS10M"][f"{year}{str(i).zfill(2)}"] for i in range(1, 13)]
|
|||
|
|
|
|||
|
|
# 检查数据完整性
|
|||
|
|
if len(temperatures) != 12 or len(wind_speeds) != 12:
|
|||
|
|
raise ValueError("NASA 数据不完整，未包含 12 个月的数据")
|
|||
|
|
|
|||
|
|
return temperatures, wind_speeds
|
|||
|
|
except requests.exceptions.HTTPError as http_err:
|
|||
|
|
raise Exception(f"HTTP 错误: {str(http_err)}\n响应内容: {response.text}")
|
|||
|
|
except Exception as e:
|
|||
|
|
raise Exception(f"从 NASA POWER API 获取数据时出错: {str(e)}")
|
|||
|
|
|
|||
|
|
def adjust_wind_speed(v_10m, h_ref=10, h_hub=100, alpha=0.143):
|
|||
|
|
"""根据风切变公式调整风速：v_hub = v_ref * (h_hub/h_ref)^alpha"""
|
|||
|
|
return v_10m * (h_hub / h_ref) ** alpha
|
|||
|
|
|
|||
|
|
def air_density(altitude, hub_height, T0):
|
|||
|
|
z = altitude + hub_height
|
|||
|
|
LR = 0.0065
|
|||
|
|
T = T0 - LR * z + 273.15
|
|||
|
|
return (353.05 / T) * math.exp(-0.034 * (z / T))
|
|||
|
|
|
|||
|
|
def wind_power_density(densities, wind_speeds):
|
|||
|
|
sum_rho_v3 = sum(rho * (v ** 3) for rho, v in zip(densities, wind_speeds))
|
|||
|
|
return (1 / (2 * 12)) * sum_rho_v3
|
|||
|
|
|
|||
|
|
def estimated_wind_power(num_turbines, rated_power):
|
|||
|
|
if not isinstance(num_turbines, int) or num_turbines < 0:
|
|||
|
|
raise ValueError("风机数量必须为非负整数")
|
|||
|
|
return rated_power * num_turbines
|
|||
|
|
|
|||
|
|
def calculate_power_output(S, w, Cp, eta,num_turbines):
|
|||
|
|
return w * S * Cp * 8760 * eta *num_turbines
|
|||
|
|
|
|||
|
|
def calculate_equivalent_hours(P, P_r):
|
|||
|
|
if P_r == 0:
|
|||
|
|
raise ValueError("额定功率不能为 0")
|
|||
|
|
#传入的P(发电量)wh,P_r(额定功率)Kw
|
|||
|
|
return (P / 1000) / P_r
|
|||
|
|
|
|||
|
|
def calculate_environmental_benefits(E_p_million_kwh):
|
|||
|
|
if E_p_million_kwh < 0:
|
|||
|
|
raise ValueError("年发电量需≥0")
|
|||
|
|
return {
|
|||
|
|
"coal_reduction": E_p_million_kwh * 0.404 * 10,
|
|||
|
|
"CO2_reduction": E_p_million_kwh * 0.977 * 10,
|
|||
|
|
"SO2_reduction": E_p_million_kwh * 0.03 * 10,
|
|||
|
|
"NOX_reduction": E_p_million_kwh * 0.015 * 10
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
def calculate_reference_yield(E_p, electricity_price, IC, q, n=20):
|
|||
|
|
def npv_equation(irr, p, w, ic, q_val, n=n):
|
|||
|
|
term1 = (1 + irr) ** (-1)
|
|||
|
|
term2 = irr * (1 + irr) ** (-1) if irr != 0 else float('inf')
|
|||
|
|
pv_revenue = p * w * (term1 / term2) * (1 - (1 + irr) ** (-n))
|
|||
|
|
pv_salvage = q_val * ic * (term1 / term2) * (1 - (1 + irr) ** (-n))
|
|||
|
|
return pv_revenue - ic + pv_salvage
|
|||
|
|
|
|||
|
|
irr_guess = 0.1
|
|||
|
|
irr = float(fsolve(npv_equation, irr_guess, args=(E_p, electricity_price, IC, q))[0])
|
|||
|
|
if not 0 <= irr <= 1:
|
|||
|
|
raise ValueError(f"IRR计算结果{irr:.4f}不合理")
|
|||
|
|
return irr * 100
|
|||
|
|
|
|||
|
|
# 获取设备信息
|
|||
|
|
rated_power, swept_area, blade_diameter = get_wind_turbine_specs(device_name, file_path)
|
|||
|
|
|
|||
|
|
# 估算风机数量
|
|||
|
|
num_turbines = estimate_wind_turbine_count(area_km2, blade_diameter)
|
|||
|
|
|
|||
|
|
# 从 NASA POWER API 获取气温和风速数据
|
|||
|
|
monthly_temps, wind_speeds = fetch_nasa_data(latitude, longitude, start_year="2023", end_year="2023")
|
|||
|
|
|
|||
|
|
# 调整风速到轮毂高度
|
|||
|
|
wind_speeds = [adjust_wind_speed(v) for v in wind_speeds]
|
|||
|
|
|
|||
|
|
# 计算空气密度
|
|||
|
|
densities = [air_density(altitude, hub_height, T0) for T0 in monthly_temps]
|
|||
|
|
avg_density = sum(densities) / len(densities)
|
|||
|
|
|
|||
|
|
# 计算风功率密度 w/m2
|
|||
|
|
wpd = wind_power_density(densities, wind_speeds)
|
|||
|
|
|
|||
|
|
# 计算装机容量 KW
|
|||
|
|
total_power = estimated_wind_power(num_turbines, rated_power)
|
|||
|
|
|
|||
|
|
# 计算初始投资成本
|
|||
|
|
IC = total_power * cost_per_kw * 1000
|
|||
|
|
|
|||
|
|
# 计算年发电量 Wh
|
|||
|
|
P_test = calculate_power_output(swept_area, wpd, Cp, eta,num_turbines)
|
|||
|
|
|
|||
|
|
# 计算等效小时数 年发电量(Wh)/额定功率(KW)
|
|||
|
|
h = calculate_equivalent_hours(P_test, rated_power)
|
|||
|
|
# 计算环境收益（转换为万 kWh）
|
|||
|
|
E_p_million_kwh = P_test / 10000000 # 转换为 万 kWh
|
|||
|
|
env_benefits = calculate_environmental_benefits(E_p_million_kwh)
|
|||
|
|
# 计算 IRR
|
|||
|
|
P_test_IRR = P_test/1000
|
|||
|
|
irr = calculate_reference_yield(P_test_IRR, electricity_price, IC, q)
|
|||
|
|
|
|||
|
|
# 返回结果
|
|||
|
|
return {
|
|||
|
|
"device": device_name,
|
|||
|
|
"rated_power": rated_power,
|
|||
|
|
"swept_area": swept_area,
|
|||
|
|
"blade_diameter": blade_diameter,
|
|||
|
|
"num_turbines": num_turbines,
|
|||
|
|
"avg_density": avg_density,
|
|||
|
|
"wpd": wpd,
|
|||
|
|
"total_power": total_power /1000, #变为了MW
|
|||
|
|
"annual_power_output": P_test/10000000 , # 万 kWh
|
|||
|
|
"equivalent_hours": h,
|
|||
|
|
"coal_reduction": env_benefits["coal_reduction"],
|
|||
|
|
"CO2_reduction": env_benefits["CO2_reduction"],
|
|||
|
|
"SO2_reduction": env_benefits["SO2_reduction"],
|
|||
|
|
"NOX_reduction": env_benefits["NOX_reduction"],
|
|||
|
|
"IRR": irr
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
# 主程序
|
|||
|
|
if __name__ == "__main__":
|
|||
|
|
file_path = r"/home/zhaojh/workspace/GreenTransPowerCalculate/wind/wind_product.xlsx"
|
|||
|
|
|
|||
|
|
device_name = 'GW165-4.0'
|
|||
|
|
area_km2 = 10
|
|||
|
|
electricity_price = 0.6
|
|||
|
|
latitude = 39
|
|||
|
|
longitude = 116
|
|||
|
|
|
|||
|
|
result = wind_farm_analysis(
|
|||
|
|
device_name=device_name,
|
|||
|
|
area_km2=area_km2,
|
|||
|
|
electricity_price=electricity_price,
|
|||
|
|
file_path=file_path,
|
|||
|
|
latitude=latitude,
|
|||
|
|
longitude=longitude
|
|||
|
|
)
|
|||
|
|
|
|||
|
|
print(f"\n设备: {result['device']}")
|
|||
|
|
print(f"额定功率: {result['rated_power']:.2f} KW")
|
|||
|
|
print(f"扫风面积: {result['swept_area']:.2f} m^2")
|
|||
|
|
print(f"叶片直径: {result['blade_diameter']:.2f} m")
|
|||
|
|
print(f"风机数量: {result['num_turbines']} 台")
|
|||
|
|
print(f"平均空气密度: {result['avg_density']:.3f} kg/m^3")
|
|||
|
|
print(f"风功率密度: {result['wpd']:.2f} W/m^2")
|
|||
|
|
print(f"项目装机容量: {result['total_power']:.2f} MW")
|
|||
|
|
print(f"年发电量: {result['annual_power_output']:.3f} 万 kWh")
|
|||
|
|
print(f"等效小时数: {result['equivalent_hours']:.2f} 小时")
|
|||
|
|
print(f"CO₂减排量：{result['CO2_reduction']:,.0f} kg")
|
|||
|
|
print(f"SO₂减排量：{result['SO2_reduction']:,.0f} kg")
|
|||
|
|
print(f"NOx减排量：{result['NOX_reduction']:,.0f} kg")
|
|||
|
|
print(f"内部收益率 IRR: {result['IRR']:.2f}%")
|
|||
|
|
|