Merge pull request #4 from mcneillj/surface-radiation

Merge surface radiation calculations

docs/api.rst

@@ -13,6 +13,15 @@ Settings
 .. automodule:: sitka.general.settings
    :members:
 
+Utilities
+=======
+
+Time Series
+~~~~~~~~
+
+.. automodule:: sitka.utils.time_series
+   :members:
+
 Input/Output
 ============
 
@@ -37,11 +46,43 @@ Solar
 .. automodule:: sitka.calculations.solar
    :members:
 
+
+Conduction
+~~~~~~
+
+.. automodule:: sitka.calculations.conduction
+   :members:
+
+
+Convection
+~~~~~~
+
+.. automodule:: sitka.calculations.solar.convection
+   :members:
+
+Fluid
+~~~~~~
+
+.. automodule:: sitka.calculations.fluids
+   :members:
+
+Radiation
+~~~~~~
+
+.. automodule:: sitka.calculations.radiation
+   :members:
+
 Components
 ==========
 
 Site
 ~~~~~~~~
 
 .. automodule:: sitka.components.site
+   :members:
+
+Surface
+~~~~~~~~
+
+.. automodule:: sitka.components.surface
    :members:

setup.py

@@ -4,7 +4,7 @@
 
 setup(
     name='sitka',
-    version='0.1.3',
+    version='0.1.4',
     description='A built environment analysis and modeling library for Python.',
     long_description="""
     Sitka is a built environment analysis and modeling library for Python.

sitka/calculations/conduction.py

@@ -0,0 +1,99 @@
+import numpy as np
+import pandas as pd
+
+class FiniteDifferenceMethod1D:
+    def __init__(self, time, surface, zone_air):
+        # Explicit solver
+
+        # General properties
+        self.time = time
+        self.surface = surface
+        self.zone_air = zone_air
+
+        # Stored variables
+        self.temperature_array = None
+        self.time_array = None
+
+        # Number of internal points
+        self.Nx = len(self.surface.layers)
+
+        # Calculate Spatial Step-Size
+        self.thickness = self.surface.thickness
+        self.dx = self.thickness/self.Nx
+        kx = self.dx/2
+
+        # Create grid-points on x axis
+        self.x = np.linspace(0,1,self.Nx+2)
+        self.x = self.x[1:-1]
+
+        # FD matrix
+        self.A = None
+        self.b = None
+        self.u = None
+
+        # Initial methods
+        self.update_calculated_values()
+
+    def update_calculated_values(self):
+        self.setup_finite_difference_matrix()
+        self.initialize_arrays()
+
+    def setup_finite_difference_matrix(self):
+        R = np.array(self.surface.thermal_resistance_array)
+        C = np.array(self.surface.thermal_capacitance_array)
+        dt = self.time.time_step
+        Nx = self.Nx
+
+        A1 = dt/(R[:-1]*C)
+        A2 = 1-(dt/C)*(1/R[1:]+1/R[:-1])
+        A3 = dt/(R[1:]*C)
+
+        A = np.zeros((Nx,Nx+2))
+        A[:,:-2] += np.diag(A1)
+        A[:,1:-1] += np.diag(A2)
+        A[:,2:] += np.diag(A3)
+
+        self.A = A
+
+    def initialize_arrays(self):
+        #self.temperature_array = np.zeros((len(self.time.time_range),len(self.x)))
+        #self.time_array = np.zeros(len(self.time.time_range))
+
+        df = pd.DataFrame()
+        for x in self.x:
+            df[str(x)] = np.ones(self.time.length)*self.zone_air.initial_zone_air_temperature
+        #df.index = self.time.time_range
+        self.temperature_array = df
+        self.time_array = pd.Series(self.time.time_range)
+
+    def run_solver(self, iteration, time_index):
+        surface = self.surface
+        #t = self.time.time_range[start_time_step:end_time_step]
+        x = self.x
+
+        #for sim_index, time_index in enumerate(self.time_array[start_time_step:end_time_step].index):
+        #array_index = start_time_step + sim_time_step
+        if time_index > 0 and iteration == 0:
+            b = np.array(self.temperature_array.loc[time_index-1])
+        else:
+            b = np.array(self.temperature_array.loc[time_index])
+
+        # Boundary conditions
+        inside_temperature = self.zone_air.zone_air_temperature.iloc[time_index]
+        outside_temperature = self.surface.weather.dry_bulb_temperature.iloc[time_index]
+
+        # Solve timestep
+        self.solve_timestep(time_index, b, outside_temperature, inside_temperature)
+
+    def solve_timestep(self, time_index, temperature_array, outside_temperature, inside_temperature):
+        # Boundary Conditions
+        T = outside_temperature
+        b = np.insert(temperature_array,0,T,axis=0)
+        b = np.append(b, inside_temperature)
+
+        # Solve
+        u = np.dot(self.A,b)
+        b = np.array(u)
+
+        # Store values in object
+        self.temperature_array.iloc[time_index] = u

sitka/calculations/convection.py

@@ -0,0 +1,146 @@
+"""Inside and outside surface convection models.
+"""
+import numpy as np
+import pandas as pd
+
+from sitka.utils.time_series import TimeSeriesComponent
+
+
+class OutsideConvection(TimeSeriesComponent):
+    v
+    def __init__(self, time, weather, surface, surface_temperature):
+        # Model parameters
+        self.heat_transfer_coefficient = None  # Convection coefficient (default 22.7) [W/m^2-K]
+
+        # Geometric parameters
+        self.surface_height = None
+
+        # Thermal properties
+        self.wind_speed = None
+        self.surface_temperature = surface_temperature
+        self.heat_transfer_rate = None
+
+        # Associated objects
+        self._weather = weather
+        self._surface = surface
+
+        # Add attributes from super class
+        super().__init__(time)
+
+        # Run method to update all calculated values
+        self.update_calculated_values()
+
+    def update_calculated_values(self):
+        print('Updating external longwave radiation calculations.')
+        self.calculate_surface_wind_speed()
+        self.calculate_heat_transfer_coefficient()
+        self.calculate_heat_transfer()
+
+    def calculate_surface_wind_speed(self):
+        """
+        Calculate the wind speed at the surface based on the height.
+
+        Yields
+        ----------
+        wind_speed : Series
+
+        References
+        --------
+        """
+        self.wind_speed = self.weather.wind_speed*((270/10)^0.14*(0.5*self.surface_height/370)^0.22)
+
+    def calculate_heat_transfer_coefficient(self):
+        """
+        Calculate the outside convection coefficient on the surface.
+
+        Yields
+        ----------
+        heat_transfer_coefficient : Series
+
+        References
+        --------
+        """
+        # Outside convection coefficient
+        D = 10.79
+        E = 4.192
+        F = 0.0
+        self.heat_transfer_coefficient = D+E*self.wind_speed+F*self.wind_speed^2
+
+    def calculate_heat_transfer(self):
+        """
+        Calculate the heat transfer rate for the surface.
+
+        Yields
+        ----------
+        heat_transfer_rate : Series
+
+        References
+        --------
+        """
+        # TODO add heat transfer rate calculation
+        self.heat_transfer_rate = 1
+
+
+class InsideConvection(TimeSeriesComponent):
+    """
+    Inside convection calculation for time-series.
+
+    Parameters
+    ----------
+    time : Time
+    weather : Weather
+    surface : Surface
+    surface_temperature : Series
+
+    Attributes
+    ----------
+    heat_transfer_coefficient : Series
+        Convection heat transfer coefficient [W/m^2-K].
+    surface_height : float
+        Surface height [m].
+    wind_speed : Series
+        Ambient wind speed [m/s].
+    surface_temperature : Series
+        Surface temperature [C].
+    heat_transfer_rate : Series
+        Heat transfer rate [W].
+    time : Time
+    weather : Weather
+    surface : Surface
+    """
+    def __init__(self):
+        self.heat_transfer_coefficient = 8.29 # Convection coefficient [W/m^2-K]
+        self.surface_tilt = 90
+        self.air_temperature = []
+        self.surface_temperature = []
+
+    def ashrae_vertical_wall(self):
+        #  ASHRAE Vertical Wall Model
+        delta_temperature = self.surface_temperature - self.air_temperature
+        h = 1.31*np.abs(delta_temperature)**(1/3)
+        h[h < 0.1] = 0.1
+        return h
+
+    def simple_natural_convection_algorithm(self):
+        # Simple Natural Convection Algorithm [Walton (1983)]
+        #inside_convection_coefficient = 1.31*np.abs(delta_temperature)**(1/3)
+        surface_tilt = self.surface_tilt
+
+        if (surface_tilt == 90):
+            h = 3.076
+            #obj.h_in(tim) = 1.31*np.abs(delta_temperature)**(1/3)
+        elif (surface_tilt  < 90) & (delta_temperature > 0):
+             h = 4.040
+             #obj.h_in(tim) = 9.482*np.abs(delta_temperature)**(1/3)/(7.283-np.abs(np.cos(rad_surface_tilt)))
+        elif (surface_tilt  < 90) & (DT < 0):
+             h = 0.948
+             #obj.h_in(tim) = 1.810*np.abs(delta_temperature)**(1/3)/(1.382+np.abs(np.cos(rad_surface_tilt)))
+        else:
+            h = 3.076
+
+        return h
+
+    def calculate_coefficient(self):
+        # inside surface natural convection coefficient
+        inside_convection_coefficient = simple_natural_convection_algorithm()
+        self.heat_transfer_coefficient = inside_convection_coefficient