"""
solar collectors
"""
from __future__ import division
import numpy as np
import pandas as pd
from math import *
import re
from cea.utilities import epwreader
from cea.utilities import solar_equations
__author__ = "Jimeno A. Fonseca"
__copyright__ = "Copyright 2015, Architecture and Building Systems - ETH Zurich"
__credits__ = ["Jimeno A. Fonseca"]
__license__ = "MIT"
__version__ = "0.1"
__maintainer__ = "Daren Thomas"
__email__ = "cea@arch.ethz.ch"
__status__ = "Production"
# SC heat generation
[docs]def calc_SC(locator, sensors_data, radiation, latitude, longitude, year, gv, weather_path):
# weather data
weather_data = epwreader.epw_reader(weather_path)
# solar properties
g, Sz, Az, ha, trr_mean, worst_sh, worst_Az = solar_equations.calc_sun_properties(latitude, longitude, weather_data,
gv)
# read radiation file
hourly_data = pd.read_csv(radiation)
# get only datapoints with production beyond min_production
Max_Isol = hourly_data.total.max()
Min_Isol = Max_Isol * gv.min_production # 80% of the local average maximum in the area
sensors_data_clean = sensors_data[sensors_data["total"] > Min_Isol]
radiation_clean = radiation.loc[radiation['sensor_id'].isin(sensors_data_clean.sensor_id)]
# Calculate the heights of all buildings for length of vertical pipes
height = locator.get_total_demand().height.sum()
# calculate optimal angle and tilt for panels
optimal_angle_and_tilt(sensors_data, latitude, worst_sh, worst_Az, trr_mean, gv.grid_side,
gv.module_lenght_SC, gv.angle_north, Min_Isol, Max_Isol)
Number_groups, hourlydata_groups, number_points, prop_observers = calc_groups(radiation_clean, sensors_data_clean)
result, Final = SC_generation(gv.type_SCpanel, hourlydata_groups, prop_observers, number_points, g, Sz, Az, ha,
latitude,
gv.Tin, height)
Final.to_csv(locator.solar_collectors_result(), index=True, float_format='%.2f')
return
[docs]def SC_generation(type_SCpanel, group_radiation, prop_observers, number_points, weather_data, g, Sz, Az, ha, latitude,
Tin, height):
# get properties of the panel to evaluate
n0, c1, c2, mB0_r, mB_max_r, mB_min_r, C_eff, t_max, IAM_d, Aratio, Apanel, dP1, dP2, dP3, dP4 = calc_properties_SC(
type_SCpanel)
Area_a = Aratio * Apanel
listgroups = number_points.count() # counter from the vector with number of points
listresults = [None] * listgroups
listresults_perarea = [None] * listgroups
listareasgroups = [None] * listgroups
Sum_mcp = np.zeros(8760)
Sum_qout = np.zeros(8760)
Sum_Eaux = np.zeros(8760)
Sum_qloss = np.zeros(8760)
Tin_array = np.zeros(8760) + Tin
Sum_Area_m = (prop_observers['area_netpv'] * number_points).sum()
lv = 2 # grid lenght module length
Le = (2 * lv * number_points.sum()) / (Sum_Area_m * Aratio)
Li = 2 * height / (Sum_Area_m * Aratio)
Leq = Li + Le # in m/m2
if type_SCpanel == 2: # for vaccum tubes
Nseg = 100 # default number of subsdivisions for the calculation
else:
Nseg = 10 # default number of subsdivisions for the calculation
for group in range(listgroups):
teta_z = prop_observers.loc[group, 'aspect'] # azimuth of paneles of group
Area_group = prop_observers.loc[group, 'area_netpv'] * number_points[group]
tilt_angle = prop_observers.loc[group, 'slope'] # tilt angle of panels
radiation = pd.DataFrame({'I_sol': group_radiation[group]}) # choose vector with all values of Isol
radiation['I_diffuse'] = weather_data.ratio_diffhout * radiation.I_sol # calculate diffuse radiation
radiation['I_direct'] = radiation['I_sol'] - radiation['I_diffuse'] # direct radaition
# calculate angle modifiers,
IAM_b = calc_anglemodifierSC(Az, g, ha, teta_z, tilt_angle, type_SCpanel,
latitude, Sz) # direct angle modifier
listresults[group] = Calc_SC_module2(radiation, tilt_angle, IAM_b, radiation.I_direct,
radiation.I_diffuse, weather_data.drybulb_C,
n0, c1, c2, mB0_r, mB_max_r, mB_min_r, C_eff, t_max, IAM_d, Area_a,
dP1, dP2, dP3, dP4, Tin, Leq, Le, Nseg)
K = Area_group / Apanel
listresults[group][5] = listresults[group][5] * K
listresults[group][0] = listresults[group][0] * K
listresults[group][1] = listresults[group][1] * K
listresults[group][2] = listresults[group][2] * K
listareasgroups[group] = Area_group
for group in range(listgroups):
mcp_array = listresults[group][5]
qloss_array = listresults[group][0]
qout_array = listresults[group][1]
Eaux_array = listresults[group][2]
Sum_qout = Sum_qout + qout_array
Sum_Eaux = Sum_Eaux + Eaux_array
Sum_qloss = Sum_qloss + qloss_array
Sum_mcp = Sum_mcp + mcp_array
Tout_group = (Sum_qout / Sum_mcp) + Tin # in C
Final = pd.DataFrame(
{'Qsc_Kw': Sum_qout, 'Tscs': Tin_array, 'Tscr': Tout_group, 'mcp_kW/C': Sum_mcp, 'Eaux_kW': Sum_Eaux,
'Qsc_l_KWH': Sum_qloss, 'Area': sum(listareasgroups)}, index=range(8760))
return listresults, Final
[docs]def calc_groups(Clean_hourly, observers_fin):
# calculate number of optima groups as number of optimal combiantions.
groups_ob = Clean_hourly.groupby(['CATB', 'CATGB', 'CATteta_z'])
hourlydata_groups = groups_ob.mean().reset_index()
hourlydata_groups = pd.DataFrame(hourlydata_groups)
Number_pointsgroup = groups_ob.size().reset_index()
number_points = Number_pointsgroup[0]
groups_ob = observers_fin.groupby(['CATB', 'CATGB', 'CATteta_z'])
prop_observers = groups_ob.mean().reset_index()
prop_observers = pd.DataFrame(prop_observers)
Number_groups = groups_ob.size().count()
hourlydata_groups = hourlydata_groups.drop({'ID', 'GB', 'grid_code', 'pointid', 'array_s', 'area_netpv', 'aspect',
'slope', 'CATB', 'CATGB', 'CATteta_z'}, axis=1).transpose().reindex(
axis=1) # vector with radiation points of group
hourlydata_groups['newindex'] = hourlydata_groups.index
hourlydata_groups['newindex'] = hourlydata_groups.newindex.apply(lambda x: re.findall('\d+', x))
hourlydata_groups.index = range(8760)
for hour in range(8760):
hourlydata_groups.loc[hour, 'newindex'] = int(hourlydata_groups.loc[hour, 'newindex'][0])
hourlydata_groups.set_index('newindex', inplace=True)
hourlydata_groups.sort_index(inplace=True)
hourlydata_groups.index = range(8760)
return Number_groups, hourlydata_groups, number_points, prop_observers
[docs]def Calc_SC_module2(radiation, tilt_angle, IAM_b_vector, I_direct_vector, I_diffuse_vector, Te_vector, n0, c1, c2,
mB0_r,
mB_max_r, mB_min_r, C_eff, t_max, IAM_d, Area_a, dP1, dP2, dP3, dP4, Tin, Leq, Le, Nseg):
# panel to store the results per flow
# method with no condensaiton gains, no wind or long-wave dependency, sky factor set to zero.
# calculate radiation part
# local variables
Cpwg = 3680 # J/kgK water grlycol specific heat
maxmsc = mB_max_r * Area_a / 3600
# Do the calculation of every time step for every possible flow condition
# get states where highly performing values are obtained.
mopt = 0 # para iniciar
temperature_out = [np.zeros(8760), np.zeros(8760), np.zeros(8760), np.zeros(8760), np.zeros(8760), np.zeros(8760)]
temperature_in = [np.zeros(8760), np.zeros(8760), np.zeros(8760), np.zeros(8760), np.zeros(8760), np.zeros(8760)]
supply_out = [np.zeros(8760), np.zeros(8760), np.zeros(8760), np.zeros(8760), np.zeros(8760), np.zeros(8760)]
supply_losses = [np.zeros(8760), np.zeros(8760), np.zeros(8760), np.zeros(8760), np.zeros(8760), np.zeros(8760)]
Auxiliary = [np.zeros(8760), np.zeros(8760), np.zeros(8760), np.zeros(8760), np.zeros(8760), np.zeros(8760)]
temperature_m = [np.zeros(8760), np.zeros(8760), np.zeros(8760), np.zeros(8760), np.zeros(8760), np.zeros(8760)]
specific_flows = [np.zeros(8760), (np.zeros(8760) + mB0_r) * Area_a / 3600,
(np.zeros(8760) + mB_max_r) * Area_a / 3600,
(np.zeros(8760) + mB_min_r) * Area_a / 3600, np.zeros(8760), np.zeros(8760)] # in kg/s
specific_pressurelosses = [np.zeros(8760), (np.zeros(8760) + dP2) * Area_a, (np.zeros(8760) + dP3) * Area_a,
(np.zeros(8760) + dP4) * Area_a, np.zeros(8760), np.zeros(8760)] # in Pa
supply_out_pre = np.zeros(8760)
supply_out_total = np.zeros(8760)
mcp = np.zeros(8760)
# calculate net radiant heat (absorbed)
tilt = radians(tilt_angle)
qrad_vector = np.vectorize(calc_qrad)(n0, IAM_b_vector, I_direct_vector, IAM_d, I_diffuse_vector,
tilt) # in W/m2 is a mean of the group
counter = 0
Flag = False
Flag2 = False
for flow in range(6):
Mo = 1
TIME0 = 0
DELT = 1 # timestep 1 hour
delts = DELT * 3600 # convert time step in seconds
Tfl = np.zeros([3, 1]) # create vector
DT = np.zeros([3, 1])
Tabs = np.zeros([3, 1])
STORED = np.zeros([600, 1])
TflA = np.zeros([600, 1])
TflB = np.zeros([600, 1])
TabsB = np.zeros([600, 1])
TabsA = np.zeros([600, 1])
qgainSeg = np.zeros([100, 1])
for time in range(8760):
Mfl = specific_flows[flow][time]
if time < TIME0 + DELT / 2:
for Iseg in range(101, 501): # 400 points with the data
STORED[Iseg] = Tin
else:
for Iseg in range(1, Nseg): # 400 points with the data
STORED[100 + Iseg] = STORED[200 + Iseg]
STORED[300 + Iseg] = STORED[400 + Iseg]
# calculate stability criteria
if Mfl > 0:
stabcrit = Mfl * Cpwg * Nseg * (DELT * 3600) / (C_eff * Area_a)
if stabcrit <= 0.5:
print 'ERROR' + str(stabcrit) + ' ' + str(Area_a) + ' ' + str(Mfl)
Te = Te_vector[time]
qrad = qrad_vector[time]
Tfl[1] = 0 # mean absorber temp
Tabs[1] = 0 # mean absorber initial tempr
for Iseg in range(1, Nseg + 1):
Tfl[1] = Tfl[1] + STORED[100 + Iseg] / Nseg
Tabs[1] = Tabs[1] + STORED[300 + Iseg] / Nseg
# first guess for DeatT
if Mfl > 0:
Tout = Tin + (qrad - (c1 + 0.5) * (Tin - Te)) / (Mfl * Cpwg / Area_a)
Tfl[2] = (Tin + Tout) / 2
else:
Tout = Te + qrad / (c1 + 0.5)
Tfl[2] = Tout # fluid temperature same as output
DT[1] = Tfl[2] - Te
# calculate qgain with the guess
qgain = calc_qgain(Tfl, Tabs, qrad, DT, Tin, Tout, Area_a, c1, c2, Mfl, delts, Cpwg, C_eff, Te)
Aseg = Area_a / Nseg
for Iseg in range(1, Nseg + 1):
TflA[Iseg] = STORED[100 + Iseg]
TabsA[Iseg] = STORED[300 + Iseg]
if Iseg > 1:
TinSeg = ToutSeg
else:
TinSeg = Tin
if Mfl > 0 and Mo == 1:
ToutSeg = ((Mfl * Cpwg * (TinSeg + 273)) / Aseg - (C_eff * (TinSeg + 273)) / (2 * delts) + qgain +
(C_eff * (TflA[Iseg] + 273) / delts)) / (Mfl * Cpwg / Aseg + C_eff / (2 * delts))
ToutSeg = ToutSeg - 273
TflB[Iseg] = (TinSeg + ToutSeg) / 2
else:
Tfl[1] = TflA[Iseg]
Tabs[1] = TabsA[Iseg]
qgain = calc_qgain(Tfl, Tabs, qrad, DT, TinSeg, Tout, Aseg, c1, c2, Mfl, delts, Cpwg, C_eff, Te)
ToutSeg = Tout
if Mfl > 0:
TflB[Iseg] = (TinSeg + ToutSeg) / 2
ToutSeg = TflA[Iseg] + (qgain * delts) / C_eff
else:
TflB[Iseg] = ToutSeg
TflB[Iseg] = ToutSeg
qfluid = (ToutSeg - TinSeg) * Mfl * Cpwg / Aseg
qmtherm = (TflB[Iseg] - TflA[Iseg]) * C_eff / delts
qbal = qgain - qfluid - qmtherm
if abs(qbal) > 1:
time = time
qgainSeg[Iseg] = qgain # in W/m2
# the resulting energy output
qout = Mfl * Cpwg * (ToutSeg - Tin)
Tabs[2] = 0
# storage of the mean temperature
for Iseg in range(1, Nseg + 1):
STORED[200 + Iseg] = TflB[Iseg]
STORED[400 + Iseg] = TabsB[Iseg]
Tabs[2] = Tabs[2] + TabsB[Iseg] / Nseg
# outputs
temperature_out[flow][time] = ToutSeg
temperature_in[flow][time] = Tin
supply_out[flow][time] = qout / 1000 # in kW
temperature_m[flow][time] = (Tin + ToutSeg) / 2 # Mean absorber temperature at present
qgain = 0
TavgB = 0
TavgA = 0
for Iseg in range(1, Nseg + 1):
qgain = qgain + qgainSeg * Aseg # W
TavgA = TavgA + TflA[Iseg] / Nseg
TavgB = TavgB + TflB[Iseg] / Nseg
# OUT[9] = qgain/Area_a # in W/m2
qmtherm = (TavgB - TavgA) * C_eff * Area_a / delts
qbal = qgain - qmtherm - qout
# OUT[11] = qmtherm
# OUT[12] = qbal
if flow < 4:
Auxiliary[flow] = np.vectorize(calc_Eaux_SC)(specific_flows[flow], specific_pressurelosses[flow],
Leq, Area_a) # in kW
if flow == 3:
q1 = supply_out[0]
q2 = supply_out[1]
q3 = supply_out[2]
q4 = supply_out[3]
E1 = Auxiliary[0]
E2 = Auxiliary[1]
E3 = Auxiliary[2]
E4 = Auxiliary[3]
specific_flows[4], specific_pressurelosses[4] = SelectminimumenergySc(q1, q2, q3, q4, E1, E2, E3, E4, 0,
mB0_r, mB_max_r, mB_min_r, 0,
dP2, dP3, dP4, Area_a)
if flow == 4:
Auxiliary[flow] = np.vectorize(calc_Eaux_SC)(specific_flows[flow], specific_pressurelosses[flow],
Leq, Area_a) # in kW
dp5 = specific_pressurelosses[flow]
q5 = supply_out[flow]
m5 = specific_flows[flow]
##poits where load is negative
specific_flows[5], specific_pressurelosses[5] = Selectminimumenergy2(m5, q5, dp5)
if flow == 5:
supply_losses[flow] = np.vectorize(Calc_qloss_net)(specific_flows[flow], Le, Area_a, temperature_m[flow],
Te_vector, maxmsc)
supply_out_pre = supply_out[flow].copy() + supply_losses[flow].copy()
Auxiliary[flow] = np.vectorize(calc_Eaux_SC)(specific_flows[flow], specific_pressurelosses[flow],
Leq, Area_a) # in kW
supply_out_total = supply_out + 0.5 * Auxiliary[flow] - supply_losses[flow]
mcp = specific_flows[flow] * (Cpwg / 1000) # mcp in kW/c
result = [supply_losses[5], supply_out_total[5], Auxiliary[5], temperature_out[flow], temperature_in[flow], mcp]
return result
[docs]def calc_qrad(n0, IAM_b, I_direct, IAM_d, I_diffuse, tilt):
qrad = n0 * IAM_b * I_direct + n0 * IAM_d * I_diffuse * (1 + cos(tilt)) / 2
return qrad
[docs]def calc_qgain(Tfl, Tabs, qrad, DT, TinSub, Tout, Aseg, c1, c2, Mfl, delts, Cpwg, C_eff, Te):
xgain = 1
xgainmax = 100
exit = False
while exit == False:
qgain = qrad - c1 * (DT[1]) - c2 * abs(DT[1]) * DT[1]
if Mfl > 0:
Tout = ((Mfl * Cpwg * TinSub) / Aseg - (C_eff * TinSub) / (2 * delts) + qgain + (
C_eff * Tfl[1]) / delts) / (Mfl * Cpwg / Aseg + C_eff / (2 * delts))
Tfl[2] = (TinSub + Tout) / 2
DT[2] = Tfl[2] - Te
qdiff = Mfl / Aseg * Cpwg * 2 * (DT[2] - DT[1])
else:
Tout = Tfl[1] + (qgain * delts) / C_eff
Tfl[2] = Tout
DT[2] = Tfl[2] - Te
qdiff = 5 * (DT[2] - DT[1])
if abs(qdiff < 0.1):
DT[1] = DT[2]
exit = True
else:
if xgain > 40:
DT[1] = (DT[1] + DT[2]) / 2
if xgain == xgainmax:
exit = True
else:
DT[1] = DT[2]
qout = Mfl * Cpwg * (Tout - TinSub) / Aseg
qmtherm = (Tfl[2] - Tfl[1]) * C_eff / delts
qbal = qgain - qout - qmtherm
if abs(qbal) > 1:
qbal = qbal
return qgain
[docs]def Calc_qloss_net(Mfl, Le, Area_a, Tm, Te, maxmsc):
qloss = 0.217 * Le * Area_a * (Tm - Te) * (Mfl / maxmsc) / 1000
return qloss # in kW
[docs]def calc_anglemodifierSC(Az_vector, g_vector, ha_vector, teta_z, tilt_angle, type_SCpanel, latitude, Sz_vector):
def calc_Teta_L(Az, teta_z, tilt, Sz):
# calculate incident angles longitudinal and trasnversally of the solar collector
teta_la = tan(Sz) * cos(teta_z - Az)
teta_l = degrees(abs(atan(teta_la) - tilt)) # longitudinal incidence angle
if teta_l < 0:
teta_l = min(89, abs(teta_T))
if teta_l >= 90:
teta_l = 89.999
return teta_l # in degrees
def calc_Teta_T(Az, Sz, teta_z): # Az is the solar azimuth
# calculate incident angles transversal
teta_ta = sin(Sz) * sin(abs(teta_z - Az))
teta_T = degrees(atan(teta_ta / cos(teta_ta))) # transversal angle modifier
if teta_T < 0:
teta_T = min(89, abs(teta_T))
if teta_T >= 90:
teta_T = 89.999
return teta_T
def calc_maxtetaL(teta_L):
if teta_L < 0:
teta_L = min(89, abs(teta_L))
if teta_L >= 90:
teta_L = 89.999
return teta_L
def Calc_IAMb(teta_l, teta_T, type_SCpanel):
if type_SCpanel == 1: # # Flat plate collector SOLEX blu SFP, 2012
IAM_b = -0.00000002127039627042 * teta_l ** 4 + 0.00000143550893550934 * teta_l ** 3 - 0.00008493589743580050 * teta_l ** 2 + 0.00041588966590833100 * teta_l + 0.99930069929920900000
if type_SCpanel == 2: # # evacuated tube Zewotherm SOL ZX-30 SFP, 2012
IAML = -0.00000003365384615386 * teta_l ** 4 + 0.00000268745143745027 * teta_l ** 3 - 0.00010196678321666700 * teta_l ** 2 + 0.00088830613832779900 * teta_l + 0.99793706293541500000
IAMT = 0.000000002794872 * teta_T ** 5 - 0.000000534731935 * teta_T ** 4 + 0.000027381118880 * teta_T ** 3 - 0.000326340326281 * teta_T ** 2 + 0.002973799531468 * teta_T + 1.000713286764210
IAM_b = IAMT * IAML # overall incidence angle modifier for beam radiation
return IAM_b
# convert to radians
teta_z = radians(teta_z)
tilt = radians(tilt_angle)
# Az_vector = np.radians(Az_vector)
g_vector = np.radians(g_vector)
ha_vector = np.radians(ha_vector)
lat = radians(latitude)
Sz_vector = np.radians(Sz_vector)
Az_vector = np.radians(Az_vector)
Incidence_vector = np.vectorize(Calc_incidenteangleB)(g_vector, lat, ha_vector, tilt,
teta_z) # incident angle in radians
# calculate incident angles
if type_SCpanel == 1: #
Teta_L = np.degrees(Incidence_vector)
Teta_T = 0 # not necessary
Teta_L = np.vectorize(calc_maxtetaL)(Teta_L)
if type_SCpanel == 2: #
Teta_L = np.vectorize(calc_Teta_L)(Az_vector, teta_z, tilt, Sz_vector) # in degrees
Teta_T = np.vectorize(calc_Teta_T)(Az_vector, Sz_vector, teta_z, Incidence_vector) # in degrees
# calculate incident angle modifier
IAM_b_vector = np.vectorize(Calc_IAMb)(Teta_L, Teta_T, type_SCpanel)
return IAM_b_vector
[docs]def Calc_incidenteangleB(g, lat, ha, tilt, teta_z):
# calculate incident angle beam radiation
part1 = sin(lat) * sin(g) * cos(tilt) - cos(lat) * sin(g) * sin(tilt) * cos(teta_z)
part2 = cos(lat) * cos(g) * cos(ha) * cos(tilt) + sin(lat) * cos(g) * cos(ha) * sin(tilt) * cos(teta_z)
part3 = cos(g) * sin(ha) * sin(tilt) * sin(teta_z)
teta_B = acos(part1 + part2 + part3)
return teta_B # in radains
[docs]def calc_properties_SC(type_SCpanel):
"""
properties of module
:param type_SCpanel:
:return:
"""
if type_SCpanel == 1: # # Flat plate collector SOLEX blu SFP, 2012
n0 = 0.775 # zero loss efficeincy
c1 = 3.91 # W/m2K
c2 = 0.0081 # W/m2K2
# specific mass flow rates
mB0_r = 57.98 # # in kg/h/m2 of aperture area
mB_max_r = 86.97 # in kg/h/m2 of aperture area
mB_min_r = 28.99 # in kg/h/m2 of aperture area
C_eff = 8000 # thermal capacitance of module J/m2K
t_max = 192 # stagnation temperature in C
IAM_d = 0.87 # diffuse incident angle considered at 50 degrees
Aratio = 0.888 # the aperture/gross area ratio
Apanel = 2.023 # m2
dP1 = 0
dP2 = 170 / (Aratio * Apanel)
dP3 = 270 / (Aratio * Apanel)
dP4 = 80 / (Aratio * Apanel)
if type_SCpanel == 2: # # evacuated tube Zewotherm SOL ZX-30 SFP, 2012
n0 = 0.721
c1 = 0.89 # W/m2K
c2 = 0.0199 # W/m2K2
# specific mass flow rates
mB0_r = 88.2 # in kg/h/m2 of aperture area
mB_max_r = 147.12 # in kg/h/m2
mB_min_r = 33.10 # in kg/h/m2
C_eff = 38000 # thermal capacitance of module anf fluid for Brine J/m2K
t_max = 196 # stagnation temperature in C
IAM_d = 0.91 # diffuse incident angle considered at 50 degrees
Aratio = 0.655 # the aperture/gross area ratio
Apanel = 4.322 # m2
dP1 = 0 # in Pa per m2
dP2 = 8000 / (Aratio * Apanel) # in Pa per m2
dP3 = 22000 / (Aratio * Apanel) # in Pa per m2
dP4 = 2000 / (Aratio * Apanel) # in Pa per m2
# Fluids Cp [kJ/kgK] Density [kg/m3] Used for
# Water-glyucol 33% 3.68 1044 Collector Loop
# Water 4.19 998 Secondary collector loop, store, loops for auxiliary
return n0, c1, c2, mB0_r, mB_max_r, mB_min_r, C_eff, t_max, IAM_d, Aratio, Apanel, dP1, dP2, dP3, dP4
[docs]def calc_Eaux_SC(qV_des, Dp_collector, Leq, Aa):
"""
auxiliary electricity solar collectort
:param qV_des:
:param Dp_collector:
:param Leq:
:param Aa:
:return:
"""
Ro = 1000 # kg/m3
dpl = 200 # pressure losses per length of pipe according to solar districts
Fcr = 1.3 # factor losses of accessories
Dp_friction = dpl * Leq * Aa * Fcr # HANZENWILIAMSN PA/M
Eaux = (qV_des / Ro) * (Dp_collector + Dp_friction) / 0.6 / 10 # energy necesary in kW from pumps
return Eaux # energy spent in kWh
# minimization of Eaux
[docs]def SelectminimumenergySc(q1, q2, q3, q4, E1, E2, E3, E4, m1, m2, m3, m4, dP1, dP2, dP3, dP4, Area_a):
mopt = np.empty(8760)
dpopt = np.empty(8760)
const = Area_a / 3600
ms = [m1 * const, m2 * const, m3 * const, m4 * const] # float points
dps = [dP1 * Area_a, dP2 * Area_a, dP3 * Area_a, dP4 * Area_a] # float points
balances = [q1 - E1 * 2, q2 - E2 * 2, q3 - E3 * 2,
q4 - E4 * 2] # vectors #quality of electricity is twice as expensive
for time in range(8760):
balances2 = [balances[0][time], balances[1][time], balances[2][time], balances[3][time]]
maxenergy = np.max(balances2)
ix_maxenergy = np.where(balances2 == maxenergy)
mopt[time] = ms[ix_maxenergy[0][0]]
dpopt[time] = dps[ix_maxenergy[0][0]]
return mopt, dpopt
[docs]def Selectminimumenergy2(m, q, dp):
for time in range(8760):
if q[time] <= 0:
m[time] = 0
dp[time] = 0
return m, dp
# optimal angle and tilt
[docs]def optimal_angle_and_tilt(observers_all, latitude, worst_sh, worst_Az, transmittivity,
grid_side, module_lenght, angle_north, Min_Isol, Max_Isol):
def Calc_optimal_angle(teta_z, latitude, transmissivity):
if transmissivity <= 0.15:
gKt = 0.977
elif 0.15 < transmissivity <= 0.7:
gKt = 1.237 - 1.361 * transmissivity
else:
gKt = 0.273
Tad = 0.98
Tar = 0.97
Pg = 0.2 # ground reflectance of 0.2
l = radians(latitude)
a = radians(teta_z) # this is surface azimuth
b = atan((cos(a) * tan(l)) * (1 / (1 + ((Tad * gKt - Tar * Pg) / (2 * (1 - gKt))))))
return degrees(b)
def Calc_optimal_spacing(Sh, Az, tilt_angle, module_lenght):
h = module_lenght * sin(radians(tilt_angle))
D1 = h / tan(radians(Sh))
D = max(D1 * cos(radians(180 - Az)), D1 * cos(radians(Az - 180)))
return D
def Calc_categoriesroof(teta_z, B, GB, Max_Isol):
if -122.5 < teta_z <= -67:
CATteta_z = 1
elif -67 < teta_z <= -22.5:
CATteta_z = 3
elif -22.5 < teta_z <= 22.5:
CATteta_z = 5
elif 22.5 < teta_z <= 67:
CATteta_z = 4
elif 67 <= teta_z <= 122.5:
CATteta_z = 2
if 0 < B <= 5:
CATB = 1 # flat roof
elif 5 < B <= 15:
CATB = 2 # tilted 25 degrees
elif 15 < B <= 25:
CATB = 3 # tilted 25 degrees
elif 25 < B <= 40:
CATB = 4 # tilted 25 degrees
elif 40 < B <= 60:
CATB = 5 # tilted 25 degrees
elif B > 60:
CATB = 6 # tilted 25 degrees
GB_percent = GB / Max_Isol
if 0 < GB_percent <= 0.25:
CATGB = 1 # flat roof
elif 0.25 < GB_percent <= 0.50:
CATGB = 2
elif 0.50 < GB_percent <= 0.75:
CATGB = 3
elif 0.75 < GB_percent <= 0.90:
CATGB = 4
elif 0.90 < GB_percent <= 1:
CATGB = 5
return CATB, CATGB, CATteta_z
# calculate values for flat roofs Slope < 5 degrees.
optimal_angle_flat = Calc_optimal_angle(0, latitude, transmittivity)
optimal_spacing_flat = Calc_optimal_spacing(worst_sh, worst_Az, optimal_angle_flat, module_lenght)
arcpy.AddField_management(observers_all, "array_s", "DOUBLE")
arcpy.AddField_management(observers_all, "area_netpv", "DOUBLE")
arcpy.AddField_management(observers_all, "CATB", "SHORT")
arcpy.AddField_management(observers_all, "CATGB", "SHORT")
arcpy.AddField_management(observers_all, "CATteta_z", "SHORT")
fields = ('aspect', 'slope', 'GB', "array_s", "area_netpv", "CATB", "CATGB", "CATteta_z")
# go inside the database and perform the changes
with arcpy.da.UpdateCursor(observers_all, fields) as cursor:
for row in cursor:
aspect = row[0]
slope = row[1]
if slope > 5: # no t a flat roof.
B = slope
array_s = 0
if 180 <= aspect < 360: # convert the aspect of arcgis to azimuth
teta_z = aspect - 180
elif 0 < aspect < 180:
teta_z = aspect - 180 # negative in the east band
elif aspect == 0 or aspect == 360:
teta_z = 180
if -angle_north <= teta_z <= angle_north and row[2] > Min_Isol:
row[0] = teta_z
row[1] = B
row[3] = array_s
row[4] = (grid_side - array_s) / cos(radians(abs(B))) * grid_side
row[5], row[6], row[7] = Calc_categoriesroof(teta_z, B, row[2], Max_Isol)
cursor.updateRow(row)
else:
cursor.deleteRow()
else:
teta_z = 0 # flat surface, all panels will be oriented towards south # optimal angle in degrees
B = optimal_angle_flat
array_s = optimal_spacing_flat
if row[2] > Min_Isol:
row[0] = teta_z
row[1] = B
row[3] = array_s
row[4] = (grid_side - array_s) / cos(radians(abs(B))) * grid_side
row[5], row[6], row[7] = Calc_categoriesroof(teta_z, B, row[2], Max_Isol)
cursor.updateRow(row)
else:
cursor.deleteRow()
# investment and maintenance costs
[docs]def calc_Cinv_SC(Area):
"""
Lifetime 35 years
"""
InvCa = 2050 * Area / 35 # [CHF/y]
return InvCa
[docs]def test_solar_collector():
import cea.inputlocator
locator = cea.inputlocator.InputLocator(r'C:\reference-case\baseline')
# for the interface, the user should pick a file out of of those in ...DB/Weather/...
weather_path = locator.get_default_weather()
radiation = locator.get_radiation()
gv = cea.globalvar.GlobalVariables()
calc_SC(locator=locator, radiation=radiation, latitude=46.95240555555556, longitude=7.439583333333333, year=2014,
gv=gv,
weather_path=weather_path)
if __name__ == '__main__':
test_solar_collector()