GreenTransPowerCalculate/Wind/wind_total.py

192 lines
8.5 KiB
Python
Raw Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

import pandas as pd
import math
from scipy.optimize import fsolve
def wind_farm_analysis(device_name, area_km2, electricity_price, file_path, velocity_path, T_path,
lateral_spacing_factor=5, longitudinal_spacing_factor=10, q=0.02, altitude=11,
hub_height=100, Cp=0.45, eta=0.8, cost_per_mw=5000):
"""
封装函数:分析风电场的风机数量及各项经济和技术指标
参数:
device_name (str): 设备名称
area_km2 (float): 风电场面积(平方公里)
electricity_price (float): 电价(元/kWh
file_path (str): 风机参数 Excel 文件路径
velocity_path (str): 风速数据文件路径12 个月平均风速)
T_path (str): 温度数据文件路径12 个月平均温度)
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] / 1000 # kW 转换为 MW
swept_area = match.iloc[0, 7] # 扫风面积
blade_diameter = match.iloc[0, 6] # 叶片直径
print(f"找到设备 '{device_name}',额定功率: {rated_power} MW, "
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 read_monthly_temperatures(file_path):
try:
with open(file_path, 'r', encoding='utf-8') as file:
temperatures = [float(line.strip()) for line in file.readlines()]
if len(temperatures) != 12:
raise ValueError(f"温度文件应包含 12 个月的数据,但实际有 {len(temperatures)}")
return temperatures
except Exception as e:
raise Exception(f"读取温度文件时出错: {str(e)}")
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, file_path):
try:
with open(file_path, 'r', encoding='utf-8') as file:
wind_speeds = [float(line.strip()) for line in file.readlines()]
if len(wind_speeds) != 12:
raise ValueError(f"风速文件应包含 12 个月的数据,但实际有 {len(wind_speeds)}")
sum_rho_v3 = sum(rho * (v ** 3) for rho, v in zip(densities, wind_speeds))
return (1 / (2 * 12)) * sum_rho_v3
except Exception as e:
raise Exception(f"读取风速文件时出错: {str(e)}")
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):
return w * S * Cp * 8760 * eta
def calculate_equivalent_hours(P, P_r):
if P_r == 0:
raise ValueError("额定功率不能为 0")
return (P / 1000) / (P_r * 1000)
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
try:
# 获取设备信息
rated_power, swept_area, blade_diameter = get_wind_turbine_specs(device_name, file_path)
# 估算风机数量
num_turbines = estimate_wind_turbine_count(area_km2, blade_diameter)
# 读取温度数据并计算空气密度
monthly_temps = read_monthly_temperatures(T_path)
densities = [air_density(altitude, hub_height, T0) for T0 in monthly_temps]
avg_density = sum(densities) / len(densities)
# 计算风功率密度
wpd = wind_power_density(densities, velocity_path)
# 计算装机容量
total_power = estimated_wind_power(num_turbines, rated_power)
# 计算初始投资成本
IC = total_power * cost_per_mw * 1000000
# 计算年发电量
P_test = calculate_power_output(swept_area, wpd, Cp, eta) * num_turbines
# 计算等效小时数
h = calculate_equivalent_hours(P_test, rated_power)
# 计算 IRR
irr = calculate_reference_yield(P_test, 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,
"annual_power_output": P_test / 10000000, # 万 kWh
"equivalent_hours": h,
"IRR": irr
}
except Exception as e:
raise Exception(f"风电场分析出错: {str(e)}")
# 主程序
if __name__ == "__main__":
file_path = r".\wind_product.xlsx"
velocity_path = r".\wind_speed.txt"
T_path = r".\temperature.txt"
device_name = input("请输入设备名称: ")
area_km2 = float(input("请输入风电场面积(平方公里): "))
electricity_price = float(input("请输入电价(元/kWh: "))
try:
result = wind_farm_analysis(
device_name=device_name,
area_km2=area_km2,
electricity_price=electricity_price,
file_path=file_path,
velocity_path=velocity_path,
T_path=T_path
)
print(f"\n设备: {result['device']}")
print(f"额定功率: {result['rated_power']:.2f} MW")
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"内部收益率 IRR: {result['IRR']:.2f}%")
except Exception as e:
print(f"错误: {str(e)}")