pv_workflow/pv_assessment/solar_calculator.py

222 lines
10 KiB
Python
Raw Normal View History

import math
import pandas as pd
from datetime import datetime
import hub.helpers.constants as cte
from hub.helpers.monthly_values import MonthlyValues
class SolarCalculator:
def __init__(self, city, tilt_angle, surface_azimuth_angle, standard_meridian=-75,
solar_constant=1366.1, maximum_clearness_index=1, min_cos_zenith=0.065, maximum_zenith_angle=87):
"""
A class to calculate the solar angles and solar irradiance on a tilted surface in the City
:param city: An object from the City class -> City
:param tilt_angle: tilt angle of surface -> float
:param surface_azimuth_angle: The orientation of the surface. 0 is North -> float
:param standard_meridian: A standard meridian is the meridian whose mean solar time is the basis of the time of day
observed in a time zone -> float
:param solar_constant: The amount of energy received by a given area one astronomical unit away from the Sun. It is
constant and must not be changed
:param maximum_clearness_index: This is used to calculate the diffuse fraction of the solar irradiance -> float
:param min_cos_zenith: This is needed to avoid unrealistic values in tilted irradiance calculations -> float
:param maximum_zenith_angle: This is needed to avoid negative values in tilted irradiance calculations -> float
"""
self.city = city
self.location_latitude = city.latitude
self.location_longitude = city.longitude
self.location_latitude_rad = math.radians(self.location_latitude)
self.surface_azimuth_angle = surface_azimuth_angle
self.surface_azimuth_rad = math.radians(surface_azimuth_angle)
self.tilt_angle = tilt_angle
self.tilt_angle_rad = math.radians(tilt_angle)
self.standard_meridian = standard_meridian
self.longitude_correction = (self.location_longitude - standard_meridian) * 4
self.solar_constant = solar_constant
self.maximum_clearness_index = maximum_clearness_index
self.min_cos_zenith = min_cos_zenith
self.maximum_zenith_angle = maximum_zenith_angle
timezone_offset = int(-standard_meridian / 15)
self.timezone = f'Etc/GMT{"+" if timezone_offset < 0 else "-"}{abs(timezone_offset)}'
self.eot = []
self.ast = []
self.hour_angles = []
self.declinations = []
self.solar_altitudes = []
self.solar_azimuths = []
self.zeniths = []
self.incidents = []
self.i_on = []
self.i_oh = []
self.times = pd.date_range(start='2023-01-01', end='2023-12-31 23:00', freq='h', tz=self.timezone)
self.solar_angles = pd.DataFrame(index=self.times)
self.day_of_year = self.solar_angles.index.dayofyear
def solar_time(self, datetime_val, day_of_year):
b = (day_of_year - 81) * 2 * math.pi / 364
eot = 9.87 * math.sin(2 * b) - 7.53 * math.cos(b) - 1.5 * math.sin(b)
self.eot.append(eot)
# Calculate Local Solar Time (LST)
lst_hour = datetime_val.hour
lst_minute = datetime_val.minute
lst_second = datetime_val.second
lst = lst_hour + lst_minute / 60 + lst_second / 3600
# Calculate Apparent Solar Time (AST) in decimal hours
ast_decimal = lst + eot / 60 + self.longitude_correction / 60
ast_hours = int(ast_decimal) % 24 # Adjust hours to fit within 023 range
ast_minutes = round((ast_decimal - ast_hours) * 60)
# Ensure ast_minutes is within valid range
if ast_minutes == 60:
ast_hours += 1
ast_minutes = 0
elif ast_minutes < 0:
ast_minutes = 0
ast_time = datetime(year=datetime_val.year, month=datetime_val.month, day=datetime_val.day,
hour=ast_hours, minute=ast_minutes)
self.ast.append(ast_time)
return ast_time
def declination_angle(self, day_of_year):
declination = 23.45 * math.sin(math.radians(360 / 365 * (284 + day_of_year)))
declination_radian = math.radians(declination)
self.declinations.append(declination)
return declination_radian
def hour_angle(self, ast_time):
hour_angle = ((ast_time.hour * 60 + ast_time.minute) - 720) / 4
hour_angle_radian = math.radians(hour_angle)
self.hour_angles.append(hour_angle)
return hour_angle_radian
def solar_altitude(self, declination_radian, hour_angle_radian):
solar_altitude_radians = math.asin(math.cos(self.location_latitude_rad) * math.cos(declination_radian) *
math.cos(hour_angle_radian) + math.sin(self.location_latitude_rad) *
math.sin(declination_radian))
solar_altitude = math.degrees(solar_altitude_radians)
self.solar_altitudes.append(solar_altitude)
return solar_altitude_radians
def zenith(self, solar_altitude_radians):
solar_altitude = math.degrees(solar_altitude_radians)
zenith_degree = 90 - solar_altitude
zenith_radian = math.radians(zenith_degree)
self.zeniths.append(zenith_degree)
return zenith_radian
def solar_azimuth_analytical(self, hourangle, declination, zenith):
numer = (math.cos(zenith) * math.sin(self.location_latitude_rad) - math.sin(declination))
denom = (math.sin(zenith) * math.cos(self.location_latitude_rad))
if math.isclose(denom, 0.0, abs_tol=1e-8):
cos_azi = 1.0
else:
cos_azi = numer / denom
cos_azi = max(-1.0, min(1.0, cos_azi))
sign_ha = math.copysign(1, hourangle)
solar_azimuth_radians = sign_ha * math.acos(cos_azi) + math.pi
solar_azimuth_degrees = math.degrees(solar_azimuth_radians)
self.solar_azimuths.append(solar_azimuth_degrees)
return solar_azimuth_radians
def incident_angle(self, solar_altitude_radians, solar_azimuth_radians):
incident_radian = math.acos(math.cos(solar_altitude_radians) *
math.cos(abs(solar_azimuth_radians - self.surface_azimuth_rad)) *
math.sin(self.tilt_angle_rad) + math.sin(solar_altitude_radians) *
math.cos(self.tilt_angle_rad))
incident_angle_degrees = math.degrees(incident_radian)
self.incidents.append(incident_angle_degrees)
return incident_radian
def dni_extra(self, day_of_year, zenith_radian):
i_on = self.solar_constant * (1 + 0.033 * math.cos(math.radians(360 * day_of_year / 365)))
i_oh = i_on * max(math.cos(zenith_radian), self.min_cos_zenith)
self.i_on.append(i_on)
self.i_oh.append(i_oh)
return i_on, i_oh
def clearness_index(self, ghi, i_oh):
k_t = ghi / i_oh
k_t = max(0, k_t)
k_t = min(self.maximum_clearness_index, k_t)
return k_t
def diffuse_fraction(self, k_t, zenith):
if k_t <= 0.22:
fraction_diffuse = 1 - 0.09 * k_t
elif k_t <= 0.8:
fraction_diffuse = (0.9511 - 0.1604 * k_t + 4.388 * k_t ** 2 - 16.638 * k_t ** 3 + 12.336 * k_t ** 4)
else:
fraction_diffuse = 0.165
if zenith > self.maximum_zenith_angle:
fraction_diffuse = 1
return fraction_diffuse
def radiation_components_horizontal(self, ghi, fraction_diffuse, zenith):
diffuse_horizontal = ghi * fraction_diffuse
dni = (ghi - diffuse_horizontal) / math.cos(math.radians(zenith))
if zenith > self.maximum_zenith_angle or dni < 0:
dni = 0
return diffuse_horizontal, dni
def radiation_components_tilted(self, diffuse_horizontal, dni, incident_angle):
beam_tilted = dni * math.cos(math.radians(incident_angle))
beam_tilted = max(beam_tilted, 0)
diffuse_tilted = diffuse_horizontal * ((1 + math.cos(math.radians(self.tilt_angle))) / 2)
total_radiation_tilted = beam_tilted + diffuse_tilted
return total_radiation_tilted
def solar_angles_calculator(self, csv_output=False):
for i in range(len(self.times)):
datetime_val = self.times[i]
day_of_year = self.day_of_year[i]
declination_radians = self.declination_angle(day_of_year)
ast_time = self.solar_time(datetime_val, day_of_year)
hour_angle_radians = self.hour_angle(ast_time)
solar_altitude_radians = self.solar_altitude(declination_radians, hour_angle_radians)
zenith_radians = self.zenith(solar_altitude_radians)
solar_azimuth_radians = self.solar_azimuth_analytical(hour_angle_radians, declination_radians, zenith_radians)
self.incident_angle(solar_altitude_radians, solar_azimuth_radians)
self.dni_extra(day_of_year=day_of_year, zenith_radian=zenith_radians)
self.solar_angles['DateTime'] = self.times
self.solar_angles['AST'] = self.ast
self.solar_angles['hour angle'] = self.hour_angles
self.solar_angles['eot'] = self.eot
self.solar_angles['declination angle'] = self.declinations
self.solar_angles['solar altitude'] = self.solar_altitudes
self.solar_angles['zenith'] = self.zeniths
self.solar_angles['solar azimuth'] = self.solar_azimuths
self.solar_angles['incident angle'] = self.incidents
self.solar_angles['extraterrestrial normal radiation (Wh/m2)'] = self.i_on
self.solar_angles['extraterrestrial radiation on horizontal (Wh/m2)'] = self.i_oh
if csv_output:
self.solar_angles.to_csv('solar_angles_new.csv')
def tilted_irradiance_calculator(self):
if self.solar_angles.empty:
self.solar_angles_calculator()
for building in self.city.buildings:
hourly_tilted_irradiance = []
roof_ghi = building.roofs[0].global_irradiance[cte.HOUR]
for i in range(len(roof_ghi)):
k_t = self.clearness_index(ghi=roof_ghi[i], i_oh=self.i_oh[i])
fraction_diffuse = self.diffuse_fraction(k_t, self.zeniths[i])
diffuse_horizontal, dni = self.radiation_components_horizontal(ghi=roof_ghi[i],
fraction_diffuse=fraction_diffuse,
zenith=self.zeniths[i])
hourly_tilted_irradiance.append(int(self.radiation_components_tilted(diffuse_horizontal=diffuse_horizontal,
dni=dni,
incident_angle=self.incidents[i])))
building.roofs[0].global_irradiance_tilted[cte.HOUR] = hourly_tilted_irradiance
building.roofs[0].global_irradiance_tilted[cte.MONTH] = (MonthlyValues.get_total_month(
building.roofs[0].global_irradiance_tilted[cte.HOUR]))
building.roofs[0].global_irradiance_tilted[cte.YEAR] = [sum(building.roofs[0].global_irradiance_tilted[cte.MONTH])]