WIP

Cargo.toml

@@ -11,6 +11,7 @@ arrayvec = { version = "0.7.4", features = ["serde"] }
 clap = { version = "4.5.4", features = ["derive"] }
 csv = "1.3.0"
 derivative = "2.2.0"
+fast_ode = "1.0.0"
 indexmap = { version = "2.2.6", features = ["serde"] }
 indicatif = "0.17.8"
 interp = "1.0.3"
@@ -25,6 +26,7 @@ polyfit-rs = "0.2.1"
 serde = { version = "1.0.197", features = ["derive", "rc"] }
 serde-enum-str = "0.4.0"
 serde_json = "1.0.115"
+serde_repr = "0.1.18"
 variants-struct = "0.1.1"
 
 [dev-dependencies]

src/core/heating_systems/common.rs

@@ -2,9 +2,11 @@ use crate::core::heating_systems::boiler::{
     BoilerServiceSpace, BoilerServiceWaterCombi, BoilerServiceWaterRegular,
 };
 use crate::core::heating_systems::heat_battery::HeatBatteryServiceWaterRegular;
-use crate::core::heating_systems::heat_network::HeatNetworkServiceWaterStorage;
+use crate::core::heating_systems::heat_network::{
+    HeatNetworkServiceSpace, HeatNetworkServiceWaterStorage,
+};
 use crate::core::heating_systems::heat_pump::{
-    HeatPumpHotWaterOnly, HeatPumpServiceSpaceWarmAir, HeatPumpServiceWater,
+    HeatPumpHotWaterOnly, HeatPumpServiceSpace, HeatPumpServiceSpaceWarmAir, HeatPumpServiceWater,
 };
 use crate::core::heating_systems::instant_elec_heater::InstantElecHeater;
 use crate::simulation_time::SimulationTimeIteration;
@@ -86,3 +88,25 @@ impl SpaceHeatSystem {
         }
     }
 }
+
+pub enum SpaceHeatingService {
+    HeatPump(HeatPumpServiceSpace),
+    Boiler(BoilerServiceSpace),
+    HeatNetwork(HeatNetworkServiceSpace),
+    HeatBattery(()),
+}
+
+// macro so accessing individual controls through the enum isn't so repetitive
+#[macro_use]
+macro_rules! per_space_heating {
+    ($val:expr, $pattern:pat => { $res:expr }) => {
+        match $val {
+            SpaceHeatingService::HeatPump($pattern) => $res,
+            SpaceHeatingService::Boiler($pattern) => $res,
+            SpaceHeatingService::HeatNetwork($pattern) => $res,
+            SpaceHeatingService::HeatBattery($pattern) => unreachable!(),
+        }
+    };
+}
+
+pub(crate) use per_space_heating;

src/core/heating_systems/emitters.rs

@@ -0,0 +1,237 @@
+use crate::compare_floats::max_of_2;
+use crate::core::heating_systems::common::SpaceHeatingService;
+use crate::core::space_heat_demand::zone::Zone;
+use crate::external_conditions::ExternalConditions;
+use crate::input::{EcoDesignController, EcoDesignControllerClass};
+use crate::simulation_time::SimulationTimeIteration;
+use bacon_sci::ivp::{RungeKuttaSolver, RK45};
+use std::sync::Arc;
+
+/// This module provides objects to represent radiator and underfloor emitter systems.
+
+pub struct Emitters {
+    thermal_mass: f64,
+    c: f64,
+    n: f64,
+    temp_diff_emit_dsgn: f64,
+    frac_convective: f64,
+    heat_source: Arc<SpaceHeatingService>,
+    zone: Arc<Zone>,
+    external_conditions: Arc<ExternalConditions>,
+    design_flow_temp: f64,
+    ecodesign_controller_class: EcoDesignControllerClass,
+    min_outdoor_temp: Option<f64>,
+    max_outdoor_temp: Option<f64>,
+    min_flow_temp: Option<f64>,
+    max_flow_temp: Option<f64>,
+    simulation_timestep: f64,
+    temp_emitter_prev: f64,
+}
+
+impl Emitters {
+    /// Construct an Emitters object
+    ///
+    /// Arguments:
+    /// * `thermal_mass` - thermal mass of emitters, in kWh / K
+    /// * `c` - constant from characteristic equation of emitters (e.g. derived from BS EN 442 tests)
+    /// * `n` - exponent from characteristic equation of emitters (e.g. derived from BS EN 442 tests)
+    /// * `temp_diff_emit_dsgn` - design temperature difference across the emitters, in deg C or K
+    /// * `frac_convective` - convective fraction for heating
+    /// * `heat_source` - reference to an object representing the system (e.g.
+    ///                       boiler or heat pump) providing heat to the emitters
+    /// * `zone` - reference to the Zone object representing the zone in which the
+    ///             emitters are located
+    /// * `simulation_timestep` - timestep length for simulation time being used in this context
+    ///
+    /// Other variables:
+    /// * `temp_emitter_prev` - temperature of the emitters at the end of the
+    ///    previous timestep, in deg C
+    pub fn new(
+        thermal_mass: f64,
+        c: f64,
+        n: f64,
+        temp_diff_emit_dsgn: f64,
+        frac_convective: f64,
+        heat_source: Arc<SpaceHeatingService>,
+        zone: Arc<Zone>,
+        external_conditions: Arc<ExternalConditions>,
+        ecodesign_controller: EcoDesignController,
+        design_flow_temp: f64,
+        simulation_timestep: f64,
+    ) -> Self {
+        let ecodesign_controller_class = ecodesign_controller.ecodesign_control_class;
+        let (min_outdoor_temp, max_outdoor_temp, min_flow_temp, max_flow_temp) = if matches!(
+            ecodesign_controller_class,
+            EcoDesignControllerClass::Class2
+                | EcoDesignControllerClass::Class3
+                | EcoDesignControllerClass::Class6
+                | EcoDesignControllerClass::Class7
+        ) {
+            (
+                ecodesign_controller.min_outdoor_temp,
+                ecodesign_controller.max_outdoor_temp,
+                ecodesign_controller.min_flow_temp,
+                Some(design_flow_temp),
+            )
+        } else {
+            (None, None, None, None)
+        };
+        Self {
+            thermal_mass,
+            c,
+            n,
+            temp_diff_emit_dsgn,
+            frac_convective,
+            heat_source,
+            zone,
+            external_conditions,
+            design_flow_temp,
+            ecodesign_controller_class,
+            min_outdoor_temp,
+            max_outdoor_temp,
+            min_flow_temp,
+            max_flow_temp,
+            simulation_timestep,
+            temp_emitter_prev: 20.0,
+        }
+    }
+
+    pub fn temp_setpnt(&self, simulation_time_iteration: &SimulationTimeIteration) -> Option<f64> {
+        match self.heat_source.as_ref() {
+            SpaceHeatingService::HeatPump(heat_pump) => {
+                heat_pump.temp_setpnt(simulation_time_iteration)
+            }
+            SpaceHeatingService::Boiler(boiler) => Some(boiler.temp_setpnt()),
+            SpaceHeatingService::HeatNetwork(heat_network) => {
+                heat_network.temperature_setpnt(simulation_time_iteration)
+            }
+            SpaceHeatingService::HeatBattery(_) => unreachable!(),
+        }
+    }
+
+    pub fn in_required_period(
+        &self,
+        simulation_time_iteration: &SimulationTimeIteration,
+    ) -> Option<bool> {
+        match self.heat_source.as_ref() {
+            SpaceHeatingService::HeatPump(heat_pump) => {
+                heat_pump.in_required_period(simulation_time_iteration)
+            }
+            SpaceHeatingService::Boiler(boiler) => Some(boiler.in_required_period()),
+            SpaceHeatingService::HeatNetwork(heat_network) => {
+                heat_network.in_required_period(simulation_time_iteration)
+            }
+            SpaceHeatingService::HeatBattery(_) => unreachable!(),
+        }
+    }
+
+    pub fn frac_convective(&self) -> f64 {
+        self.frac_convective
+    }
+
+    pub fn temp_flow_return(
+        &self,
+        simulation_time_iteration: &SimulationTimeIteration,
+    ) -> (f64, f64) {
+        let flow_temp = match self.ecodesign_controller_class {
+            EcoDesignControllerClass::Class2
+            | EcoDesignControllerClass::Class3
+            | EcoDesignControllerClass::Class6
+            | EcoDesignControllerClass::Class7 => {
+                // A heater flow temperature control that varies the flow temperature of
+                // water leaving the heat dependant upon prevailing outside temperature
+                // and selected weather compensation curve.
+                //
+                // They feature provision for manual adjustment of the weather
+                // compensation curves and therby introduce a technical risk that optimal
+                // minimised flow temperatures are not always achieved.
+
+                // use weather temperature at the timestep
+                let outside_temp = self.external_conditions.air_temp(simulation_time_iteration);
+
+                let min_flow_temp = self.min_flow_temp.unwrap();
+                let max_flow_temp = self.max_flow_temp.unwrap();
+                let min_outdoor_temp = self.min_outdoor_temp.unwrap();
+                let max_outdoor_temp = self.max_outdoor_temp.unwrap();
+
+                // set outdoor and flow temp limits for weather compensation curve
+                if outside_temp < min_outdoor_temp {
+                    max_flow_temp
+                } else if outside_temp > max_outdoor_temp {
+                    min_flow_temp
+                } else {
+                    // Interpolate
+                    // Note: A previous version used numpy interpolate, but this
+                    //        seemed to be giving incorrect results, so interpolation
+                    //        is implemented manually here.
+                    min_flow_temp
+                        + (outside_temp - max_outdoor_temp)
+                            * ((max_flow_temp - min_flow_temp)
+                                / (min_outdoor_temp - max_outdoor_temp))
+                }
+            }
+            _ => self.design_flow_temp,
+        };
+
+        let return_temp = if flow_temp >= 70.0 {
+            60.0
+        } else {
+            flow_temp * 6.0 / 7.0
+        };
+
+        (flow_temp, return_temp)
+    }
+
+    /// Calculate emitter output at given emitter and room temp
+    ///
+    /// Power output from emitter (eqn from 2020 ASHRAE Handbook p644):
+    ///            power_output = c * (T_E - T_rm) ^ n
+    ///        where:
+    ///            T_E is mean emitter temperature
+    ///            T_rm is air temperature in the room/zone
+    ///            c and n are characteristic of the emitters (e.g. derived from BS EN 442 tests)
+    pub fn power_output_emitter(&self, temp_emitter: f64, temp_rm: f64) -> f64 {
+        self.c * max_of_2(0., temp_emitter - temp_rm).powf(self.n)
+    }
+
+    /// Calculate emitter temperature that gives required power output at given room temp
+    ///
+    /// Power output from emitter (eqn from 2020 ASHRAE Handbook p644):
+    ///            power_output = c * (T_E - T_rm) ^ n
+    ///        where:
+    ///            T_E is mean emitter temperature
+    ///            T_rm is air temperature in the room/zone
+    ///            c and n are characteristic of the emitters (e.g. derived from BS EN 442 tests)
+    ///        Rearrange to solve for T_E
+    pub fn temp_emitter_req(&self, power_emitter_req: f64, temp_rm: f64) -> f64 {
+        (power_emitter_req / self.c).powf(1. / self.n) + temp_rm
+    }
+
+    fn func_temp_emitter_change_rate(
+        &self,
+        power_input: f64,
+    ) -> impl FnOnce(f64, &[f64]) -> f64 + '_ {
+        let Self {
+            c, n, thermal_mass, ..
+        } = self;
+
+        move |t, temp_diff: &[f64]| {
+            (power_input - c * max_of_2(0., temp_diff[0]).powf(*n)) / thermal_mass
+        }
+    }
+
+    /// Calculate emitter temperature after specified time with specified power input
+    // pub fn temp_emitter(
+    //     &self,
+    //     time_start: f64,
+    //     time_end: f64,
+    //     temp_emitter_start: f64,
+    //     temp_rm: f64,
+    //     power_input: f64,
+    //     temp_emitter_max: Option<f64>,
+    // ) -> Result<(f64, f64), &'static str> {
+    //     // Calculate emitter temp at start of timestep
+    //
+    //     let temp_diff_start = temp_emitter_start - temp_rm;
+    // }
+}

src/core/heating_systems/mod.rs

@@ -5,6 +5,7 @@ pub mod heat_network;
 pub mod point_of_use;
 
 pub mod common;
+pub mod emitters;
 pub mod heat_battery;
 pub mod heat_pump;
 pub mod instant_elec_heater;