Development of a Space Heating Model Suitable for the Automated Model Generation of Existing Multifamily Buildings—A Case Study in Nordic Climate

Abstract

Building energy performance modeling is essential for energy planning, management, and efficiency. This paper presents a space heating model suitable for auto-generating baseline models of existing multifamily buildings. Required data and parameter input are kept within such a level of detail that baseline models can be auto-generated from, and calibrated by, publicly accessible data sources. The proposed modeling framework consists of a thermal network, a typical hydronic radiator heating system, a simulation procedure, and data handling procedures. The thermal network is a lumped and simplified version of the ISO 52016-1:2017 standard. The data handling consists of procedures to acquire and make use of satellite-based solar radiation data, meteorological reanalysis data (air temperature, ground temperature, wind, albedo, and thermal radiation), and pre-processing procedures of boundary conditions to account for impact from shading objects, window blinds, wind- and stack-driven air leakage, and variable exterior surface heat transfer coefficients. The proposed model was compared with simulations conducted with the detailed building energy simulation software IDA ICE. The results show that the proposed model is able to accurately reproduce hourly energy use for space heating, indoor temperature, and operative temperature patterns obtained from the IDA ICE simulations. Thus, the proposed model can be expected to be able to model space heating, provided by hydronic heating systems, of existing buildings to a similar degree of confidence as established simulation software. Compared to IDA ICE, the developed model required one-thousandth of computation time for a full-year simulation of building model consisting of a single thermal zone. The fast computation time enables the use of the developed model for computation time sensitive applications, such as Monte-Carlo-based calibration methods.

1. Introduction

Building energy performance modeling is essential for energy planning, management, and efficiency. To improve energy efficiency of building stocks, knowledge on how each building is currently performing is needed, often referred to as baselines. These baselines can for example be used to compare and rank building performance, estimate energy and cost saving potential of energy efficiency measures, and detect faults when systems are malfunctioning. Energy use of a building depends on aspects such as physical and thermal properties, climate, control and operational schemes, occupants’ presence, and behavior. Information regarding these aspects is often uncertain and scattered, sometimes missing altogether, and establishing the energy performance of buildings can therefore be a tedious task [1,2].

Building energy modeling is often categorized into two major approaches: law-driven (forward) and data-driven (inverse) [1,3]. Law-driven models apply a given set of physical laws that govern the system. Conversely, data-driven models use system behavior as a predictor for system performance. An advantage of data-driven models is that they generally require less input, while law-driven models can offer more flexibility regarding the fine adjustment of a building’s subsystems and components [1]. Examples of building energy simulation software that utilize a law-driven approach include EnergyPlus [4] and IDA ICE (https://www.equa.se/en/ida-ice). These types of software calculate building energy consumption based on detailed building and environmental information such as climate, operation schedules, building geometrics, and the shading/sheltering impact of the surroundings [1,2]. Examples of using law-driven simulation software for parameter identification can be found in Kang and Krarti [5], Tian et al. [6]. Data-driven modeling approaches have gained much research attention in recent years. The review article by Amasyali and El-Gohary [7] identified more than 50 articles using data-driven approaches. Weaknesses with data-driven approaches are that these may not perform well outside their training range and that statistical models yield limited knowledge of under-laying aspects that govern energy use. The main weaknesses of law-driven approaches are that they require detailed input and tend to model how the building is designed to operate and not how it actually operates.

Gray-box modeling combines law-driven and data-driven approaches, thereby leveraging the advantages and minimizing the disadvantages of both approaches, see for example Amasyali and El-Gohary [7], Bacher and Madsen [8], Boodi et al. [9]. In gray-box models, internal parameter and equations are physically interpretable and thus provide knowledge about the modeled building. Gray-box models are calibrated to better match a building’s actual operation performance, but can also provide physically reasonable results outside the training data range. The authors in Boodi et al. [9] recognized a need for more research efforts regarding model structures for gray box models for residential buildings.

The ISO 52016-1:2017 [10] standard, which replaces the older ISO 13790:2008 [11] standard, presents a set of calculation methods for a building energy needs and internal air temperatures. The hourly method described in the standard proposes a system of linear equations that model heat transfer through opaque and transparent components of the envelope and air exchange between the internal and external environments. The calculation result is hourly internal air and component temperatures and heating and cooling loads. Each construction component (e.g., roof, windows, and walls) is modeled as serially connected RC (resistance and capacitance) thermal networks. Compared to the old ISO 13790:2008 [11] standard, the new hourly method is a much closer representation of algorithms used in whole building simulation software, such as IDA ICE or EnergyPlus. The simplistic hourly method of the now deprecated ISO 13790:2008 [11] has been employed and explored in several scientific works, for example [12,13,14,15]. To the best of our knowledge, there has not yet been any scientific works employing the newer ISO 52016-1:2017 [10] standard.

Law-driven approaches require detailed input regarding the physical properties of the building. Most existing buildings were built before the era of modern BIM (building information modeling) and thus lack such information in digitized form. However, 3D shape representation of existing building can be achieved through the use of aerial images or LiDAR (light detection and ranging) or a combination of both [16]. Such 3D shape representations allow for the determination of buildings’ footprint, height, orientation, and surface areas. The Swedish Land Survey (Lantmäteriet) has collected LiDAR data for almost all of Sweden [17]. Parameters describing buildings’ thermal properties and technical system are, however, usually not readily accessible in a digitized format. Gray-box modeling approaches have proven successful for estimating such uncertain or unknown parameters [7,8,9].

Energy meter readings of high resolution are required for gray-box modeling approaches. In Sweden, since 2015, energy bills have been required to be based on actual energy consumption [18]. Therefore, most Swedish district heating system operators have automatic meter reading systems installed that gather hourly or sub-hourly readings to centralized databases. Most Swedish cities have extensive district heating networks, and district heating operators have an approximately 90% share of the heating market for multifamily buildings [19]. Approximately 94% of Swedish multifamily buildings have hydronic heating systems, and approximately 94% of these have thermostatic control valves installed [20].

From a Swedish multifamily building portfolio owner’s perspective (as energy strategist at Eskilstuna Kommunfastighet AB and corresponding author responsible for the energy performance of a 7000 apartment portfolio), it would be beneficial to have an energy performance model of each building—a model revealing both current and potential energy performance of the building. Such a model needs to be law-driven and detailed enough so that its parameters reveal meaningful and actionable information. However, the model also needs to be simple enough, so that basic models of existing buildings can be auto-generated from readily available data sources. It also needs to have a structure that allows calibration with actual meter readings, to ensure it represents the actual performance of the building and not only the intended/designed performance. A practical model fulfilling the above-mentioned requirements is currently missing: detailed simulation software (e.g., IDA ICE or EnergyPlus) requires too detailed input to be readily auto-generated and are too slow and complex to be readily calibrated, while the data-driven and gray-box models proposed in the literature are not detailed enough and too far from physics to reveal meaningful interpretable information. This paper presents a space heating model aiming at filling the identified gap.

Methodology and Outline

The ISO 52016-1:2017 [10] standard employs an RC network to calculate hourly internal air temperatures, hourly heating and cooling loads, and space heating/cooling needs by summing the loads over certain time periods (e.g., monthly). This article presents the development of a thermal RC network that is based on the ISO 52016-1:2017 standard. However, our implementation is further simplified by lumping building elements by type and decreasing the number of temperature nodes in opaque building elements. Since the ISO 52016-1:2017 omits the conditions of use, the standard calculates the minimum energy need when assuming perfect temperature control and no system losses. To enable the calculation of actual energy use for space heating, comparable to metered values, a typical Swedish hydronic radiator heating system model [21,22] was implemented. Moreover, the following boundary condition procedures are proposed: a procedure to acquire satellite-based solar radiation data and conversion to a suitable format; a procedure to obtain and utilize meteorological reanalysis data (air temperature, ground temperature, wind, and thermal radiation and surroundings); procedures to account for impact from shading objects, window blinds, and wind- and stack-driven air infiltration, and variable exterior surface heat transfer coefficients. The developed model (including the RC network, heating system, and pre-processing steps) is from here on called ISO14N, where “ISO” denotes that it is based on the ISO 52016-1:2017 [10] standard and “14N” denotes the 14 temperature nodes (and 14 system equations). Figure 1 shows an overview of the proposed ISO14N modeling framework: from data input to pre-processing the boundary conditions data and constructing the RC network and to the simulation procedure.

As a validation, the simulation results were compared with results obtained with the well established and validated (https://www.equa.se/en/ida-ice/validation-certifications) simulation software IDA ICE. The main motivation behind using IDA ICE instead of the verification case provided by the ISO 52016-1:2017 [10] was that it also enabled testing the procedures that are not derived from the standard. Although the proposed RC network is generic, the use of a continuously operating hydronic radiator heating system and a constant air flow ventilation system limits the application of the proposed model to mainly multifamily buildings in Nordic climates. Therefore, the work is delimited to energy use for space heating for continuously heated multifamily buildings. Domestic hot water use, internal heat gains, and variable air flow ventilation are left for further work. Parallel work developing the calibration procedure, beyond the scope of this article, has been conducted and has influenced the design choices. The construction of the linear equation system of the RC network, the heating system, and the simulation procedure was implemented using the probabilistic programming language Stan [23], which parses to C++ code before compiling. The pre-processing procedures and the data input handling were implemented in the R statistical language based package dplyr [24]. Code can be requested from the corresponding author. The ISO 52010:2017 [25] standard (solar irradiance at an arbitrarily tilted and oriented surface) calculation procedure was implemented in C++ and published as an open source R-package (https://github.com/lukas-rokka/solarCalcISO52010).

The main features of the proposed ISO14N model are introduced in Section 2, while the pre-processing of the boundary conditions data is presented in Section 3. Section 4 describes the required data input and used data sources. Section 5.1, Section 5.2, Section 5.3, Section 5.4, Section 5.5, Section 5.6, Section 5.8, and Section 5.9 compare simulation results of the developed ISO14N model with results obtained with the IDA ICE. Section 5.7 show results from the proposed shading reduction. Details on constructing the linear equation system of the RC network are given in Appendix A.

2. The ISO14N Space Heating Model

In this paper, a one-zone model representation for the whole building is illustrated. Figure 2 displays the proposed one-zone RC network. All external walls $ew$ are modeled as a lumped 3-node element, all window glazing $gl$ are modeled as a lumped 2-node element, and the external roof $rf$ and the ground floor $gf$ are modeled with 3-node elements. The remaining internal mass $im$ (internal walls, intermediate floors, and adiabatic external walls) are represented with a 2-node element. Thus, the thermal network consists of 14 unknown temperature nodes (including the internal air node ${\theta }_{int}$). The thermal model can be represented as a system of linear equations:

$$\mathbf{A}\times \mathbf{X}=\mathbf{b}$$

(1)

where $\mathbf{A}$ is a $14\times 14$ square matrix holding system coefficients, $\mathbf{X}$ is the state vector holding the 14 node temperatures, and $\mathbf{b}$ is a vector holding determined terms (node temperatures from the previous time step and known boundary condition values). Appendix A shows the details of the model matrices. For every time step $\mathrm{t}$, the system of equations needs to be solved to compute the new states:

$${\mathbf{X}}_{\mathrm{t}}=\mathbf{A}/{\mathbf{b}}_{\mathrm{t}}={\mathbf{A}}^{-1}\times {\mathbf{b}}_{\mathrm{t}}.$$

(2)

Matrix inversion is computationally intensive. However, if matrix $\mathbf{A}$ only consists of constants, this computation step is only required once before the iterative simulation procedure. Vector $\mathbf{b}$ holds terms that can be treated as pre-processed data input (e.g., solar irradiance) and terms that need to be calculated at each simulation step due to dependence on the state of the node temperatures (e.g., thermal storage in nodes). Figure 1c shows the simulation procedure in pseudo code.

2.1. Geometrical and Thermal Properties

In this paper, the calculation is conducted normalized to per square meter floor (a deviation from the standard). One benefit is that geometrical data can be given as ratios (r) between the total interior surface area of each building element type and the total floor area (unit ${\mathrm{m}}_{\mathrm{s}}/{\mathrm{m}}_{\mathrm{fl}}$). These ratios can be seen as weighting factors describing how large an impact one type of building element has on the average thermal balance of the whole building. It is straightforward to estimate these ratios from readily available building data:

$$\begin{array}{cc}\hfill r_{rf}& =r_{gf}=1/N_{fl}\hfill \\ \hfill r_{ew}& =P·H_{fl}·N_{fl}/A_{fl}-r_{gl}\hfill \\ \hfill r_{im}& =(2-2/N_{fl})+1.5\hfill \end{array}$$

(3)

where $N_{fl}$ is the number of floors of the building, $A_{fl}$ is the total floor area of the building, P is the building perimeter, and $H_{fl}$ is the average internal floor height. The term $(2-2/N_{fl})$ in $r_{im}$ calculates the surface area of internal floors and ceilings per total floor area and will result in a value between 0 and 2, depending on the number of floors. The additional 1.5 constant represents the internal walls, and it is based on the old ISO 13790:2008 [11] standard (where total interior surface area per total floor area is given as 4.5). The glazing-to-floor area ratio $r_{gl}$ is a property that is often regulated in standards and building codes, typically ranging between 0.1 and 0.2 for Swedish multi-family buildings.

Weighted surface heat transfer coefficients for each building element type are as follows: $el=[rf,ew,gl,im,gf]$, which are abbreviations of [roof, external walls, glazing, internal mass, ground floor]:

$$\begin{array}{ccccc}{H}_{se;rf;\mathrm{t}}=r_{rf}·(h_{ce;rf;\mathrm{t}}+(1-F_{sky;hor})·4.14),\hfill & \hfill & H_{ci;rf}=r_{rf}·0.7,\hfill & \hfill & H_{ri;rf}=r_{rf}·5.13\hfill \\ H_{se;ew;\mathrm{t}}=r_{ew}(·h_{se;ew;\mathrm{t}}+(1-F_{sky;ver})·4.14),\hfill & \hfill & H_{ci;ew}=r_{ew}·2.5,\hfill & \hfill & H_{ri;ew}=r_{ew}·5.13\hfill \\ H_{se;gl;\mathrm{t}}=r_{gl}·(h_{se;gl;\mathrm{t}}+(1-F_{sky;ver})·4.14),\hfill & \hfill & H_{ci;gl}=r_{gl}·2.5,\hfill & \hfill & H_{ri;gl}=r_{gl}·5.13\hfill \\ & \hfill & H_{ci;im}=r_{im}·2.5,\hfill & \hfill & H_{ri;im}=r_{im}·5.13\hfill \\ H_{se;gf}=r_{gf}·h_{gr;vi},\hfill & \hfill & H_{ci;gf}=r_{gf}·5.0,\hfill & \hfill & H_{ri;gf}=r_{gf}·5.13\hfill \end{array}$$

(4)

where $H_{se;el;\mathrm{t}}$ [$\mathrm{W}/(\mathrm{K}·{\mathrm{m}}_{\mathrm{fl}}^2)$] represents normalized variable heat transfer coefficients (convective + radiative) between the exterior surface of the elements and external air at time step t, $h_{ce;el;\mathrm{t}}$ [$\mathrm{W}/(\mathrm{K}·{\mathrm{m}}_{\mathrm{s}}^2)$] represents exterior surface convective heat transfer coefficients at time step t (see Section 3.1), $F_{sky;hor}$ and $F_{sky;ver}$ are the sky view factors for horizontal and vertical elements respectively, $h_{re}$ is the surface exterior radiative heat transfer coefficient, $h_{gr;vi}$ is the heat transfer coefficient towards a virtual ground layer (see Equation (8)), $H_{ci;el}$ represents weighted convective heat transfer coefficients between the interior surface side of the elements and the internal air node, and $H_{ri;el}$ represents weighted radiative heat transfer coefficients of interior surface nodes. The constants are taken from ISO 52016-1:2017 [10] (Table 25) and describe conventional convective and radiative surface heat transfer coefficients for interior or exterior surfaces oriented upwards, horizontally, or downward.

In ISO 52016-1:2017 [10], opaque building elements (walls, roof, etc.) are by default described as 5-node and window glazing as 2-node elements. The distribution of the thermal resistances are slightly weighted towards the center of the building elements, while heat capacity nodes are weighted depending on chosen class (mass concentrated externally, internally, inside, or equally). In this paper, 3-node elements are used for opaque elements; therefore, there is less flexibility in distributing the thermal properties over the elements. The thermal resistance is equally distributed:

$$\begin{array}{ccc}{H}_{1;rf}=H_{2;rf}=2/R_{rf},\hfill & \hfill & H_{1;gl}=1/R_{gl}\hfill \\ H_{1;ew}=H_{2;ew}=2/R_{ew},\hfill & \hfill & H_{1;im}=1/R_{im}\hfill \end{array}$$

(5)

where the floor area normalized thermal resistance R [${\mathrm{m}}_{\mathrm{fl}}^2·{}^{\circ }\mathrm{C}/\mathrm{W}$] is calculated as

$$R_{el}=1/(r_{el}·U_{el})-1/(H_{ci;el}+H_{ri;el})-1/H_{se;el}$$

(6)

where $U_{el}$ [$\mathrm{W}/{(}^{\circ }\mathrm{C}·{\mathrm{m}}_{\mathrm{s}}^2)$] is the thermal transmittance of building element $el$. The heat capacity of 3-node building elements roof $rf$ and external walls $ew$ is distributed according to

Table 1. Node heat capacity distribution for the the 3-node building elements roof $rf$ and external walls $ew$. Class I: concentrated on interior side; Class E: concentrated on exterior side; Class IE: divided on interior and exterior side; Class D: equally distributed; Class M: inside/centered.

Node	Class I	Class E	Class IE	Class D	Class M
1	0.10	0.50	0.40	0.33	0.10
2	0.40	0.40	0.20	0.33	0.80
3	0.50	0.10	0.40	0.33	0.10

. The table is based on the ISO 52016-1:2017 [10] but altered to suit 3-node elements. Table B.14 of ISO 52016-1:2017 [10] specifies ${\kappa }_m$ values according to five classes (from very light to very heavy, 14–70 $\mathrm{Wh}/{(}^{\circ }\mathrm{C}·{\mathrm{m}}_{\mathrm{s}}^2)$). The simulation procedure is conducted with normalized values, which are acquired by multiplying with surface ratios: $C_{el}=r_{el}·{\kappa }_{m;el}$.

The internal mass $im$ building element is configured for 24 h thermal admittance cycles:

$$\begin{array}{ccccc}{H}_{1;im}=1/1.0,\hfill & \hfill & {\kappa }_{1;im}=0.85·{\kappa }_{m;im},\hfill & \hfill & {\kappa }_{2;im}=0.15·{\kappa }_{m;im}.\hfill \end{array}$$

(7)

The ground floor element needs to account for effects of the ground, and its thermal resistance and heat capacity are distributed as follows:

$$\begin{array}{ccccc}{H}_{2;gf}=2/R_{gf},\hfill & \hfill & H_{1;gf}=1/(R_{gf}/2+R_{gr}),\hfill & \hfill & H_{se;gf}=r_{gr}·h_{gr;vi}\hfill \\ {\kappa }_{3;gf}={\kappa }_{m;gf}/2,\hfill & \hfill & {\kappa }_{2;gf}={\kappa }_{m;gf}/2,\hfill & \hfill & {\kappa }_{1;gf}={\kappa }_{gr}\hfill \end{array}$$

(8)

where $R_{gr}$ and ${\kappa }_{gr}$ represent the thermal resistance and heat capacity of a $0.5$ m thick soil layer (0.25 ${\mathrm{m}}_{\mathrm{s}}^2·\mathrm{K}/\mathrm{W}$ and 280 $\mathrm{Wh}/(\mathrm{K}·{\mathrm{m}}_{\mathrm{s}}^2)$ are used as default values), $R_{gf}$ is the thermal resistance of the ground floor including the effects of the ground, and $h_{gr;vi}$ is the thermal transmittance of a virtual ground layer. $R_{gf}$ and $h_{gr;vi}$ are calculated according to ISO 13370:2017 [26].

2.2. Hydronic Heating System

Heat dissipation from, in Sweden typically, a hydronic heating system with radiator panels, at time step t, can be approximated with the following non-linear equation [21,22]:

$${\varphi }_{hyd;\mathrm{t}}=H_{hyd}·{({\theta }_{hyd;lmtd;\mathrm{t}})}^{n_{hyd}}·u_{trv;\mathrm{t}}$$

(9)

where $H_{hyd}$ is the per floor square meter normalized radiator constant, $n_{hyd}$ is the radiator exponent (this depends on the design and type, but a value of 1.3 is commonly used for typical Swedish radiator panels [21,22]), $u_{trv;\mathrm{t}}$ is the control signal from the local thermostatic radiator valve(s) (TRV), and ${\theta }_{hyd;lmtd;\mathrm{t}}$ is the logarithmic mean temperature difference between the radiator and the internal air temperature, which is calculated as

$${\theta }_{hyd;lmtd;\mathrm{t}}=\frac{{\theta }_{hyd;sup;\mathrm{t}}-{\theta }_{hyd;ret;\mathrm{t}}}{ln\left({\theta }_{hyd;sup;\mathrm{t}}-{\theta }_{int;\mathrm{t}-1}\right)-ln\left({\theta }_{hyd;ret;\mathrm{t}}-{\theta }_{int;\mathrm{t}-1}\right)}$$

(10)

where ${\theta }_{hyd;sup;\mathrm{t}}$ and ${\theta }_{hyd;ret;\mathrm{t}}$ are the supply/return temperatures to/from the radiators, ${\theta }_{int;\mathrm{t}-1}$ is the internal air temperature from the previous time step (internal air temperature is unknown at time step t until the linear equations system has been solved). The supply temperature ${\theta }_{hyd;sup}$ is usually centrally controlled. For Swedish multifamily buildings, ${\theta }_{hyd;sup}$ is usually controlled by linear interpolation from a look-up table. Such a look-up table consists of value pairs specifying supply temperature set-points ${\theta }_{hyd;sup;set}$ at certain external air temperature ${\theta }_e$. an example.

For a hydronic heating system without local TRVs, the following empirical equation was derived to estimate the return temperature:

$$\begin{array}{cc}\hfill {\theta }_{hyd;ret;\mathrm{t}}& ={\theta }_{hyd;sup;\mathrm{t}}-b·{({\theta }_{hyd;sup;\mathrm{t}}-{\theta }_{int;\mathrm{t}-1})}^a\hfill \\ \hfill a& =n_{hyd}-\Delta {\theta }_{hyd;d}/200,b=\frac{\Delta {\theta }_{hyd;d}}{{({\theta }_{hyd;sup;d}-{\theta }_{int;set})}^a}\hfill \end{array}$$

(11)

where $\Delta {\theta }_{hyd;d}$ is the temperature drop between supply and return temperatures at design power output, ${\theta }_{hyd;sup;d}$ is the supply temperature at design power output, and ${\theta }_{int;set}$ is the internal air set-point. Equation (11) was derived empirically and chosen to suit common Swedish radiator configurations. Validation is provided as Supplementary Materials. The TRV(s) are modeled as proportional controllers (P-controllers):

$$u_{trv;\mathrm{t}}=0.5·u_{trv;\mathrm{t}-1}+0.5·max(0,min(1,({\theta }_{int;set}-{\theta }_{int;\mathrm{t}-1})/{\theta }_{trv;pb}))$$

(12)

where ${\theta }_{int;set}$ is the set-point for desired indoor temperature, ${\theta }_{trv;pb}$ is the proportional band of the TRV (typically ranging between 0.5–2.0 ${}^{\circ }\mathrm{C}$). Taking the average of current step and previous time step dampens oscillations in the control signal.

2.3. Ventilation

Heat transfer due to active ventilation is modeled as

$${\varphi }_{ve;\mathrm{t}}=\kappa {\rho }_a·Q_{ve;\mathrm{t}}·(1-{\eta }_{ve})·({\theta }_{int;\mathrm{t}-1}-{\theta }_{e;\mathrm{t}})$$

(13)

where $\kappa {\rho }_a$ is the heat capacity of air per volume 1.21 $\mathrm{W}·\mathrm{s}/(\mathrm{l}·\mathrm{K})$, $Q_{ve;\mathrm{t}}$ is the specific air flow rate [$\mathrm{l}/(\mathrm{s}·{\mathrm{m}}_{\mathrm{fl}}^2$)], ${\eta }_{ve}$ is the temperature transfer efficiency of the heat recovery unit, ${\theta }_{int;\mathrm{t}-1}$ is the internal air temperature from the previous time step t, and ${\theta }_{e;\mathrm{t}}$ is the external air temperature at current time step t.

2.4. Infiltration

Infiltration is driven by pressure differences across the building envelope caused by unbalanced mechanical ventilation, wind, and air temperature difference between internal and external air (known as the stack or buoyancy effect). One of the most established single zone models for estimating air infiltration rates is the Alberta air infiltration model (AIM-2), originally presented by Walker and Wilson [27] and referred to as the “enhanced model” in the ASHRAE Handbook—Fundamentals [3]. The AIM-2 model was developed for dwellings and has proven to be fairly accurate at estimating air infiltration rates for dwellings [28]. Hayati et al. [29] successfully used and evaluated the AIM-2 model for Swedish churches, modeled as large single zone models. The proposed thermal model of this paper is implemented as a single zone model, but this single zone is a representation of the average of a whole building (that could be a height of several floors). Thus, the AIM-2 model is implemented with modifications as to make the model more suitable for multi-floor buildings. Compared to earlier representations [3,27,28,29], equations are presented such that parameters are separated from parts that can be treated as pre-calculable data input. Thermal power losses due to infiltration are calculated as

$${\varphi }_{inf;\mathrm{t}}=C_{inf}·Q_{inf;\mathrm{t}}^{*}·\kappa {\rho }_a·({\theta }_{int;\mathrm{t}-1}-{\theta }_{e;\mathrm{t}})$$

(14)

where $C_{inf}$ [$\mathrm{l}/(\mathrm{s}·{\mathrm{Pa}}^{\mathrm{n}}·{\mathrm{m}}_{\mathrm{fl}}^2)$] is the infiltration coefficient (often referred to as the flow coefficient), n is the flow exponent, $\kappa {\rho }_a$ is the heat capacity of air per volume [1.21 $\mathrm{Ws}/(\mathrm{l}·{}^{\circ }\mathrm{C})$], and $Q_{inf}^{*}$ is the potential specific air flow rate [${\mathrm{Pa}}^{\mathrm{n}}$] due to air infiltration. $Q_{inf}^{*}$ is pre-calculated (see Section 3.2), the flow exponent n is treated as a constant ($n=0.73$ [30]), while $C_{inf}$ needs to be provided as an input parameter.

2.5. Radiant Heat Gains Exterior Surfaces

Net radiant heat gains on exterior surfaces due to thermal radiation exchange with the sky dome and absorption of solar irradiance are expressed as

$$\begin{array}{cc}\hfill {\varphi }_{re;rf;\mathrm{t}}& =r_{rf}·({\alpha }_{sol}·I_{tot;hor;sh;\mathrm{t}}-{\varphi }_{sky;rf;\mathrm{t}})\hfill \\ \hfill {\varphi }_{re;ew;\mathrm{t}}& =r_{ew}·({\alpha }_{sol}·I_{tot;ver;sh;\mathrm{t}}-{\varphi }_{sky;ew;\mathrm{t}})\hfill \\ \hfill {\varphi }_{re;gl;\mathrm{t}}& =r_{gl}·(-{\varphi }_{sky;gl;\mathrm{t}})\hfill \end{array}$$

(15)

where r is the weighting ratio of the element, ${\alpha }_{sol}$ is the solar absorption coefficient of the exterior surface (0.5 is used as default), $I_{tot}$ is the calculated total hemispherical solar irradiance (see Section 3.3), ${\varphi }_{sky}$ is the extra thermal radiation to the sky, subscript $ver$ denotes vertical surfaces, subscript $hor$ denotes horizontal surfaces, and subscript $sh$ denotes shading. Thermal radiation exchange with the surrounding environment is accounted for in Section 3.1. The extra thermal radiation to the sky variable follows the definition used in [10] (positive values denote heat loss; negative values denote heat gains). However, an alternative calculation procedure is proposed so that it can be derived from climate model data sources:

$${\varphi }_{sky;el;\mathrm{t}}=F_{sky;el}·{ϵ}_{el}·(\sigma ·{({\theta }_{1;el;\mathrm{t}-1}+273.15)}^4-{\varphi }_{strd})$$

(16)

where $F_{sky;el}$ is the sky view factor for element $el$ (0.5 used for vertical elements external walls and glazing and 1.0 used for the horizontal roof element), $ϵ$ is the emissivity of the surfaces (0.9 used as default), ${\theta }_{1;el}$ is the exterior surface temperature of the element $el$, $\sigma$ is the Stefan–Boltzmann constant (5.67 × 10${}^{-8}$ $\mathrm{W}/({\mathrm{m}}^2·{\mathrm{K}}^4)$), 273.15 is the Kelvin to Celsius conversion constant, and ${\varphi }_{strd}$ is the “surface thermal radiation downwards” variable retrieved from the climate reanalysis dataset ERA5 (see Section 4.2).

2.6. Radiant Heat Gains Interior Surfaces

Interior surfaces are subjected to solar and thermal radiation, expressed as

$$\begin{array}{cc}\hfill {\varphi }_{ri;el;\mathrm{t}}& =f_{el}·(f_{r;sol}·{\varphi }_{sol;gl;\mathrm{t}}+f_{r;int}·{\varphi }_{int;\mathrm{t}}+f_{r;hyd}·{\varphi }_{hyd;\mathrm{t}})\hfill \end{array}$$

(17)

where $f_{el}$ is the surface area fraction of the element $el$, and $f_r$ is the fraction of radiative heat ($f_r=1-f_c$). Default convective fractions from the ISO 52016-1:2017 [10] standard: $f_{c;sol}=0.1$, $f_{c;int}=0.4$, and $f_{c;hyd}=0.4$. The surface area fractions $f_{el}$ are simply calculated by dividing the ratio (r) of each building element with the sum of all ratios of the five building element types.

2.7. Operative Temperature

The ISO 52016-1:2017 [10] calculates operative temperature as follows:

$$\begin{array}{cc}\hfill {\theta }_{op;\mathrm{t}}& =0.5·({\theta }_{int;\mathrm{t}}+{\theta }_{int;r;mn;\mathrm{t}})\hfill \end{array}$$

(18)

where ${\theta }_{int;r;mn}$ is the mean radiant temperature, calculated as the area weighted mean of interior surface temperatures of all building elements. However, Equation (18) omits the impact of radiative heating from the hydronic heating system. With radiator panels of the known surface area and surface temperatures, their contribution to the operative temperature can be area-weighted in the same manner as interior walls. An equation that does not require such detailed knowledge of the heating system is proposed:

$$\begin{array}{cc}\hfill {\theta }_{op;\mathrm{t}}& =0.5·({\theta }_{int;\mathrm{t}}+{\theta }_{int;r;mn;\mathrm{t}}+C_{hyd;r;mn}·f_{r;hyd}·{\varphi }_{hyd;\mathrm{t}})\hfill \end{array}$$

(19)

where $C_{hyd;r;mn}$ is a weighting and conversion coefficient expressing the impact of the radiative heating part of the hydronic heating system on the operative temperature. The unit of the $C_{hyd;r;mn}$ parameter is the same as that of a floor area normalized thermal resistance, ${\mathrm{m}}_{\mathrm{fl}}^2·{}^{\circ }\mathrm{C}/\mathrm{W}$.

3. Pre-Processing of Boundary Conditions Data

3.1. Convective Heat Transfer Coefficients for Weather-Exposed Exterior Surfaces

The ISO 52016-1:2017 [10] standard provides constant conventional surface heat transfer coefficients for exterior surfaces ($h_{ce}$). However, $h_{ce}$ varies strongly depending on weather conditions [31,32]. Here, we implement the correlation equations by Liu et al. [32], derived by computational fluid dynamic calculations. Compared to earlier works [31], Liu et al. [32] accounts for wind sheltering effects, and separate correlations factors are provided for vertical walls and the roof. Additionally, our implementation accounts for surface roughness by adopting parts of the DOE-2 convection model [4] method.

$$\begin{array}{cc}\hfill h_{ce;rf;\mathrm{t}}& =h_{c;n}+R_{f;rf}·\left(\sqrt[]{h_{c;n}^2+h_{c;f;hor;\mathrm{t}}^2}-h_{c;n}\right)\hfill \\ \hfill h_{ce;ew;\mathrm{t}}& =h_{c;n}+R_{f;ew}·\left(\sqrt[]{h_{c;n}^2+h_{c;f;ver;\mathrm{t}}^2}-h_{c;n}\right)\hfill \\ \hfill h_{ce;gl;\mathrm{t}}& =\sqrt[]{h_{c;n}^2+h_{c;f;ver;\mathrm{t}}^2}\hfill \end{array}$$

(20)

where $h_{c;n}$ is the natural convective heat transfer coefficient, $h_{c;f}$ is the forced (wind effects) convective heat transfer coefficient, and $R_f$ is the surface roughness multiplier. The surface roughness multiplier is set as smooth $R_{f;rf}=1.11$ for the roof and as medium rough $R_{f;ew}=1.52$ for external walls [4]. The authors in Liu et al. [32] calculated the natural convective heat transfer coefficient as depending on temperature difference between external air and the exterior surface of the building element: $h_{c;n;\mathrm{t}}=1.52·\Delta {\theta }_{\mathrm{t}}^{0.36}$ $\mathrm{W}/(\mathrm{K}·{\mathrm{m}}_{\mathrm{s}}^2)$. Here we assume a constant temperature difference of $\Delta \theta =5{}^{\circ }\mathrm{C}$, which results in a constant $h_{c;n}=2.7$ $\mathrm{W}/(\mathrm{K}·{\mathrm{m}}_{\mathrm{s}}^2)$. Assuming $h_{c;n}$ as a constant allows treating Equation (20) as pre-calculable data input. Wind has the strongest impact on the exterior surface convective heat transfer, and the forced heat transfer coefficients are calculated as

$$\begin{array}{cc}\hfill h_{c;f;ver;\mathrm{t}}& =F_{ww}·(3.39-5.03·{\lambda }_p)·U_{loc;ver}^{0.94}+(1-F_{ww})·(1.15+0.82·{\lambda }_p)·U_{loc;ver}^{0.94}\hfill \\ \hfill h_{c;f;hor;\mathrm{t}}& =(3.57+1.72·{\lambda }_p)·U_{loc;hor}^{0.84}\hfill \end{array}$$

(21)

where $F_{ww}$ is the fraction of windward-oriented exterior surfaces, ${\lambda }_p$ is the plan area density, and $U_{loc}$ is the local wind speed. $F_{ww}$ is assumed as a constant of 0.5 but could also be calculated as a variable depending on wind orientation using an approach similar to that taken for solar irradiance in Equation (29). The plan area density ${\lambda }_p$ represents the projected built area viewed from above to the total area in consideration [32]. Conceptually it relates to the wind sheltering factor used for the wind infiltration calculation (see Equation (26)), and here we present them in the same table,

Table 2. Plan area densities ${\lambda }_p$ [32], wind sheltering factors ${\lambda }_w$ [3], and solar shielding factors ${\lambda }_{sol}$ for five shelter classes: ranging from (1) no shelter to (5) a highly dense urban area (neighboring buildings closer than one building height).

Shelter Class	1	2	3	4	5
${\lambda }_p$	0.00	0.04	0.11	0.25	0.44
${\lambda }_w$	1.00	0.90	0.70	0.50	0.30
${\lambda }_{sol}$	3.00	2.50	1.90	1.30	0.70

, and categorize the ${\lambda }_p$ values used in [32] into the sheltering classes defined in Table 5, Chapter 16 ASHRAE [3]. The local wind speeds are calculated according to Equation (22), where the height H is set as the building height for horizontal surfaces and half the building height for vertical surfaces.

$$U_{loc}=U_{10\mathrm{m}}·1.59·{\left(\frac{H_{loc}}{\delta }\right)}^{\alpha }$$

(22)

where $U_{10\mathrm{m}}$ is the meteorological wind speed at 10 m height, the constant 1.59 describes the terrain conditions of the meteorological wind speed measurement site, $H_{loc}$ is the height above ground for local wind calculation point, and $\delta$ and $\alpha$ are atmospheric boundary layer coefficients for different terrain categories, provided, e.g., in Table 1 of Chapter 24 of the ASHRAE Handbook—Fundamentals [3] ($\delta =370$ and $\alpha =0.22$ for urban and suburban terrain).

3.2. Infiltration Potential

Potential air infiltration from wind-driven and stack-effect-driven infiltration is not additive, which in the AIM-2 [3] is modeled with the following superposition formula:

$$\begin{array}{cc}\hfill Q_{inf;\mathrm{t}}^{*}& ={\left({(Q_{s;\mathrm{t}}^{*})}^{\frac{1}{n}}+{(Q_{w;\mathrm{t}}^{*})}^{\frac{1}{n}}-0.33·{(Q_{s;\mathrm{t}}^{*}·Q_{w;\mathrm{t}}^{*})}^{\frac{1}{2n}}\right)}^n\hfill \end{array}$$

(23)

where the flow exponent n is treated as a constant ($n=0.73$ [30]), and $Q_s^{*}$ and $Q_w^{*}$ are the contributions of the stack and the wind effects. Compared to the original formula [3,27], Equation (23) is modified such that $Q^{*}=Q/C_{inf}$ [${\mathrm{Pa}}^{\mathrm{n}}$]. The contribution due to stack effect is calculated as

$$Q_{s;\mathrm{t}}^{*}=C_s·{\left({\rho }_a·g_0·H^{*}·\frac{|{\theta }_{int;set}-{\theta }_{e;\mathrm{t}}|}{{\theta }_{int;set}+273.15}\right)}^n$$

(24)

where $C_s$ [${(\mathrm{Pa}/\mathrm{K})}^{\mathrm{n}}$] is the stack coefficient, ${\rho }_a$ [$\mathrm{kg}/{\mathrm{m}}^3$] is the density of air, $g_0$ [$\mathrm{m}/{\mathrm{s}}^2$] is the acceleration of gravity, $H^{*}$ [$\mathrm{m}$] is the modified building’s ceiling height, ${\theta }_{int;set}$ is the internal temperature set point, and ${\theta }_{e;\mathrm{t}}$ is the external air temperature at time step t. The set point temperature is used instead of the actual internal air temperature so that Equation (24) can be treated as a pre-calculable data variable. The AIM-2 model was originally designed as a single zone model for residential building with up to three floor levels, assuming free air contact between the floor levels. Using the model as such for high rise buildings would probably overestimate the stack effect. Thus, a modified building height parameter $H^{*}$ is introduced:

$$H^{*}=\left\{\begin{array}{cc}{N}_{fl}·0.5·H_{fl}\hfill & \mathrm{if}{N}_{fl}<3\hfill \\ N_{fl}·0.5·H_{fl}/3\hfill & \mathrm{if}{N}_{fl}\ge 3\hfill \end{array}\right.$$

(25)

where $N_{fl}$ is the number of floor levels, and $0.5·H_{fl}$ is half the average floor height. The rationale behind dividing by 3 for non-dwelling buildings with three or more floor levels is that the stack effect of such buildings will be a mix of the per zone/floor level height-induced stack effect and the whole building stack stack effect due to elevator shafts and stairways. The contribution due to wind effect is calculated as

$$Q_{w;\mathrm{t}}^{*}=C_w·{\left(\frac{{\rho }_a·{\left({\lambda }_w·U_{loc}\right)}^2}{2}\right)}^n$$

(26)

where the $C_w$ [${(\mathrm{Pa}·{\mathrm{s}}^2/{\mathrm{m}}^2)}^{\mathrm{n}}$] is the wind coefficient, ${\lambda }_w$ is the wind sheltering factor given in Table 2, and $U_{loc}$ is the local wind speed at the building’s uppermost eaves height calculated with Equation (22).

In the AIM-2 air infiltration model [27], the stack and wind coefficients ($C_s$ and $C_w$) can be calculated as functions of leakage distribution, occurrence of flue, and the type of foundation (crawlspace, basement, or slab-on-grade) [27]. Assuming evenly distributed leakage, no flue, and slab-on-grade foundation, $C_s=0.25\left[{(\mathrm{Pa}/\mathrm{K})}^{\mathrm{n}}\right]$ and $C_w=0.22\left[{(\mathrm{Pa}·{\mathrm{s}}^2/{\mathrm{m}}^2)}^{\mathrm{n}}\right]$. The infiltration coefficient $C_{inf}$ is given in liters instead of cubic meters and normalized to square meter floor. Otherwise, it is the same coefficient as described in [3,27]. Thus, it can be derived from building air leakage databases or from measurements (such as blower door pressurization tests).

3.3. Solar Heat Gains

Solar heat gains transmitted through windows glazing ($gl$) are calculated as

$$\begin{array}{cc}\hfill {\varphi }_{sol;gl;\mathrm{t}}& =r_{gl}·g_{gl}·g_{bl}·I_{tot;ver;sh;\mathrm{t}}\hfill \end{array}$$

(27)

where $I_{tot;ver;sh;\mathrm{t}}$ [$\mathrm{W}/{\mathrm{m}}_{\mathrm{s}}^2$] is the (by surface orientation) weighted hemispherical solar irradiance on vertical surfaces including shading effects, $g_{bl}$ is the total solar energy transmittance of the window blinds (see Section 3.5), $g_{gl}$ is the total solar energy transmittance of the glazing, and $r_{gl}$ is the glazing-to-floor ratio. $I_{tot;ver;sh;\mathrm{t}}$ is calculated as

$$I_{tot;ver;sh;\mathrm{t}}=I_{dif;ver;\mathrm{t}}+I_{dir;ver;\mathrm{t}}·F_{sh;ver;\mathrm{t}}$$

(28)

where $I_{dif;ver;\mathrm{t}}$ and $I_{dir;ver;\mathrm{t}}$ are the total diffuse (inclusive ground reflectance and exclusive circum-solar) and the total direct (inclusive circum-solar) irradiance received on one square meter of weighted vertical surface area at time step t, and $F_{sh;ver;\mathrm{t}}$ is a shading factor (see Section 3.4). Building elements are lumped according to their type, so solar irradiance variables also need to be allocated and weighted accordingly:

$$I_{dif;ver;\mathrm{t}}=\sum _{k=1}^{n_{\gamma }}{F}_{{\gamma }_k}·I_{dif;ver;{\gamma }_k;\mathrm{t}}{I}_{dir;ver;\mathrm{t}}=\sum _{k=1}^{n_{\gamma }}{F}_{{\gamma }_k}·I_{dir;ver;{\gamma }_k;\mathrm{t}}$$

(29)

where $I_{dir;ver;{\gamma }_k;\mathrm{t}}$ is the total direct solar irradiance at a vertical surface with azimuth $\gamma$ and time step t, $F_{{\gamma }_k}$ is the surface area fraction of external walls at azimuth ${\gamma }_k$, and $n_{\gamma }$ is the total number of surfaces of different azimuths to loop through. The same procedure is used for $I_{dif;ver;\mathrm{t}}$. For example, solar irradiance for the cardinally oriented cube would be calculated using $F_{\gamma }=[0.25,0.25,0.25,0.25]$ and $\gamma =[0^{\circ },{90}^{\circ },-{90}^{\circ },{180}^{\circ }]$.

The solar irradiance absorbed by the roof needs its own variable. The calculation procedure is the same as that for the vertical surfaces. However, there is no need for area weighting as the roof is assumed to be representable with one horizontal surface:

$$I_{tot;hor;sh;\mathrm{t}}=I_{dif;hor;\mathrm{t}}+I_{dir;hor;\mathrm{t}}·F_{sh;hor;\mathrm{t}}.$$

(30)

The solar irradiance at an arbitrarily tilted and oriented surface is calculated accordingly to the ISO 52010:2017 [25] standard. The calculation procedure was implemented in C++/R and validated against the accompanying test data of the the standard. The code with validation/unit tests is published under a public domain license and is located at the first author’s GitHub code repository (https://github.com/lukas-rokka/solarCalcISO52010).

3.4. Shading Reduction Factor

Distant obstacles (hills, trees, buildings etc) and obstacles on or nearby the building (balconies, rebates etc.) block parts of the solar radiation reaching a buildings surface. For best accuracy, these are calculated according to the actual surroundings of the building site. However, one of the goals with this work was to make reasonable estimates without actual information of the surroundings. Therefore, the shading factors are calculated based on assumptions about what the surroundings of a typical Swedish urban multifamily building typically look like. The resulting estimation can be expected to be more accurate on average than if shading were not accounted for at all. The shading reducing factor for the vertical envelope is calculated as

$$F_{sh;ver;\mathrm{t}}=({\alpha }_{sol;\mathrm{t}}>0)·(F_{sh;obst;\mathrm{t}}+F_{sh;ovh;\mathrm{t}}-1)$$

(31)

where $F_{sh;obst}$ is the shading from obstacles, and $F_{sh;ovh}$ is the shading from overhangs on the buildings. The shading factor for infinite length overhangs is calculated as described by Method 2 in Annex F of ISO 52016-1:2017 [10]:

$$F_{sh;ovh;\mathrm{t}}=1-max(0,min(1,(D_{ovh}·tan({\alpha }_{sol;\mathrm{t}})-L_{ovh})/H_{gl}))$$

(32)

where $D_{ovh}$ is the depth of the overhang, ${\alpha }_{sol}$ is the solar altitude, $L_{ovh}$ is the vertical distance between the glazing and the overhang, and $H_{gl}$ is the height of the window glazing. Many buildings have balconies; to account for this to at least some extent, the following values were used $D_{ovh}=2\mathrm{m}$, $L_{ovh}=1\mathrm{m}$, and $H_{gl}=6\mathrm{m}$. Setting $H_{gl}=6\mathrm{m}$ is approximately the same as having overhangs above every 4th window of a more typical window glazing height of 1.5 $\mathrm{m}$.

The shading reduction factor calculation for distant obstacles is calculated as a mix of two semitransparent obstacles:

$$\begin{array}{cc}\hfill F_{sh;obst;\mathrm{t}}=1-& min(H_b,max(0,H_{obst;1}-{\lambda }_{sol}·L_{obst;1}·tan({\alpha }_{sol;\mathrm{t}})))/H_b+\hfill \\ & min(H_b,max(0,H_{obst;2}-{\lambda }_{sol}·L_{obst;2}·tan({\alpha }_{sol;\mathrm{t}})))/H_b\hfill \end{array}$$

(33)

where $L_{obst}$ is the distance between the shading obstacle and the shaded object, $H_{obst}$ is the height of the obstacle, $H_b$ is the height of the shaded building, and ${\lambda }_{sol}$ is a shading factor set based on the categorization of the surroundings of the building (see Table 2). Obstacle 1 models nearby objects such as trees, and its parameters are set as $L_{obst;1}=15\mathrm{m}$ and $H_{obst;1}=20\mathrm{m}$. Obstacle 2 models shading from other buildings, and its parameters are set as a function of the building height $L_{obst;2}=H_b$ and $H_{obst;2}=1.3·H_b$. Shading at the roof level is calculated as

$$F_{sh;hor;\mathrm{t}}=min(3,H_b-H_{obst;1})/3·0.5+0.5.$$

(34)

3.5. Window Blinds

Window blinds exist in most Swedish residential buildings. However, there is no comprehensive consensus about the way people operate blinds or the motivating factors that influence their decisions [33]. A Swedish survey by Sandberg and Engvall [34] showed a clear correlation between exposure to solar radiation and to both the occurrence and position of window blinds. The g-value of the window blinds at time interval $\mathrm{t}$ is calculated as

$$g_{bl;\mathrm{t}}=1+(g_{bl;max}-1)·u_{bl;\mathrm{t}}$$

(35)

where $g_{bl;max}$ is a parameter describing the g-value of the blinds when at maximum blocking capacity, and $u_{bl;\mathrm{t}}$ is the control signal to the blinds. If the window blinds were automatically operated, the $u_{bl;\mathrm{t}}$ could for example be modeled as an on/off signal depending on solar radiation and a threshold value. The following function attempts to deterministically estimate the control signal of occupant operated window blinds:

$$\begin{array}{cc}\hfill u_{bl;\mathrm{t}}& =max(0.3,min(0.7,I_{tot;wma;\mathrm{t}}/200))\hfill \\ \hfill I_{tot;wma;\mathrm{t}}& =\sum _{k=0\mathrm{h}}^{23\mathrm{h}}\left(\frac{24-k}{300}·I_{tot;\mathrm{t}-k}\right)\hfill \end{array}$$

(36)

where the max and min functions limit the control signal between 30 and 70%, and $I_{tot;wma;\mathrm{t}}$ is a weighted moving average of solar radiation exposure from the last 24 h. The motivation behind using weighted moving average was to model the position of the blinds as depending both on the current (active control) and recent (reactive control) solar radiation exposure. The use of a 24 h moving average will also result in a seasonal effect, as there are more hours with no solar radiation in winter than in the summer (modeling occupants as slightly more aware of the positions of the blinds in the summer season seems to be based on a reasonable assumption). The motivation behind limiting the control signal between the somewhat arbitrary 30 and 70% was that it can be expected that some blinds are also drawn in a situation of no solar exposure, that it is unlikely that all windows have blinds, and/or that these are in full position in a situation of high solar exposure.

4. Data Input

Section 4.1 and Section 4.2 describe the environmental data sources proposed for the ISO14N modeling framework. For the comparison study, environmental data for Norrköping, Sweden ($58.575$ ${}^{\circ }$ N, $16.15$ ${}^{\circ }$ E) and the year 2016 was used. Section 4.3 describes the example model building as implemented in IDA ICE and ISO14N.

4.1. Solar Irradiance

Direct normal and diffuse horizontal irradiance data were acquired from CAMS [35] (http://www.soda-pro.com/web-services/radiation/cams-radiation-service). This satellite-based solar irradiance data are available at a horizontal resolution of 5 km and 15 min time steps from 2004 until the present time, and covers Europe, Africa, and the Middle East.

4.2. Meteorological Reanalysis Data

Table 3. Environmental variables derived from ERA5. The ISO14N column shows variable names as they are used in this paper, while the ERA5 column shows the variables sourced from ERA5 according to the notation convention of ERA5.

Description	ISO14N	ERA5	Transformation
External air temperature at 2 m	${\theta }_e$	$2t$	Kelvin to centigrade
Ground temperature at 1.0–2.89 m	${\theta }_{gr}$	$stl4$	Kelvin to centigrade
Wind speed at 10 m	$U_{10m}$	$10u$, $10v$	$U_{10m}=\sqrt{10u^2+10v^2}$
Ground albedo	${\alpha }_{gr}$	$ssr$, $ssrd$	${\alpha }_{gr}=1-ssr/ssrd$
Surface thermal radiation downwards	${\varphi }_{strd}$	$strd$	Joule to Watt-hours

shows which variables were acquired from the reanalysis climate dataset ERA5 of the Copernicus Climate Change Service. ERA5 data are available at a global horizontal resolution of 31 km and hourly time steps [36] (https://cds.climate.copernicus.eu/cdsapp). Reanalysis uses numerical weather prediction schemes to assimilate historical observational data.

Apparent sky temperature was compared with results from IDA ICE but was not otherwise used in the ISO14N modeling framework. It is calculated as

$${\theta }_{sky;\mathrm{t}}={({\varphi }_{strd}/\sigma )}^{0.25}-273.15.$$

(37)

4.3. Example Building Model Data Input

This section describes the example building model used for comparison. Section 4.3.2 describes the example building model as constructed in IDA ICE. Section 4.3.1 describes the parameter input needed to construct the same building model in ISO14N format. The IDA ICE model was originally [37] based on an aerated concrete, four-floor high, multifamily building, typical for the early Swedish million homes programme era of the 1960s and 1970s. The material properties of walls, windows, and the roof were kept according to the original model, but the shape was simplified to a 20 × 10 ${\mathrm{m}}^2$ rectangular, one-floor-high building. The simpler shape was chosen so that thermal zoning would not impact the results and to facilitate logging node temperatures in IDA ICE.

4.3.1. ISO14N

No shading was accounted for. Window blinds were assumed fully drawn. The supply temperature look-up table of

Table 4. Supply temperature look-up table, used for both IDA ICE and ISO14N models.

${\theta }_e$ [${}^{\circ }C$]	−20	0	12	18
${\theta }_{hyd;sup;set}$ [${}^{\circ }\mathrm{C}$]	60	43	28	18

was used. Used parameter inputs are given in

Table 5. Parameter input values used for the ISO14N model simulation.

Parameter	Description	Value(s)	Unit
$N_{fl}$	Number of floors	1	-
$H_{fl}$	Floor height	2.6	$\mathrm{m}$
P	Building perimeter	60	$\mathrm{m}$
$A_{fl}$	Floor area	2.6	${\mathrm{m}}_{\mathrm{fl}}^2$
$\gamma$	Azimuth angles for the external walls	[0, 90, −90, 180]	${}^{\circ }$
$F_{\gamma }$	Surface area fraction for the external walls	[0.17, 0.33, 0.33, 0.17]	-
$r_{gl}$	Glazing to floor ratio	0.15	-
$g_{gl}$	Solar transmittance of the glazing	0.76	-
$g_{bl}$	Solar transmittance of the window blinds	0.53	-
$U_{rf}$	Thermal transmittance of the roof	0.20	$\mathrm{W}/{(}^{\circ }\mathrm{C}·{\mathrm{m}}_{\mathrm{s}}^2)$
$U_{ew}$	Thermal transmittance of the external walls	0.72	$\mathrm{W}/{(}^{\circ }\mathrm{C}·{\mathrm{m}}_{\mathrm{s}}^2)$
$U_{gl}$	Thermal transmittance of the glazing	2.9	$\mathrm{W}/{(}^{\circ }\mathrm{C}·{\mathrm{m}}_{\mathrm{s}}^2)$
$U_{gf}$	Thermal transmittance of the ground floor	0.23	$\mathrm{W}/{(}^{\circ }\mathrm{C}·{\mathrm{m}}_{\mathrm{s}}^2)$
${\kappa }_{m;rf}$	Areal heat capacity of the external walls	56	$\mathrm{Wh}/{(}^{\circ }\mathrm{C}·{\mathrm{m}}_{\mathrm{s}}^2)$
${\kappa }_{m;ew}$	Areal heat capacity of the external walls	15	$\mathrm{Wh}/{(}^{\circ }\mathrm{C}·{\mathrm{m}}_{\mathrm{s}}^2)$
${\kappa }_{m;gf}$	Areal heat capacity of the external walls	56	$\mathrm{Wh}/{(}^{\circ }\mathrm{C}·{\mathrm{m}}_{\mathrm{s}}^2)$
$Q_{ve}$	Specific air flow rate of the ventilation system	0.35	$\mathrm{l}/(\mathrm{s}·{\mathrm{m}}_{\mathrm{fl}}^2)$
${\eta }_{ve}$	Efficiency of the ventilation heat recovery	0.0	-
$C_{inf}$	Normalized infiltration coefficient	0.0	$\mathrm{l}/(\mathrm{s}·{\mathrm{Pa}}^{\mathrm{n}}·{\mathrm{m}}_{\mathrm{fl}}^2)$
$H_{hyd}$	Radiator constant, normalized per floor area	0.66	$\mathrm{W}/{(}^{\circ }\mathrm{C}·{\mathrm{m}}_{\mathrm{fl}}^2)$
${\theta }_{int;set}$	Internal air set-point temperature	21.0	${}^{\circ }\mathrm{C}$
$H_{tb}$	Heat transfer due to thermal bridges	0.0	$\mathrm{W}/{(}^{\circ }\mathrm{C}·{\mathrm{m}}_{\mathrm{fl}}^2)$
${\varphi }_{int}$	Internal heat gains	3.0	$\mathrm{W}/({\mathrm{m}}_{\mathrm{fl}}^2)$
$\Delta t$	Calculations interval	0.5

.

4.3.2. IDA ICE

The model was constructed using IDA ICE version 4.8. The modeled building is a simple 20 × 10 ${\mathrm{m}}^2$ rectangular, one-floor-high building, with the long sides oriented in an east–west direction and a window-to-floor ratio of 0.15. The following constructions were used (layers from outside to inside):

• Roof: 0.20 m insulation + 0.15 m concrete, U-value of 0.20 $\mathrm{W}/({\mathrm{m}}_{\mathrm{s}}^2·\mathrm{K})$;
• External walls: 0.20 m aerated concrete, U-value of 0.72 $\mathrm{W}/({\mathrm{m}}_{\mathrm{s}}^2·\mathrm{K})$;
• Glazing: 2-pane, U-value of 2.9 $\mathrm{W}/({\mathrm{m}}_{\mathrm{s}}^2·\mathrm{K})$, g-value of 0.76;
• Ground floor: 0.1 m virtual ground layer + 0.5 m soil + 0.1 m insulation + 0.1 m concrete, U-value of 0.22 $\mathrm{W}/({\mathrm{m}}_{\mathrm{s}}^2·\mathrm{K})$(0.33 if not including the ground);
• Internal walls: 150 ${\mathrm{m}}^2$ of 0.15 m thick aerated concrete;
• Furniture: 130 ${\mathrm{m}}^2$ of 0.01 m thick default furniture material.

The heating system was modeled using the IDA ICE Water Radiator model (CeWatHet) with a nominal heat output of 50 $\mathrm{W}/{\mathrm{m}}_{\mathrm{fl}}^2$ (at a supply temperature of 60 ${}^{\circ }\mathrm{C}$ and a return at 40 ${}^{\circ }\mathrm{C}$), and the geometry of the radiator was set so that approximately half of the heat output was convective. The thermostatic control was modeled with a P-controller with a proportional band of 2 ${}^{\circ }\mathrm{C}$, a time constant of 20 min, and the set-point set to 21 ${}^{\circ }\mathrm{C}$. Table 4 shows the used look-up table for supply temperatures to the radiators. Ventilation was set to a constant 0.35 $\mathrm{l}/(\mathrm{s}·{\mathrm{m}}_{\mathrm{fl}}^2)$ exhaust flow. Window blinds g-value was set to 0.53 and always fully drawn, internal heat gain was set to a constant of 3 $\mathrm{W}/{\mathrm{m}}_{\mathrm{fl}}^2$, and thermal bridges were set to zero. Infiltration was nominally set to zero; however, for the case of studying the air infiltration (Section 5.4), the building envelope airtightness was set to 1 ACH at a 50 Pa pressure difference, and the surface pressure coefficient was semi-exposed.

5. Results

The developed ISO14N model and its procedures were validated by comparing it with simulations conducted with the detailed building energy simulation software IDA ICE.

5.1. Full-Year Simulation Comparison with IDA ICE

In this section, simulation results from the developed ISO14N model are compared with results from an IDA ICE simulation. As seen in Figure 3, the ISO14N is capable of reproducing results of the detailed IDA ICE model (see Appendix B for comparisons on monthly and daily scales). The first half of January 2016 had two successive cold spells. During this period, the models show larger deviations, suggesting that they differ in thermal capacity behavior in the lower frequencies. A perfectly sized heating system would result in internal air temperature close to the set-point during the heating season. The heating system was on purpose undersized to cause more deviations from the set-point temperature, which makes differences in dynamical behavior more apparent.

Table 6. Root mean square deviation (RMSD) for hourly, daily, and monthly averages on energy use (${\varphi }_{hyd}$) and internal air (${\theta }_{int}$) and operative temperatures (${\theta }_{op}$) variables. For ${\varphi }_{hyd}$, CV(RMSD) is shown within parenthesis.

Variable	Hourly	Daily	Monthly
${\varphi }_{hyd}$, $\mathrm{W}/{\mathrm{m}}_{\mathrm{fl}}^2$	0.94 (2.6%)	0.74 (2.0%)	0.68 (1.8%)
${\theta }_{int}$, ${}^{\circ }\mathrm{C}$	0.30	0.26	0.17
${\theta }_{op}$, ${}^{\circ }\mathrm{C}$	0.34	0.28	0.19

shows the root mean square deviation (RMSD) [38] and the coefficient of variation (CV) of the RMSD [38], calculated on hourly, daily, and monthly averages. As can be seen from the table, the ISO14N is able to replicate simulation results of IDA ICE on an hourly basis, which is much more sensitive to dynamics than daily or monthly averages.

5.2. Node Temperature Profiles of External Wall Elements

For the developed ISO14N model, the external walls are represented by one element with three temperature nodes, while the IDA ICE model has four external wall elements (one in each cardinal direction) with five temperatures nodes each. For better comparability, the external wall temperatures of the IDA ICE model is averaged by node position. As seen in Figure 4, the node temperature profiles of the two models behave quite differently, which is mostly due to the different number of nodes and different approaches to distribute thermal properties between the nodes. Still, the estimated internal air temperatures match well.

5.3. Operative Temperature

Operative temperature was calculated in two different ways and compared with output from IDA ICE. “ISO14N alt. 1” in Figure 5 was calculated according to Equation (18), while “ISO14N alt. 2” was calculated using Equation (19) with the weighting/conversion factor $C_{hyd;r;mn}$ (expressing the impact of the hydronic heating system) set to $C_{hyd;r;mn}=0.02$. As seen in Figure 5, the operative temperature of “ISO14N alt. 1” deviates from the operative temperature acquired from IDA ICE, especially during cold weather, which was expected as IDA ICE weights in the surface temperatures of the radiator panels in the calculation of operative temperatures. Using Equation (19) and setting $C_{hyd;r;mn}=0.02$ results in a closer match. The $C_{hyd;r;mn}$ parameter could be derived from typical radiator panel configurations, but here it has been chosen by empirically matching it to the results of IDA ICE.

From Figure 5, it can be observed that IDA ICE appears to be more responsive to higher frequencies. This can partly be explained by the fact that IDA ICE splits the building elements into more nodes and uses surface nodes without mass. The other part of the explanation is that IDA ICE uses a fully dynamic simulation procedure (dynamic state transitions and variable heat transfer coefficients), while ISO14N linearizes the state transitions and uses constants for the interior surface heat transfer coefficients.

5.4. Air Infiltration

The air infiltration implementation described in Section 2.4 (with shelter class set to 3 and the flow coefficient set to $C_{inf}=0.08\mathrm{l}/(\mathrm{s}·{\mathrm{Pa}}^{\mathrm{n}}·{\mathrm{m}}_{\mathrm{fl}}^2)$ is here compared with output from IDA ICE. As seen in Figure 6, the implemented infiltration model produces similar results as the single-zone, one-floor-high, IDA ICE building model.

5.5. Sky Temperature and Thermal Radiation to the Sky

Apparent sky temperature was calculated from the ERA5 variable “surface thermal radiation downwards” using Equation (37) and compared to output from IDA ICE. In IDA ICE, the apparent sky temperature was empirically derived from dry and dew point temperatures and cloud cover using the Walton–Clark–Allen model [39]. As seen in Figure 7, the proposed apparent sky temperature calculation differs substantially from that of IDA ICE, resulting in an hourly RMSD of $5.5$ ${}^{\circ }\mathrm{C}$. Using the defaults suggested in the ISO 52016-1:2017 [10] standard, the extra thermal radiation to the sky variable ${\varphi }_{sky}$ results in a constant of 45 $\mathrm{W}/{\mathrm{m}}_{\mathrm{s}}^2$; deriving ${\varphi }_{sky}$ from the “surface thermal radiation downwards” variable of ERA5 (as proposed in Section 2.5), this constant varies between −15 and 295 $\mathrm{W}/{\mathrm{m}}_{\mathrm{s}}^2$, with an average of 61 $\mathrm{W}/{\mathrm{m}}_{\mathrm{s}}^2$ for the used climate file.

5.6. Exterior Surface Convective Heat Transfer Coefficients

Exterior surface convective heat transfer coefficients calculated according to Section 3.1 are presented and compared with IDA ICE. IDA ICE refers to Clarke [40] for its calculation procedure, a method that accounts for wind direction but does not distinguish between building element type. As seen in Figure 8, our Liu et al. [32]-based method results in the roof element being more sensitive to wind speed than other elements. The lower values for the glazing element, compared to the external walls, is explained by its surface roughness multiplier being set to smooth. IDA ICE’s method shows much less sensitivity to wind speed. For IDA ICE, the effect of wind direction is removed by showing smoothed lines of the average of the four external walls.

5.7. Shading Reduction Factor

Results of the proposed simplified shading calculation procedure (see Section 3.4) based on classification of the surrounding environment are presented.

Table 7. Annual reduction of total solar irradiance on vertical surfaces, per sheltering class and building height.

Shading Class	1	1	2	2	3	3	4	4	5	5
Building height, m	18	6	18	6	18	6	18	6	18	6
Annual solar reduction, %	6	9	11	17	19	28	29	40	38	49

shows the impact on annual solar irradiance received on vertical surfaces. The climate file described in Section 4 was used. The results can be compared to study conducted by Romero Rodríguez et al. [41], which showed an average of 35% reduction on solar irradiance for facades in the urban parts of Ludwigsburg. This part of the proposed modeling framework would benefit from further development, for example by employing aerial LiDAR data to model shading from the surrounding environment more accurately, see for example Lingfors et al. [17].

5.8. Window Blinds

Calculating the impact from window blinds according to Section 3.5, with the window blinds maximum g-value set to $g_{bl;max}=0.4$, results in a 30% reduction of annual solar heat gains through window glazing. According to Van Den Wymelenberg [33], there is no comprehensive consensus about how to model occupant-operated window blinds. Nonetheless, the proposed procedure is anticipated to model the impact from window blinds more accurately than for example using a constant value. IDA ICE uses an on/off controller, which controls the window blinds position as a function of the outside solar radiation level (100 W/m² default set-point). This control scheme resulted in a 31% reduction of annual solar heat gains for the studied building and the climate file.

5.9. Computation Benchmark

This section show computation timings from conducted speed tests with the proposed ISO14N model, variants of it, and IDA ICE. The proposed ISO14N model pre-inversed its matrix $\mathbf{A}$ (the RC network) before the actual simulation run. For comparison, timing comparisons for model variants, where the matrix $\mathbf{A}$ of the RC network was inverted at each time-step, were also conducted. To compare speed in relation to a full ISO 52016-1:2017 implementation, an emulation was constructed and denoted “ISO41N” (due to its 41 temperature nodes).

Table 8. Computation benchmark, timing in milliseconds.

Model (Variant)	Computation Time, ms
ISO14N (matrix pre-inverted)	19
ISO41N (matrix pre-inverted)	56
ISO14N (matrix inversion each time-step)	83
ISO41N (matrix inversion each time-step)	755
IDA ICE	18,000

shows the results. The largest speed gain was achieved from the matrix pre-inversion. Matrix inversion is computationally intensive and has a computation cost relationship of O($n^2$). The use of a pre-inverted matrix $\mathbf{A}$ breaks the O($n^2$) computational cost relationship and results in an approximately linear relationship between the RC network size and computation time. For the pre-inverted cases, it can be seen that the proposed simplification compared to a full ISO 52016-1:2017 implementation results in an approximately three-fold speed gain. The computation time difference between IDA ICE and the ISO14N is approximately 900 times (or an order of 1 magnitude).

The benchmark was conducted on a standard laptop computer with an Intel Core i7-7500U processor using one thread. All models were simulated using one year of climate, plus 14 days for the adaptation period. IDA ICE simulation was performed within the user interface, using initial default settings and all variable logging turned off; timing was conducted manually with a stopwatch. Various types of the ISO14N model used a calculation interval of 0.5 h (18,384 calculation time-steps in total), and logged and returned energy use and internal air temperature.

6. Discussion

It has been shown that the proposed ISO14N model was able to reproduce the hourly energy use pattern obtained with the detailed building energy simulation software IDA ICE. The use of a hydronic radiator heating system limits the application to mainly multifamily buildings in Nordic climates. However, the proposed RC network itself is generic and could probably replace a full ISO 52016-1:2017 implementation in other applications, such as generic energy need calculations or as a basis for energy use calculation for cooling.

In the development of the proposed ISO14N model, effort was made to achieve fast simulation times without a significant loss of accuracy. The optimization for computational efficiency is motivated by that calibration approaches typically use algorithmic differentiation which requires the model to be simulated several times (10k or more for Monte Carlo simulation approaches). Speed gains were achieved by moving parts of the calculation load to the pre-processing, i.e., model nodal reduction, and by a simulation procedure that uses pre-inverted matrices for the RC network. The proposed simplification of the RC network consists of lumping building elements by type and decreasing the number of temperature nodes in opaque building elements (see Figure 2), resulting in a system of 14 temperature nodes. The largest speed gain was achieved from pre-inverting the matrix $\mathbf{A}$ of the RC network, breaking the O($n^2$) cost relationship and resulting in an approximately linear relationship between the RC network size and computation time. Pre-inversion requires the use of a constant matrix $\mathbf{A}$ element values, which results in a less flexible modeling structure. A constant $\mathbf{A}$ matrix, in turn, requires the use of node temperatures from the earlier time step for certain calculations (e.g., for the heating system control), which can result in time-shifts and numerical instabilities. However, our results show that the proposed ISO14N model achieves stable results without noticeable time-shifts.

During the development stage, a five-temperature node version (not included in this paper) was also tested and showed only limited performance gains compared to the proposed three-node version. Using five nodes adds flexibility and is likely to better model thermal admittance and more complex building structures such as sandwich walls (Akander [42]). However, for the intended use case of auto-generation of baseline models in large building stocks, it is anticipated that detailed information about construction will seldom be available in digital form. In the case of access to more detailed information, the proposed thermal network can be extended by adding more temperature nodes and/or by dividing the modeled building into more thermal zones.

Procedures to acquire and make use of satellite-based solar radiation data were presented, see Section 3.3. Solar radiation has a large impact on the thermal balance of buildings, and thermostatic radiators valves are often specifically used to avoid overheating due to solar heat gains. Thus, being able to model solar heat gains on an hourly basis (both directly through glazing and transmitted through opaque building elements) is anticipated to facilitate parameter estimation with hybrid/gray-box calibration methods. Reanalysis climate model data from ERA5 (see Section 4.2) was utilized. Conventional climate variables wind and air temperature were used as well as the more unconventional variables of ground temperature, ground albedo, and surface thermal radiation downward. Some of the benefits of using reanalysis data sources such as ERA5 are that the data are homogeneous and harmonized, covering the entire world, there is no need to deal with missing data, there is access to variables that are seldom measured locally, and the data are readily and publicly accessible. The temporal and spatial resolution of ERA5 is relatively rough, and in some cases it might be better to use local measurements (e.g., air temperature is readily measured locally at the building site).

7. Conclusions

A space heating model suitable for auto-generating baseline models of existing multifamily buildings has been proposed. As demonstrated in the result section, the proposed ISO14N model was able to reproduce the hourly energy use of space heating, indoor temperature, and operative temperature patterns obtained from the detailed building energy simulation software IDA ICE. Thus, the proposed model can be expected to model existing multifamily buildings, where space heating is provided by hydronic heating systems, to a similar degree of confidence as established simulation software.

The conducted computation timing showed that the proposed ISO14N was in the approximately 900 times faster than the detailed IDA ICE simulation. The largest speed gain was achieved from pre-inverting the matrix $\mathbf{A}$ of the RC network, resulting in an approximately linear cost relationship between the RC network size and computation time. Despite the use of constant pre-inverted $\mathbf{A}$ matrix, the proposed ISO14N model can achieve stable results without noticeable time-shifts. The achieved fast computation time enables using time-sensitive applications, such as Monte Carlo based calibration methods.

Solar heat gains have large impact on the thermal balance of buildings, and hydronic heating system often use thermostatic radiators valves to avoid overheating due to solar heat gains. Thus, it is necessary to model actual solar heat gains. We have presented procedures for converting sub-hourly satellite-based solar irradiation data to solar heat gains, while considering effects from building orientation, shading, and window blinds.

Future work will be conducted on calibrating the proposed model with actual energy meter readings, to ensure it represents the actual performance of the building and not only intended/designed performance. Hourly meter readings from district heating substations usually do not distinguish between domestic hot water and space heating. Thus, energy use for the occupant dependent domestic hot water needs to be modeled to be able to apply the proposed model on real buildings. Also, internal heat gains should be either measured (domestic electricity is always measured but not always readily available) or modeled to account for its occupant behavior dependent variations. The solar heat gain calculation could also be improved by employing aerial LiDAR data to model shading from the surrounding environment more accurately.