From 6032b89189a48ae607f20a9d3ff2ce59dacfac9e Mon Sep 17 00:00:00 2001 From: Douglas Greenshields Date: Wed, 12 Jul 2023 18:40:09 +0100 Subject: [PATCH] =?UTF-8?q?add=20initial=20fast=20solver=20implementation?= =?UTF-8?q?=20=F0=9F=A4=9E?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- Cargo.lock | 142 +++++ Cargo.toml | 1 + .../space_heat_demand/building_element.rs | 105 ++++ src/core/space_heat_demand/mod.rs | 2 + src/core/space_heat_demand/thermal_bridge.rs | 145 +++++ src/core/space_heat_demand/zone.rs | 525 +++++++++++++++--- src/core/units.rs | 5 + src/external_conditions.rs | 2 +- src/input.rs | 4 +- 9 files changed, 863 insertions(+), 68 deletions(-) create mode 100644 src/core/space_heat_demand/building_element.rs create mode 100644 src/core/space_heat_demand/thermal_bridge.rs diff --git a/Cargo.lock b/Cargo.lock index edb365e..79defbb 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -57,6 +57,15 @@ dependencies = [ "windows-sys", ] +[[package]] +name = "approx" +version = "0.5.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "cab112f0a86d568ea0e627cc1d6be74a1e9cd55214684db5561995f6dad897c6" +dependencies = [ + "num-traits", +] + [[package]] name = "autocfg" version = "1.1.0" @@ -69,6 +78,12 @@ version = "2.3.3" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "630be753d4e58660abd17930c71b647fe46c27ea6b63cc59e1e3851406972e42" +[[package]] +name = "bytemuck" +version = "1.13.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "17febce684fd15d89027105661fec94afb475cb995fbc59d2865198446ba2eea" + [[package]] name = "cc" version = "1.0.79" @@ -162,6 +177,7 @@ dependencies = [ "clap", "csv", "itertools", + "nalgebra", "rstest", "serde", "serde_json", @@ -341,18 +357,100 @@ version = "0.4.3" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "09fc20d2ca12cb9f044c93e3bd6d32d523e6e2ec3db4f7b2939cd99026ecd3f0" +[[package]] +name = "matrixmultiply" +version = "0.3.7" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "090126dc04f95dc0d1c1c91f61bdd474b3930ca064c1edc8a849da2c6cbe1e77" +dependencies = [ + "autocfg", + "rawpointer", +] + [[package]] name = "memchr" version = "2.5.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "2dffe52ecf27772e601905b7522cb4ef790d2cc203488bbd0e2fe85fcb74566d" +[[package]] +name = "nalgebra" +version = "0.32.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "307ed9b18cc2423f29e83f84fd23a8e73628727990181f18641a8b5dc2ab1caa" +dependencies = [ + "approx", + "matrixmultiply", + "nalgebra-macros", + "num-complex", + "num-rational", + "num-traits", + "simba", + "typenum", +] + +[[package]] +name = "nalgebra-macros" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "91761aed67d03ad966ef783ae962ef9bbaca728d2dd7ceb7939ec110fffad998" +dependencies = [ + "proc-macro2", + "quote", + "syn 1.0.109", +] + +[[package]] +name = "num-complex" +version = "0.4.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "02e0d21255c828d6f128a1e41534206671e8c3ea0c62f32291e808dc82cff17d" +dependencies = [ + "num-traits", +] + +[[package]] +name = "num-integer" +version = "0.1.45" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "225d3389fb3509a24c93f5c29eb6bde2586b98d9f016636dff58d7c6f7569cd9" +dependencies = [ + "autocfg", + "num-traits", +] + +[[package]] +name = "num-rational" +version = "0.4.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0638a1c9d0a3c0914158145bc76cff373a75a627e6ecbfb71cbe6f453a5a19b0" +dependencies = [ + "autocfg", + "num-integer", + "num-traits", +] + +[[package]] +name = "num-traits" +version = "0.2.15" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "578ede34cf02f8924ab9447f50c28075b4d3e5b269972345e7e0372b38c6cdcd" +dependencies = [ + "autocfg", +] + [[package]] name = "once_cell" version = "1.18.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "dd8b5dd2ae5ed71462c540258bedcb51965123ad7e7ccf4b9a8cafaa4a63576d" +[[package]] +name = "paste" +version = "1.0.13" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b4b27ab7be369122c218afc2079489cdcb4b517c0a3fc386ff11e1fedfcc2b35" + [[package]] name = "pin-project-lite" version = "0.2.10" @@ -407,6 +505,12 @@ dependencies = [ "proc-macro2", ] +[[package]] +name = "rawpointer" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "60a357793950651c4ed0f3f52338f53b2f809f32d83a07f72909fa13e4c6c1e3" + [[package]] name = "rstest" version = "0.17.0" @@ -461,6 +565,15 @@ version = "1.0.14" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "fe232bdf6be8c8de797b22184ee71118d63780ea42ac85b61d1baa6d3b782ae9" +[[package]] +name = "safe_arch" +version = "0.7.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "62a7484307bd40f8f7ccbacccac730108f2cae119a3b11c74485b48aa9ea650f" +dependencies = [ + "bytemuck", +] + [[package]] name = "same-file" version = "1.0.6" @@ -507,6 +620,19 @@ dependencies = [ "serde", ] +[[package]] +name = "simba" +version = "0.8.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "061507c94fc6ab4ba1c9a0305018408e312e17c041eb63bef8aa726fa33aceae" +dependencies = [ + "approx", + "num-complex", + "num-traits", + "paste", + "wide", +] + [[package]] name = "slab" version = "0.4.8" @@ -544,6 +670,12 @@ dependencies = [ "unicode-ident", ] +[[package]] +name = "typenum" +version = "1.16.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "497961ef93d974e23eb6f433eb5fe1b7930b659f06d12dec6fc44a8f554c0bba" + [[package]] name = "unicode-ident" version = "1.0.10" @@ -586,6 +718,16 @@ dependencies = [ "winapi-util", ] +[[package]] +name = "wide" +version = "0.7.11" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "aa469ffa65ef7e0ba0f164183697b89b854253fd31aeb92358b7b6155177d62f" +dependencies = [ + "bytemuck", + "safe_arch", +] + [[package]] name = "winapi" version = "0.3.9" diff --git a/Cargo.toml b/Cargo.toml index 07d88e0..ffcb4b2 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -9,6 +9,7 @@ edition = "2021" clap = { version = "4.3.10", features = ["derive"] } csv = "1.2.2" itertools = "0.11.0" +nalgebra = "0.32.3" rstest = "0.17.0" serde = { version = "1.0.164", features = ["derive"] } serde_json = "1.0.100" diff --git a/src/core/space_heat_demand/building_element.rs b/src/core/space_heat_demand/building_element.rs new file mode 100644 index 0000000..7ad8be1 --- /dev/null +++ b/src/core/space_heat_demand/building_element.rs @@ -0,0 +1,105 @@ +use crate::input::{BuildingElement, MassDistributionClass}; + +pub fn area_for_building_element_input(element: &BuildingElement) -> f64 { + match element { + &BuildingElement::Opaque { area: a, .. } => a, + &BuildingElement::Transparent { area: a, .. } => match a { + Some(a) => a, + None => 0.0, // just to give some nominal value + }, + &BuildingElement::Ground { area: a, .. } => a, + &BuildingElement::AdjacentZTC { area: a, .. } => a, + &BuildingElement::AdjacentZTUSimple { area: a, .. } => a, + } +} + +/// Calculate the number of nodes for a building element +/// (this is a shortcut to properly implementing building elements!) +pub fn number_of_building_element_nodes(element: &BuildingElement) -> u32 { + let k_pli = match element { + &BuildingElement::Opaque { + mass_distribution_class: ref dist_class, + ref k_m, + .. + } => match dist_class { + MassDistributionClass::I => vec![0.0, 0.0, 0.0, 0.0, *k_m], + MassDistributionClass::E => vec![*k_m, 0.0, 0.0, 0.0, 0.0], + MassDistributionClass::IE => { + let k_ie = k_m / 2.0; + vec![k_ie, 0.0, 0.0, 0.0, k_ie] + } + MassDistributionClass::D => { + let k_inner = *k_m / 4.0; + let k_outer = *k_m / 8.0; + vec![k_outer, k_inner, k_inner, k_inner, k_outer] + } + MassDistributionClass::M => vec![0.0, 0.0, *k_m, 0.0, 0.0], + }, + &BuildingElement::AdjacentZTC { + mass_distribution_class: ref dist_class, + ref k_m, + .. + } => match dist_class { + MassDistributionClass::I => vec![0.0, 0.0, 0.0, 0.0, *k_m], + MassDistributionClass::E => vec![*k_m, 0.0, 0.0, 0.0, 0.0], + MassDistributionClass::IE => { + let k_ie = *k_m / 2.0; + vec![k_ie, 0.0, 0.0, 0.0, k_ie] + } + MassDistributionClass::D => { + let k_inner = *k_m / 4.0; + let k_outer = *k_m / 8.0; + vec![k_outer, k_inner, k_inner, k_inner, k_outer] + } + MassDistributionClass::M => vec![0.0, 0.0, *k_m, 0.0, 0.0], + }, + &BuildingElement::AdjacentZTUSimple { + mass_distribution_class: ref dist_class, + ref k_m, + .. + } => match dist_class { + MassDistributionClass::I => vec![0.0, 0.0, 0.0, 0.0, *k_m], + MassDistributionClass::E => vec![*k_m, 0.0, 0.0, 0.0, 0.0], + MassDistributionClass::IE => { + let k_ie = k_m / 2.0; + vec![k_ie, 0.0, 0.0, 0.0, k_ie] + } + MassDistributionClass::D => { + let k_inner = *k_m / 4.0; + let k_outer = *k_m / 8.0; + vec![k_outer, k_inner, k_inner, k_inner, k_outer] + } + MassDistributionClass::M => vec![0.0, 0.0, *k_m, 0.0, 0.0], + }, + &BuildingElement::Ground { + mass_distribution_class: ref dist_class, + ref k_m, + .. + } => { + let k_gr = THICKNESS_GROUND_LAYER * HEAT_CAPACITY_PER_VOLUME_OF_GROUND; + match dist_class { + MassDistributionClass::I => vec![0.0, k_gr, 0.0, 0.0, *k_m], + MassDistributionClass::E => vec![0.0, k_gr, *k_m, 0.0, 0.0], + MassDistributionClass::IE => { + let k_ie = *k_m / 2.0; + vec![0.0, k_gr, k_ie, 0.0, k_ie] + } + MassDistributionClass::D => { + let k_inner = *k_m / 2.0; + let k_outer = *k_m / 4.0; + vec![0.0, k_gr, k_outer, k_inner, k_outer] + } + MassDistributionClass::M => vec![0.0, k_gr, 0.0, *k_m, 0.0], + } + } + &BuildingElement::Transparent { .. } => vec![0.0, 0.0], + }; + + k_pli.len() as u32 +} + +// # Thermal properties of ground from BS EN ISO 13370:2017 Table 7 +// # Use values for clay or silt (same as BR 443 and SAP 10) +const THERMAL_CONDUCTIVITY_OF_GROUND: f64 = 1.5; // in W/(m.K) +const HEAT_CAPACITY_PER_VOLUME_OF_GROUND: f64 = 3000000.0; // in J/(m3.K) +const THICKNESS_GROUND_LAYER: f64 = 0.5; // in m. Specified in BS EN ISO 52016-1:2017 section 6.5.8.2 diff --git a/src/core/space_heat_demand/mod.rs b/src/core/space_heat_demand/mod.rs index e837a7c..720bd5f 100644 --- a/src/core/space_heat_demand/mod.rs +++ b/src/core/space_heat_demand/mod.rs @@ -1,2 +1,4 @@ +mod building_element; +mod thermal_bridge; mod ventilation_element; pub mod zone; diff --git a/src/core/space_heat_demand/thermal_bridge.rs b/src/core/space_heat_demand/thermal_bridge.rs new file mode 100644 index 0000000..96e11d3 --- /dev/null +++ b/src/core/space_heat_demand/thermal_bridge.rs @@ -0,0 +1,145 @@ +use serde_json::Value; +use std::collections::HashMap; + +#[derive(Debug, PartialEq)] +pub enum ThermalBridging { + Number(f64), + Bridges(HashMap), +} + +#[derive(Copy, Clone, Debug, PartialEq)] +pub enum ThermalBridge { + Linear { + linear_thermal_transmittance: f64, + length: f64, + }, + Point { + heat_transfer_coefficient: f64, + }, +} + +pub fn heat_transfer_coefficient_for_thermal_bridge(thermal_bridge: &ThermalBridge) -> f64 { + match thermal_bridge { + &ThermalBridge::Linear { + linear_thermal_transmittance: t, + length: l, + } => t * l, + &ThermalBridge::Point { + heat_transfer_coefficient: h, + } => h, + } +} + +/// Converts the serde_json value for ThermalBridging into structs as have not been able to handle these variants +/// Ultimately this field in the JSON should be divided into different attributes (one for bridging elements, one for a number) +/// and then this would not be necessary (as serde_json would handle it) and could be removed! +pub fn thermal_bridging_from_input(input: Value) -> ThermalBridging { + match input { + Value::Object(map) => { + let mut bridges = HashMap::from([]); + for (name, bridge_object) in map.clone().iter() { + bridges.insert( + (*name).clone(), + match bridge_object.get("type") { + Some(Value::String(s)) if s == &String::from("ThermalBridgeLinear") => { + ThermalBridge::Linear { + linear_thermal_transmittance: bridge_object + .get("linear_thermal_transmittance") + .unwrap() + .as_f64() + .unwrap(), + length: bridge_object.get("length").unwrap().as_f64().unwrap(), + } + } + Some(Value::String(s)) if s == &String::from("ThermalBridgePoint") => { + ThermalBridge::Point { + heat_transfer_coefficient: bridge_object + .get("heat_transfer_coeff") + .unwrap() + .as_f64() + .unwrap(), + } + } + _ => panic!("unknown type of thermal bridging sent"), + }, + ); + } + + ThermalBridging::Bridges(bridges) + } + Value::Number(n) => ThermalBridging::Number(n.as_f64().unwrap()), + _ => panic!("thermal bridging input was not in an expected format"), + } +} + +#[cfg(test)] +mod test { + use super::*; + use rstest::*; + use serde_json::json; + + #[rstest] + pub fn should_convert_serde_json_value_into_thermal_bridging() { + assert_eq!( + thermal_bridging_from_input(json!(2.4)), + ThermalBridging::Number(2.4) + ); + assert_eq!( + thermal_bridging_from_input(json!( + {"TB1": { + "type": "ThermalBridgeLinear", + "linear_thermal_transmittance": 1.0, + "length": 3.0 + }, + "TB2": { + "type": "ThermalBridgeLinear", + "linear_thermal_transmittance": 0.1, + "length": 2.0 + }, + "TB3": { + "type": "ThermalBridgePoint", + "heat_transfer_coeff": 2.0 + }} + )), + ThermalBridging::Bridges(HashMap::from([ + ( + "TB1".to_string(), + ThermalBridge::Linear { + linear_thermal_transmittance: 1.0, + length: 3.0 + } + ), + ( + "TB2".to_string(), + ThermalBridge::Linear { + linear_thermal_transmittance: 0.1, + length: 2.0 + } + ), + ( + "TB3".to_string(), + ThermalBridge::Point { + heat_transfer_coefficient: 2.0 + } + ) + ])) + ); + } + + #[rstest] + pub fn should_have_correct_heat_transfer_coefficient_for_thermal_bridge() { + assert_eq!( + heat_transfer_coefficient_for_thermal_bridge(&ThermalBridge::Linear { + linear_thermal_transmittance: 1.2, + length: 4.0, + }), + 4.8 + ); + assert_eq!( + heat_transfer_coefficient_for_thermal_bridge(&ThermalBridge::Point { + heat_transfer_coefficient: 6.7 + }), + 6.7 + ); + } +} diff --git a/src/core/space_heat_demand/zone.rs b/src/core/space_heat_demand/zone.rs index b67d347..8f08576 100644 --- a/src/core/space_heat_demand/zone.rs +++ b/src/core/space_heat_demand/zone.rs @@ -1,15 +1,31 @@ +use crate::core::space_heat_demand::building_element::{ + area_for_building_element_input, number_of_building_element_nodes, +}; +use crate::core::space_heat_demand::thermal_bridge::{ + heat_transfer_coefficient_for_thermal_bridge, ThermalBridging, +}; use crate::core::space_heat_demand::ventilation_element::{ VentilationElement, WindowOpeningForCooling, }; +use crate::core::units::{kelvin_to_celsius, SECONDS_PER_HOUR}; use crate::input::BuildingElement; +use nalgebra::{vector, DMatrix, DVector, Vector}; use serde::Deserialize; -use serde_json::Value; use std::collections::HashMap; +// Convective fractions +// (default values from BS EN ISO 52016-1:2017, Table B.11) +const F_INT_C: f64 = 0.4; // Can be different for each source of internal gains +const F_SOL_C: f64 = 0.1; + +// Areal thermal capacity of air and furniture +// (default value from BS EN ISO 52016-1:2017, Table B.17) +const K_M_INT: f64 = 10000.0; // J / (m2.K) + pub struct Zone { useful_area: f64, volume: f64, - building_elements: HashMap, + building_elements: Vec, thermal_bridging: ThermalBridging, vent_elements: Vec>, vent_cool_extra: Option, @@ -18,7 +34,7 @@ pub struct Zone { area_el_total: f64, /// internal thermal capacity of the zone, in J / K c_int: f64, - /// dictionary where key is building element and + /// dictionary where key is building element (name) and /// values are 2-element tuples storing matrix row and /// column numbers (both same) where the first element /// of the tuple gives the position of the heat @@ -48,19 +64,19 @@ impl Zone { /// /// ## Arguments /// - /// `area` - useful floor area of the zone, in m2 - /// `volume` - total volume of the zone, m3 - /// `building_elements` - list of BuildingElement objects (walls, floors, windows etc.) - /// `thermal_bridging` - Either: + /// * `area` - useful floor area of the zone, in m2 + /// * `volume` - total volume of the zone, m3 + /// * `building_elements` - list of BuildingElement objects (walls, floors, windows etc.) + /// * `thermal_bridging` - Either: /// - overall heat transfer coefficient for thermal /// bridges in the zone, in W / K /// - list of ThermalBridge objects for this zone - /// `vent_elements` - list of ventilation elements (infiltration, mech vent etc.) - /// `vent_cool_extra` - element providing additional ventilation in response to high + /// * `vent_elements` - list of ventilation elements (infiltration, mech vent etc.) + /// * `vent_cool_extra` - element providing additional ventilation in response to high /// internal temperature - /// `temp_ext_air_init` - external air temperature to use during initialisation, in Celsius - /// `temp_setpnt_init` - setpoint temperature to use during initialisation, in Celsius - pub fn new( + /// * `temp_ext_air_init` - external air temperature to use during initialisation, in Celsius + /// * `temp_setpnt_init` - setpoint temperature to use during initialisation, in Celsius + pub fn new<'a>( area: f64, volume: f64, building_elements: HashMap, @@ -70,78 +86,457 @@ impl Zone { temp_ext_air_init: f64, temp_setpnt_init: f64, ) -> Zone { + let tb_heat_trans_coeff = match thermal_bridging { + ThermalBridging::Number(heat_coeff) => heat_coeff, + ThermalBridging::Bridges(ref bridges) => bridges + .values() + .map(|bridge| heat_transfer_coefficient_for_thermal_bridge(bridge)) + .sum::(), + }; + + let area_el_total = building_elements + .values() + .map(|el| area_for_building_element_input(el)) + .sum::(); + let c_int = K_M_INT * area; + + // Calculate: + // - size of required matrix/vectors (total number of nodes across all + // # building elements + 1 for internal air) + // - positions of heat balance eqns and temperatures in matrix for each node + let mut element_positions = vec![]; + let mut n = 0; + let mut named_building_elements = vec![]; + for (name, building_element) in building_elements.iter() { + let start_idx = n; + n = n + number_of_building_element_nodes(&building_element); + let end_idx = n - 1; + element_positions.push((start_idx as usize, end_idx as usize)); + named_building_elements.push(NamedBuildingElement { + name: name.clone(), + element: building_element.clone(), + }) + } + let zone_idx = n; + let no_of_temps = n + 1; + // sample to start (unimplemented!!) Zone { - useful_area: 0.0, + useful_area: area, volume, - building_elements, + building_elements: named_building_elements, thermal_bridging, vent_elements, vent_cool_extra, - tb_heat_trans_coeff: 0.0, - area_el_total: 0.0, - c_int: 0.0, - element_positions: vec![], - zone_idx: 0, - no_of_temps: 0, + tb_heat_trans_coeff, + area_el_total, + c_int, + element_positions, + zone_idx, + no_of_temps, temp_prev: vec![], } } } -pub enum ThermalBridging { - Number(f64), - Bridges(HashMap), +/// Initialise temperatures of heat balance nodes +// +/// ## Arguments +/// * `temp_ext_air_init` - external air temperature to use during initialisation, in Celsius +/// * `temp_setpnt_init` - setpoint temperature to use during initialisation, in Celsius +/// * `no_of_temps` - number of unknown temperatures (each node in each +/// building element + 1 for internal air) to be +/// solved for +pub fn init_node_temps( + temp_ext_air_init: f64, + temp_setpnt_init: f64, + no_of_temps: u32, +) -> Vec { + // Set starting point for all node temperatures (elements of + // # self.__temp_prev) as average of external air temp and setpoint. This + // is somewhat arbitrary, but of all options for a uniform initial + // temperature, this should lead to relatively fast stabilisation of + // fabric temperatures, which are expected to be close to the external + // air temperature towards the external surface nodes and close to the + // setpoint temperature towards the internal surface nodes, and therefore + // node temperatures on average should be close to the average of the + // external air and setpoint temperatures. + let temp_start = (temp_ext_air_init + temp_setpnt_init) / 2.0; + let mut temp_prev = vec![temp_start; no_of_temps as usize]; + + // Iterate over space heating calculation and meet all space heating + // demand until temperatures stabilise, under steady-state conditions + // using specified constant setpoint and external air temperatures. + loop { + break; + } + + temp_prev } -#[derive(Copy, Clone)] -pub enum ThermalBridge { - Linear { - linear_thermal_transmittance: f64, - length: f64, - }, - Point { - heat_transfer_coefficient: f64, - }, +const DELTA_T_H: u32 = 8760; // hours in a non leap year +const DELTA_T: u32 = DELTA_T_H * SECONDS_PER_HOUR; + +// # Assume default convective fraction for heating/cooling suggested in +// # BS EN ISO 52016-1:2017 Table B.11 +const FRAC_CONVECTIVE: f64 = 0.4; + +/// Calculate heating and cooling demand in the zone for the current timestep +/// +/// According to the procedure in BS EN ISO 52016-1:2017, section 6.5.5.2, steps 1 to 4. +/// +/// # Arguments +/// * `delta_t_h` - calculation timestep, in hours +/// * `temp_ext_air` - temperature of the external air for the current timestep, in deg C +/// * `gains_internal` - internal gains for the current timestep, in W +/// * `gains_solar` - directly transmitted solar gains, in W +/// * `frac_convective_heat` - convective fraction for heating +/// * `frac_convective_cool` - convective fraction for cooling +/// * `temp_setpnt_heat` - temperature setpoint for heating, in deg C +/// * `temp_setpnt_cool` - temperature setpoint for cooling, in deg C +/// * `throughput_factor` - proportional increase in ventilation rate due to +/// over-ventilation requirement +/// * `vent_cool_extra` - (optional) window cooling +pub fn space_heat_cool_demand( + delta_t_h: f64, + temp_ext_air: f64, + gains_internal: f64, + gains_solar: f64, + frac_convective_heat: f64, + frac_convective_cool: f64, + temp_setpnt_heat: f64, + temp_setpnt_cool: f64, + throughput_factor: f64, + vent_cool_extra: Option, +) -> (f64, f64, f64) { + assert!( + temp_setpnt_cool >= temp_setpnt_heat, + "ERROR: Cooling setpoint is below heating setpoint." + ); + + let mut temp_setpnt_cool_vent: Option = None; + if let Some(window_cooling) = vent_cool_extra { + // TODO: get temp_setpnt from window opening control if it exists + // (there's some unimplemented logic here in the Python 🚨) + // + // Set cooling setpoint to Planck temperature to ensure no cooling demand + let temp_setpnt_cool_vent = Some(kelvin_to_celsius(1.4e32)); + assert!( + temp_setpnt_cool_vent.is_none() || (temp_setpnt_cool_vent.unwrap() >= temp_setpnt_heat), + "ERROR: Setpoint for additional ventilation is below heating setpoint." + ); + } + + // Calculate timestep in seconds + let delta_t = delta_t_h * SECONDS_PER_HOUR as f64; + + // For calculation of demand, set heating/cooling gains to zero + let gains_heat_cool = 0.0; + + // Calculate node and internal air temperatures with heating/cooling gains of zero + + (0.0, 0.0, 0.0) } -pub fn thermal_bridging_from_input<'a>(input: Value) -> ThermalBridging { - match input { - Value::Object(map) => { - let mut bridges = HashMap::from([]); - for (name, bridge_object) in map.clone().iter() { - bridges.insert( - (*name).clone(), - match bridge_object.get("type") { - Some(Value::String(s)) if s == &String::from("ThermalBridgeLinear") => { - ThermalBridge::Linear { - linear_thermal_transmittance: bridge_object - .get("linear_thermal_transmittance") - .unwrap() - .as_f64() - .unwrap(), - length: bridge_object.get("length").unwrap().as_f64().unwrap(), - } - } - Some(Value::String(s)) if s == &String::from("ThermalBridgePoint") => { - ThermalBridge::Point { - heat_transfer_coefficient: bridge_object - .get("heat_transfer_coeff") - .unwrap() - .as_f64() - .unwrap(), - } - } - _ => panic!("unknown type of thermal bridging sent"), - }, - ); +/// Calculate temperatures according to procedure in BS EN ISO 52016-1:2017, section 6.5.6 +// +/// ## Arguments +/// * `delta_t` - calculation timestep, in seconds +/// * `temp_prev` - temperature vector X (see below) from previous timestep +/// * `temp_ext_air` - temperature of external air, in deg C +/// * `gains_internal` - total internal heat gains, in W +/// * `gains_solar` - directly transmitted solar gains, in W +/// * `gains_heat_cool` - gains from heating (positive) or cooling (negative), in W +/// * `f_hc_c` - convective fraction for heating/cooling +/// * `vent_extra_h_ve` - additional ventilation heat transfer coeff in response +/// to high internal temperature +/// * `throughput_factor` - proportional increase in ventilation rate due to +/// over-ventilation requirement +/// * `print_heat_balance` - flag to record whether to return the heat balance outputs +/// +/// Temperatures are calculated by solving (for X) a matrix equation A.X = B, where: +/// A is a matrix of known coefficients +/// X is a vector of unknown temperatures +/// B is a vector of known quantities +/// +/// Each row in vector X is a temperature variable - one for each node in each +/// building element plus the internal air temperature in the zone. +/// +/// Each row of matrix A contains the coefficients from the heat balance equations +/// for each of the nodes in each building element, plus one row for the heat +/// balance equation of the zone. +/// +/// Each column of matrix A contains the coefficients for a particular temperature +/// variable (in same order that they appear in vector X). Where the particular +/// temperature does not appear in the equation this coefficient will be zero. +/// +/// Note that for this implementation, the columns and rows will be in corresponding +/// order, so the heat balance equation for node i will be in row i and the +/// coefficients in each row for the temperature at node i will be in column i. +/// +/// Each row of vector B contains the other quantities (i.e. those that are not +/// coefficients of the temperature variables) from the heat balance equations +/// for each of the nodes in each building element, plus one row for the heat +/// balance equation of the zone, in the same order that the rows appear in matrix +/// A. +// fn calc_temperatures( +// delta_t: f64, +// temp_prev: Vec, +// temp_ext_air: f64, +// gains_internal: f64, +// gains_solar: f64, +// gains_heat_cool: f64, +// f_hc_c: f64, +// vent_extra_h_ve: Option, +// throughput_factor: Option, +// ) -> (Vec, HeatBalance) { +// } + +/// Optimised heat balance solver +// +/// ## Arguments +/// * `coeffs` - full matrix of coefficients for the heat balance eqns +/// * `rhs` - full vector of values that are not temperatures or coefficients +/// (i.e. terms on right hand side of heat balance eqns) +/// * `no_of_temps` - number of unknown temperatures (each node in each +/// building element + 1 for internal air) to be +/// solved for +/// * `building_elements` - the building elements of the zone in question +/// * `passed_zone_idx` - +/// +/// The heat balance equations from BS EN ISO 52016-1:2017 are expressed as a matrix equation and +/// solved simultaneously. While this provides a generic calculation procedure that works for an +/// arbitrary number of nodes (N), it also has a runtime proportional to N^3 which means that more +/// complex buildings can take a long time to simulate. However, many of the nodes are known not to +/// interact (e.g. the node in the middle of one wall has no heat transfer with the node in another +/// wall) and therefore we do not require the full flexibility of the matrix approach to solve for +/// every node temperature. The only part of the heat balance calculation where this flexibility is +/// needed is in the interaction between internal air and internal surfaces, so the calculation of the +/// other node temperatures can be removed from the matrix equation using algebraic substitution. +/// +/// Consider generic heat balance eqns for a 4-node element: +/// +/// A1a + B1b = Z1 # Heat balance at node a (external surface node) +/// A2a + B2b + C2c = Z2 # Heat balance at node b (inside node) +/// B3b + C3c + D3d = Z3 # Heat balance at node c (inside node) +/// C4c + D4d + J4j + K4k + Y4y = Z4 # Heat balance at node d (internal surface node) +/// where: +/// - a, b, c and d are the node temperatures in the building element to be solved for +/// - j and k are the node temperatures for the internal surfaces of other elements +/// - y is the internal air temperature +/// - A1 is the coefficient for temperature a in equation 1, A2 is the coefficient for temperature a in +/// equation 2, etc. +/// - Z1, Z2, etc. are the terms in equation 1, 2, etc. that are not the temperatures to be solved for +/// or their coefficients (i.e. Z1, Z2 etc. are the terms on the RHS of the equations) +/// +/// The heat balance equation for node a (external surface node) can be rearranged to solve for a: +/// +/// A1a + B1b = Z1 +/// A1a = Z1 - B1b +/// a = (Z1 - B1b) / A1 +/// +/// Using the rearranged heat balance equation for node a, we can substitute a in the heat balance +/// equation for node b (next inside node) to eliminate a as a variable: +/// +/// A2a + B2b + C2c = Z2 +/// A2 * (Z1 - B1b) / A1 + B2b + C2c = Z2 +/// +/// Rearranging to consolidate the occurances of b gives a new heat balance equation for b: +/// +/// A2 * Z1 / A1 + A2 * (- B1b) / A1 + B2b + C2c = Z2 +/// b * (B2 - A2 * B1 / A1) + A2 * Z1 / A1 + C2c = Z2 +/// (B2 - A2 * B1 / A1) * b + C2c = Z2 - A2 * Z1 / A1 +/// +/// This new heat balance equation can then be expressed in terms of modified versions of B2 and Z2: +/// +/// B2'b + C2c = Z2' +/// where: +/// B2' = B2 - B1 * A2 / A1 +/// Z2' = Z2 - Z1 * A2 / A1 +/// +/// The process can then be repeated, rearranging this new heat balance equation to solve for b and +/// then substituting into the heat balance equation for c: +/// +/// b = (Z2' - C2c) / B2' +/// +/// C3'c + D3d = Z3' +/// where: +/// C3' = C3 - C2 * B3 / B2' +/// Z3' = Z3 - Z2' * B3 / B2' +/// +/// And repeated again to generate a new equation for node d: +/// +/// c = (Z3' - D3d) / C3' +/// +/// D4'd + J4j + K4k + Y4y = Z4' +/// where: +/// D4' = D4 - D3 * C4 / C3' +/// Z4' = Z4 - Z3' * C4 / C3' +/// +/// At this point, we have reached the internal surface and we need the flexibility of the matrix +/// approach, but we have reduced the number of nodes to be solved for by 3. The process can be +/// repeated for each building element, and once the matrix solver has solved for the internal surface +/// temperatures, we can then go back through the other nodes, using the rearranged heat balance +/// equations that solve for (in this case) c, b and a. +/// +/// In order to deal with different building elements having different numbers of nodes, we can express +/// the above relationships generically: +/// +/// temperature[i] = (Z_adjusted[i] - coeff[i][i+1] * temperature[i+1]) / coeff_adjusted[i][i] +/// +/// coeff_adjusted[i][i] = coeff[i][i] - coeff[i-1][i] * coeff[i][i-1] / coeff_adjusted[i-1][i-1] +/// Z_adjusted[i] = Z[i] - Z_adjusted[i-1] * coeff[i][i-1] / coeff_adjusted[i-1][i-1] +/// +/// where i is the number of the node and its heat balance equation (counting from the external surface +/// to the internal surface), e.g. coeff[i-1][i] would be the coeffient for the temperature of the +/// current node in the heat balance equation for the previous node. +/// +/// The optimised calculation procedure is therefore: +/// - Loop over nodes, from external surface to internal surface, and calculate adjusted coeffs and RHS +/// for each heat balance eqn +/// - Construct matrix eqn for inside and air nodes only +/// - Solve heat balance eqns for inside and air nodes using normal matrix solver +/// - Loop over nodes, from internal inside node (i.e. inside node nearest to the internal surface) to +/// external surface, and calculate temperatures in sequence +fn fast_solver( + coeffs: DMatrix, + rhs: DMatrix, + no_of_temps: u32, + building_elements: &Vec, + element_positions: Vec<(usize, usize)>, + passed_zone_idx: usize, +) -> Vec { + // Init matrix with zeroes + // Number of rows in matrix = number of columns + // = total number of nodes + 1 for overall zone heat balance (and internal air temp) + let mut coeffs_adj: DMatrix = DMatrix::zeros(no_of_temps as usize, no_of_temps as usize); + + // Init vector_b with zeroes (length = number of nodes + 1 for overall zone heat balance) + let mut rhs_adj: DMatrix = DMatrix::zeros(1usize, no_of_temps as usize); + + // Init matrix with zeroes + // Number of rows in matrix = number of columns + // = total number of internal surface nodes + 1 for internal air node + let num_rows_cols_optimised = building_elements.len() + 1; + let zone_idx = num_rows_cols_optimised - 1; + let mut matrix_a: DMatrix = + DMatrix::zeros(num_rows_cols_optimised, num_rows_cols_optimised); + + // Init vector_b with zeroes (length = number of internal surfaces + 1 for air node) + let mut vector_b: DMatrix = DMatrix::zeros(1usize, num_rows_cols_optimised); + + for (el_idx, eli) in building_elements.iter().enumerate() { + let (idx_ext_surface, idx_int_surface) = element_positions[el_idx]; + + // No adjusted coeffs and RHS for external surface heat balance eqn + coeffs_adj[(idx_ext_surface, idx_int_surface)] = coeffs[(idx_ext_surface, idx_int_surface)]; + rhs_adj[idx_ext_surface] = rhs[idx_ext_surface]; + + // Loop over nodes, from inside node adjacent to external surface, to internal surface + for idx in (idx_ext_surface + 1)..(idx_int_surface + 1) { + // Calculate adjusted coeffs and RHS for each heat balance eqn + coeffs_adj[(idx, idx)] = coeffs[(idx, idx)] + - coeffs[(idx - 1, idx)] * coeffs[(idx, idx - 1)] / coeffs_adj[(idx - 1, idx - 1)]; + rhs_adj[idx] = rhs[idx] + - rhs_adj[idx - 1] * coeffs[(idx, idx - 1)] / coeffs_adj[(idx - 1, idx - 1)]; + } + + // Construct matrix eqn for internal surface nodes only (and air node, after this loop) + matrix_a[(el_idx, el_idx)] = coeffs_adj[(idx_int_surface, idx_int_surface)]; + vector_b[el_idx] = rhs_adj[idx_int_surface]; + + for (el_idx_other, _) in building_elements.iter().enumerate() { + if el_idx == el_idx_other { + continue; } - ThermalBridging::Bridges(bridges) + let idx_other_int_surface = element_positions[el_idx_other].1; + matrix_a[(el_idx, el_idx_other)] = coeffs[(idx_int_surface, idx_other_int_surface)]; + } + + // Add coeff for air temperature to this element's internal surface heat balance eqn + matrix_a[(el_idx, zone_idx)] = coeffs[(idx_int_surface, passed_zone_idx)]; + // Add coeff for this element's internal surface temp to the air node heat balance eqn + matrix_a[(zone_idx, el_idx)] = coeffs[(passed_zone_idx, idx_int_surface)]; + } + + // Add rest of air node heat balance eqn to matrix + // Coeffs for temperatures other than the air temp are added in the loop above + matrix_a[(zone_idx, zone_idx)] = coeffs[(passed_zone_idx, passed_zone_idx)]; + vector_b[zone_idx] = rhs[passed_zone_idx]; + + // Solve heat balance eqns for inside and air nodes using normal matrix solver + // Solve matrix eqn A.X = B to calculate vector_x (temperatures) + let vector_x = matrix_a.cholesky().unwrap().solve(&vector_b); // use Cholesky solver - may be able to speed this up with static matrices + + // Init temperature with zeroes (length = number of nodes + 1 for overall zone heat balance) + let mut temperatures = vec![0.0; no_of_temps as usize]; + temperatures[passed_zone_idx] = vector_x[zone_idx]; + + // Populate node temperature results for each building element + for (el_idx, _) in building_elements.iter().enumerate() { + let (idx_ext_surface, idx_int_surface) = element_positions[el_idx]; + + // Populate internal surface temperature result + temperatures[idx_int_surface] = vector_x[el_idx]; + + // Loop over nodes, from internal inside node (i.e. inside node nearest to the + // internal surface) to external surface, and calculate temperatures in sequence + for idx in (idx_ext_surface..idx_int_surface).rev() { + temperatures[idx] = (rhs_adj[idx] - coeffs[(idx, idx + 1)] * temperatures[idx + 1]) + / coeffs_adj[(idx, idx)]; } - Value::Number(n) => ThermalBridging::Number(n.as_f64().unwrap()), - _ => panic!("thermal bridging input was not in an expected format"), } + + temperatures +} + +#[derive(Clone)] +struct NamedBuildingElement { + name: String, + element: BuildingElement, +} + +struct HeatBalanceAirNode { + solar_gains: f64, + internal_gains: f64, + heating_or_cooling_system_gains: f64, + energy_to_change_internal_temperature: f64, + heat_loss_through_thermal_bridges: f64, + heat_loss_through_infiltration: f64, + heat_loss_through_ventilation: f64, + fabric_heat_loss: f64, +} + +struct HeatBalanceInternalBoundary { + fabric_int_air_convective: f64, + fabric_int_sol: f64, + fabric_int_int_gains: f64, + fabric_int_heat_cool: f64, +} + +struct HeatBalanceExternalBoundary { + solar_gains: f64, + internal_gains: f64, + heating_or_cooling_system_gains: f64, + thermal_bridges: f64, + ventilation: f64, + infiltration: f64, + fabric_ext_air_convective: f64, + fabric_ext_air_radiative: f64, + fabric_ext_sol: f64, + fabric_ext_sky: f64, + opaque_fabric_ext: f64, + transparent_fabric_ext: f64, + ground_fabric_ext: f64, + ztc_fabric_ext: f64, + ztu_fabric_ext: f64, +} + +struct HeatBalance { + air_node: HeatBalanceAirNode, + internal_boundary: HeatBalanceInternalBoundary, + external_boundary: HeatBalanceExternalBoundary, } #[cfg(test)] diff --git a/src/core/units.rs b/src/core/units.rs index 1479958..35bf9b7 100644 --- a/src/core/units.rs +++ b/src/core/units.rs @@ -2,3 +2,8 @@ pub const WATTS_PER_KILOWATT: u32 = 1000; pub const LITRES_PER_CUBIC_METRE: u32 = 1000; pub const SECONDS_PER_HOUR: u32 = 3600; pub const DAYS_IN_MONTH: [u32; 12] = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]; + +pub fn kelvin_to_celsius(temp_k: f64) -> f64 { + assert!(temp_k >= 0.0); + temp_k - 273.15 +} diff --git a/src/external_conditions.rs b/src/external_conditions.rs index 2a5a061..1d0fc5f 100644 --- a/src/external_conditions.rs +++ b/src/external_conditions.rs @@ -31,7 +31,7 @@ pub struct ShadingObject { distance: f64, } -#[derive(Clone, Debug, Deserialize)] +#[derive(Copy, Clone, Debug, Deserialize)] pub struct WindowShadingObject { #[serde(rename(deserialize = "type"))] object_type: WindowShadingObjectType, diff --git a/src/input.rs b/src/input.rs index 6d7c0d9..f8c8305 100644 --- a/src/input.rs +++ b/src/input.rs @@ -684,7 +684,7 @@ pub struct ZoneLighting { efficacy: f64, } -#[derive(Debug, Deserialize)] +#[derive(Clone, Debug, Deserialize)] #[serde(tag = "type", deny_unknown_fields)] pub enum BuildingElement { #[serde(rename(deserialize = "BuildingElementOpaque"))] @@ -753,7 +753,7 @@ pub enum BuildingElement { }, } -#[derive(Debug, Deserialize)] +#[derive(Copy, Clone, Debug, Deserialize)] pub enum MassDistributionClass { D, E,