diff --git a/hub/city_model_structure/building_demand/surface.py b/hub/city_model_structure/building_demand/surface.py index 2cd42755..f54d21e1 100644 --- a/hub/city_model_structure/building_demand/surface.py +++ b/hub/city_model_structure/building_demand/surface.py @@ -46,6 +46,7 @@ class Surface: self._vegetation = None self._percentage_shared = None self._solar_collectors_area_reduction_factor = None + self._global_irradiance_tilted = {} @property def name(self): @@ -384,3 +385,19 @@ class Surface: :param value: float """ self._solar_collectors_area_reduction_factor = value + + @property + def global_irradiance_tilted(self) -> dict: + """ + Get global irradiance on a tilted surface in J/m2 + :return: dict + """ + return self._global_irradiance_tilted + + @global_irradiance_tilted.setter + def global_irradiance_tilted(self, value): + """ + Set global irradiance on a tilted surface in J/m2 + :param value: dict + """ + self._global_irradiance_tilted = value diff --git a/hub/city_model_structure/city_object.py b/hub/city_model_structure/city_object.py index c18a4aa0..0d281b65 100644 --- a/hub/city_model_structure/city_object.py +++ b/hub/city_model_structure/city_object.py @@ -41,9 +41,10 @@ class CityObject: self._ground_temperature = {} self._global_horizontal = {} self._diffuse = {} - self._beam = {} + self._direct_normal = {} self._sensors = [] self._neighbours = None + self._beam = {} @property def level_of_detail(self) -> LevelOfDetail: @@ -238,20 +239,20 @@ class CityObject: self._diffuse = value @property - def beam(self) -> dict: + def direct_normal(self) -> dict: """ Get beam radiation surrounding the city object in J/m2 :return: dict{dict{[float]}} """ - return self._beam + return self._direct_normal - @beam.setter - def beam(self, value): + @direct_normal.setter + def direct_normal(self, value): """ Set beam radiation surrounding the city object in J/m2 :param value: dict{dict{[float]}} """ - self._beam = value + self._direct_normal = value @property def lower_corner(self): @@ -302,3 +303,19 @@ class CityObject: Set the list of neighbour_objects and their properties associated to the current city_object """ self._neighbours = value + + @property + def beam(self) -> dict: + """ + Get beam radiation surrounding the city object in J/m2 + :return: dict{dict{[float]}} + """ + return self._beam + + @beam.setter + def beam(self, value): + """ + Set beam radiation surrounding the city object in J/m2 + :param value: dict{dict{[float]}} + """ + self._beam = value diff --git a/hub/city_model_structure/energy_systems/pv_generation_system.py b/hub/city_model_structure/energy_systems/pv_generation_system.py index 13409c7e..ddbc847a 100644 --- a/hub/city_model_structure/energy_systems/pv_generation_system.py +++ b/hub/city_model_structure/energy_systems/pv_generation_system.py @@ -26,6 +26,10 @@ class PvGenerationSystem(GenerationSystem): self._width = None self._height = None self._electricity_power = None + self._tilt_angle = None + self._surface_azimuth = None + self._solar_altitude_angle = None + self._solar_azimuth_angle = None @property def nominal_electricity_output(self): @@ -202,3 +206,35 @@ class PvGenerationSystem(GenerationSystem): :param value: float """ self._electricity_power = value + + @property + def tilt_angle(self): + """ + Get tilt angle of PV system in degrees + :return: float + """ + return self._tilt_angle + + @tilt_angle.setter + def tilt_angle(self, value): + """ + Set PV system tilt angle in degrees + :param value: float + """ + self._tilt_angle = value + + @property + def surface_azimuth(self): + """ + Get surface azimuth angle of PV system in degrees. 0 is North + :return: float + """ + return self._surface_azimuth + + @surface_azimuth.setter + def surface_azimuth(self, value): + """ + Set PV system tilt angle in degrees + :param value: float + """ + self._surface_azimuth = value diff --git a/hub/exports/formats/simplified_radiosity_algorithm.py b/hub/exports/formats/simplified_radiosity_algorithm.py index f9ea7f1d..25189419 100644 --- a/hub/exports/formats/simplified_radiosity_algorithm.py +++ b/hub/exports/formats/simplified_radiosity_algorithm.py @@ -67,7 +67,7 @@ class SimplifiedRadiosityAlgorithm: i = (total_days + day - 1) * 24 + hour - 1 representative_building = self._city.buildings[0] _global = representative_building.diffuse[cte.HOUR][i] / cte.WATTS_HOUR_TO_JULES - _beam = representative_building.beam[cte.HOUR][i] / cte.WATTS_HOUR_TO_JULES + _beam = representative_building.direct_normal[cte.HOUR][i] / cte.WATTS_HOUR_TO_JULES content += f'{day} {month} {hour} {_global} {_beam}\n' with open(file, 'w', encoding='utf-8') as file: file.write(content) diff --git a/hub/imports/energy_systems/montreal_future_energy_systems_parameters.py b/hub/imports/energy_systems/montreal_future_energy_systems_parameters.py index 8612f845..27c2f1cd 100644 --- a/hub/imports/energy_systems/montreal_future_energy_systems_parameters.py +++ b/hub/imports/energy_systems/montreal_future_energy_systems_parameters.py @@ -82,8 +82,7 @@ class MontrealFutureEnergySystemParameters: return _generic_energy_systems - @staticmethod - def _create_generation_systems(archetype_system): + def _create_generation_systems(self, archetype_system): _generation_systems = [] archetype_generation_systems = archetype_system.generation_systems if archetype_generation_systems is not None: @@ -107,6 +106,7 @@ class MontrealFutureEnergySystemParameters: _generation_system.cell_temperature_coefficient = archetype_generation_system.cell_temperature_coefficient _generation_system.width = archetype_generation_system.width _generation_system.height = archetype_generation_system.height + _generation_system.tilt_angle = self._city.latitude _generic_storage_system = None if archetype_generation_system.energy_storage_systems is not None: _generic_storage_system = ElectricalStorageSystem() diff --git a/hub/imports/results/simplified_radiosity_algorithm.py b/hub/imports/results/simplified_radiosity_algorithm.py index 9b09b5f8..3824e7ce 100644 --- a/hub/imports/results/simplified_radiosity_algorithm.py +++ b/hub/imports/results/simplified_radiosity_algorithm.py @@ -34,7 +34,7 @@ class SimplifiedRadiosityAlgorithm: for key in self._results: _irradiance = {} header_name = key.split(':') - result = [x for x in self._results[key]] + result = [x * cte.WATTS_HOUR_TO_JULES for x in self._results[key]] city_object_name = header_name[1] building = self._city.city_object(city_object_name) surface_id = header_name[2] diff --git a/hub/imports/weather/epw_weather_parameters.py b/hub/imports/weather/epw_weather_parameters.py index c2ed2e8f..d8f3a002 100644 --- a/hub/imports/weather/epw_weather_parameters.py +++ b/hub/imports/weather/epw_weather_parameters.py @@ -114,10 +114,14 @@ class EpwWeatherParameters: for x in self._weather_values['global_horizontal_radiation_wh_m2']] building.diffuse[cte.HOUR] = [x * cte.WATTS_HOUR_TO_JULES for x in self._weather_values['diffuse_horizontal_radiation_wh_m2']] - building.beam[cte.HOUR] = [x * cte.WATTS_HOUR_TO_JULES - for x in self._weather_values['direct_normal_radiation_wh_m2']] + building.direct_normal[cte.HOUR] = [x * cte.WATTS_HOUR_TO_JULES + for x in self._weather_values['direct_normal_radiation_wh_m2']] + building.beam[cte.HOUR] = [building.global_horizontal[cte.HOUR][i] - + building.diffuse[cte.HOUR][i] + for i in range(len(building.global_horizontal[cte.HOUR]))] building.cold_water_temperature[cte.HOUR] = wh().cold_water_temperature(building.external_temperature[cte.HOUR]) + # create the monthly and yearly values out of the hourly for building in self._city.buildings: building.external_temperature[cte.MONTH] = \ diff --git a/main.py b/main.py index 4af9643a..86ffeb24 100644 --- a/main.py +++ b/main.py @@ -1,7 +1,6 @@ from scripts.geojson_creator import process_geojson from pathlib import Path import subprocess -from scripts.ep_run_enrich import energy_plus_workflow from hub.imports.geometry_factory import GeometryFactory from hub.helpers.dictionaries import Dictionaries from hub.imports.construction_factory import ConstructionFactory @@ -18,8 +17,11 @@ from scripts.costs.cost import Cost from scripts.costs.constants import SKIN_RETROFIT_AND_SYSTEM_RETROFIT_AND_PV, SYSTEM_RETROFIT_AND_PV import hub.helpers.constants as cte from hub.exports.exports_factory import ExportsFactory +import csv +from scripts.solar_angles import CitySolarAngles +from scripts.radiation_tilted import RadiationTilted # Specify the GeoJSON file path -# geojson_file = process_geojson(x=-73.5953602192335, y=45.492414530022515, diff=0.001) +geojson_file = process_geojson(x=-73.5681295982132, y=45.49218262677643, diff=0.001) file_path = (Path(__file__).parent / 'input_files' / 'output_buildings.geojson') # Specify the output path for the PDF file output_path = (Path(__file__).parent / 'out_files').resolve() @@ -30,9 +32,8 @@ city = GeometryFactory('geojson', year_of_construction_field='year_of_construction', function_field='function', function_to_hub=Dictionaries().montreal_function_to_hub_function).city -# Enrich city data +# # Enrich city data ConstructionFactory('nrcan', city).enrich() - UsageFactory('nrcan', city).enrich() WeatherFactory('epw', city).enrich() energy_plus_workflow(city) diff --git a/scripts/pv_sizing_and_simulation.py b/scripts/pv_sizing_and_simulation.py new file mode 100644 index 00000000..cd8e1673 --- /dev/null +++ b/scripts/pv_sizing_and_simulation.py @@ -0,0 +1,21 @@ +from scripts.radiation_tilted import RadiationTilted +from scripts.solar_angles import CitySolarAngles +import hub.helpers.constants as cte + + +class PVSizingSimulation: + def __init__(self, city, tilt_angle): + self.city = city + self.tilt_angle = tilt_angle + self.solar_angles = CitySolarAngles(file_name=self.city.name, + location_latitude=self.city.latitude, + location_longitude=self.city.longitude, + tilt_angle=self.tilt_angle, + surface_azimuth_angle=180, + standard_meridian=-75).calculate + self.enrich_radiation_data() + + def enrich_radiation_data(self): + for building in self.city.buildings: + roof_horizontal = [x / cte.WATTS_HOUR_TO_JULES for x in building.roofs[0].global_irradiance[cte.HOUR]] + RadiationTilted(building, self.solar_angles, tilt_angle=self.tilt_angle, ghi=roof_horizontal).enrich() diff --git a/scripts/radiation_tilted.py b/scripts/radiation_tilted.py new file mode 100644 index 00000000..60c8d7de --- /dev/null +++ b/scripts/radiation_tilted.py @@ -0,0 +1,109 @@ +import pandas as pd +import math +import hub.helpers.constants as cte +from hub.helpers.monthly_values import MonthlyValues + + +class RadiationTilted: + def __init__(self, building, solar_angles, tilt_angle, ghi, solar_constant=1366.1, maximum_clearness_index=1, + min_cos_zenith=0.065, maximum_zenith_angle=87): + self.building = building + self.ghi = ghi + self.tilt_angle = tilt_angle + self.zeniths = solar_angles['zenith'].tolist()[:-1] + self.incidents = solar_angles['incident angle'].tolist()[:-1] + self.date_time = solar_angles['DateTime'].tolist()[:-1] + data = {'DateTime': self.date_time, 'zenith': self.zeniths, 'incident angle': self.incidents, 'ghi': self.ghi} + self.df = pd.DataFrame(data) + self.df['DateTime'] = pd.to_datetime(self.df['DateTime']) + self.df.set_index('DateTime', inplace=True) + 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 + self.i_on = [] + self.i_oh = [] + self.k_t = [] + self.fraction_diffuse = [] + self.diffuse_horizontal = [] + self.beam_horizontal = [] + self.dni = [] + self.beam_tilted = [] + self.diffuse_tilted = [] + self.total_radiation_tilted = [] + self.calculate() + + def dni_extra(self): + for i in range(len(self.df)): + self.i_on.append(self.solar_constant * (1 + 0.033 * math.cos(math.radians(360 * self.df.index.dayofyear[i] / 365)))) + + self.df['extraterrestrial normal radiation (Wh/m2)'] = self.i_on + + def clearness_index(self): + for i in range(len(self.df)): + self.i_oh.append(self.i_on[i] * max(math.cos(math.radians(self.zeniths[i])), self.min_cos_zenith)) + self.k_t.append(self.ghi[i] / self.i_oh[i]) + self.k_t[i] = max(0, self.k_t[i]) + self.k_t[i] = min(self.maximum_clearness_index, self.k_t[i]) + self.df['extraterrestrial radiation on horizontal (Wh/m2)'] = self.i_oh + self.df['clearness index'] = self.k_t + + def diffuse_fraction(self): + for i in range(len(self.df)): + if self.k_t[i] <= 0.22: + self.fraction_diffuse.append(1 - 0.09 * self.k_t[i]) + elif self.k_t[i] <= 0.8: + self.fraction_diffuse.append(0.9511 - 0.1604 * self.k_t[i] + 4.388 * self.k_t[i] ** 2 - + 16.638 * self.k_t[i] ** 3 + 12.336 * self.k_t[i] ** 4) + else: + self.fraction_diffuse.append(0.165) + if self.zeniths[i] > self.maximum_zenith_angle: + self.fraction_diffuse[i] = 1 + + self.df['diffuse fraction'] = self.fraction_diffuse + + def radiation_components_horizontal(self): + for i in range(len(self.df)): + self.diffuse_horizontal.append(self.ghi[i] * self.fraction_diffuse[i]) + self.beam_horizontal.append(self.ghi[i] - self.diffuse_horizontal[i]) + self.dni.append((self.ghi[i] - self.diffuse_horizontal[i]) / math.cos(math.radians(self.zeniths[i]))) + if self.zeniths[i] > self.maximum_zenith_angle or self.dni[i] < 0: + self.dni[i] = 0 + + self.df['diffuse horizontal (Wh/m2)'] = self.diffuse_horizontal + self.df['dni (Wh/m2)'] = self.dni + self.df['beam horizontal (Wh/m2)'] = self.beam_horizontal + + def radiation_components_tilted(self): + for i in range(len(self.df)): + self.beam_tilted.append(self.dni[i] * math.cos(math.radians(self.incidents[i]))) + self.beam_tilted[i] = max(self.beam_tilted[i], 0) + self.diffuse_tilted.append(self.diffuse_horizontal[i] * ((1 + math.cos(math.radians(self.tilt_angle))) / 2)) + self.total_radiation_tilted.append(self.beam_tilted[i] + self.diffuse_tilted[i]) + + self.df['beam tilted (Wh/m2)'] = self.beam_tilted + self.df['diffuse tilted (Wh/m2)'] = self.diffuse_tilted + self.df['total radiation tilted (Wh/m2)'] = self.total_radiation_tilted + + def calculate(self) -> pd.DataFrame: + self.dni_extra() + self.clearness_index() + self.diffuse_fraction() + self.radiation_components_horizontal() + self.radiation_components_tilted() + return self.df + + def enrich(self): + tilted_radiation = self.total_radiation_tilted + print(len(tilted_radiation)) + self.building.roofs[0].global_irradiance_tilted[cte.HOUR] = [x * cte.WATTS_HOUR_TO_JULES for x in + tilted_radiation] + self.building.roofs[0].global_irradiance_tilted[cte.HOUR] = [x * cte.WATTS_HOUR_TO_JULES for x in + tilted_radiation] + self.building.roofs[0].global_irradiance_tilted[cte.MONTH] = ( + MonthlyValues.get_total_month(self.building.roofs[0].global_irradiance_tilted[cte.HOUR])) + self.building.roofs[0].global_irradiance_tilted[cte.YEAR] = \ + [sum(self.building.roofs[0].global_irradiance_tilted[cte.MONTH])] + + + diff --git a/scripts/solar_angles.py b/scripts/solar_angles.py new file mode 100644 index 00000000..d65fcf58 --- /dev/null +++ b/scripts/solar_angles.py @@ -0,0 +1,146 @@ +import math +import pandas as pd +from datetime import datetime +from pathlib import Path + + +class CitySolarAngles: + def __init__(self, file_name, location_latitude, location_longitude, tilt_angle, surface_azimuth_angle, + standard_meridian=-75): + self.file_name = file_name + self.location_latitude = location_latitude + self.location_longitude = location_longitude + self.location_latitude_rad = math.radians(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 = (location_longitude - standard_meridian) * 4 + self.timezone = 'Etc/GMT+5' + + self.eot = [] + self.ast = [] + self.hour_angles = [] + self.declinations = [] + self.solar_altitudes = [] + self.solar_azimuths = [] + self.zeniths = [] + self.incidents = [] + self.beam_tilted = [] + self.factor = [] + self.times = pd.date_range(start='2023-01-01', end='2024-01-01', freq='H', tz=self.timezone) + self.df = pd.DataFrame(index=self.times) + self.day_of_year = self.df.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) + 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 + + @property + def calculate(self) -> pd.DataFrame: + 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) + incident_angle_radian = self.incident_angle(solar_altitude_radians, solar_azimuth_radians) + + self.df['DateTime'] = self.times + self.df['AST'] = self.ast + self.df['hour angle'] = self.hour_angles + self.df['eot'] = self.eot + self.df['declination angle'] = self.declinations + self.df['solar altitude'] = self.solar_altitudes + self.df['zenith'] = self.zeniths + self.df['solar azimuth'] = self.solar_azimuths + self.df['incident angle'] = self.incidents + + return self.df + + + + + + + + diff --git a/tests/test_construction_factory.py b/tests/test_construction_factory.py index 710894bd..a340594a 100644 --- a/tests/test_construction_factory.py +++ b/tests/test_construction_factory.py @@ -81,7 +81,7 @@ class TestConstructionFactory(TestCase): self.assertEqual(len(building.external_temperature), 0, 'building external temperature is calculated') self.assertEqual(len(building.global_horizontal), 0, 'building global horizontal is calculated') self.assertEqual(len(building.diffuse), 0, 'building diffuse is calculated') - self.assertEqual(len(building.beam), 0, 'building beam is calculated') + self.assertEqual(len(building.direct_normal), 0, 'building beam is calculated') self.assertIsNotNone(building.lower_corner, 'building lower corner is none') self.assertEqual(len(building.sensors), 0, 'building sensors are assigned') self.assertIsNotNone(building.internal_zones, 'no internal zones created') diff --git a/tests/test_geometry_factory.py b/tests/test_geometry_factory.py index 3b5bd8f8..d66c815f 100644 --- a/tests/test_geometry_factory.py +++ b/tests/test_geometry_factory.py @@ -52,7 +52,7 @@ class TestGeometryFactory(TestCase): self.assertEqual(len(building.external_temperature), 0, 'building external temperature is calculated') self.assertEqual(len(building.global_horizontal), 0, 'building global horizontal is calculated') self.assertEqual(len(building.diffuse), 0, 'building diffuse is calculated') - self.assertEqual(len(building.beam), 0, 'building beam is calculated') + self.assertEqual(len(building.direct_normal), 0, 'building beam is calculated') self.assertIsNotNone(building.lower_corner, 'building lower corner is none') self.assertEqual(len(building.sensors), 0, 'building sensors are assigned') self.assertIsNotNone(building.internal_zones, 'no internal zones created') diff --git a/tests/test_usage_factory.py b/tests/test_usage_factory.py index 6b6d821a..49cb45db 100644 --- a/tests/test_usage_factory.py +++ b/tests/test_usage_factory.py @@ -44,7 +44,7 @@ class TestUsageFactory(TestCase): self.assertEqual(len(building.external_temperature), 0, 'building external temperature is calculated') self.assertEqual(len(building.global_horizontal), 0, 'building global horizontal is calculated') self.assertEqual(len(building.diffuse), 0, 'building diffuse is calculated') - self.assertEqual(len(building.beam), 0, 'building beam is calculated') + self.assertEqual(len(building.direct_normal), 0, 'building beam is calculated') self.assertIsNotNone(building.lower_corner, 'building lower corner is none') self.assertEqual(len(building.sensors), 0, 'building sensors are assigned') self.assertIsNotNone(building.internal_zones, 'no internal zones created')