Laplace and State-Space Methods for Calculating the Heat Losses in Case of Heavyweight Building Elements and Short Sampling Times

Abstract

Reducing heat losses through the building envelope is one of the most important aspects to be met if the targets set by the European Union are to be achieved. In order to obtain a more realistic energy demand, dynamic heat transfer simulations are used to calculate the energy consumption of buildings, since steady-state calculations do not take into account the thermal mass in buildings. These dynamic simulations employ methods based on analytical models since numerical models are unsuitable for longer time periods. The analytical models used herein fall into the category of conduction transfer functions (CTFs). Two methods for computing CTFs that are addressed in this research are the Laplace method and the State-Space method. The objective of this paper is to verify the efficiency of the Laplace and State-Space methods for calculating the energy demand of a building in the case of heavyweight building elements and shorter sampling times, and to provide a means for improving the algorithms used by these methods. The Laplace and State-Space method algorithms were implemented in Mathematica, and the results were compared to EnergyPlus and TRNSYS, which use similar algorithms to calculate energy demand. It was shown in this paper that for the heavyweight wall element and a time step of 0.25 h, the difference between the total energy transferred through the inner surface was about 31% for EnergyPlus and 78% for TRNSYS compared to the reference solution. For the lightweight wall element, the results were stable for the time step of 0.25 h, but for the time step of 0.1 h, the differences were 45.64% and 303% between EnergyPlus, TRNSYS and the reference solution, respectively, compared to the State-Space and Laplace methods for which the maximum difference was 12.03% with a time step of 0.1 h. While dynamic heat transfer simulations are better than calculations based on steady-state boundary conditions, they also have their limitations and could lead to unsatisfactory results for short sampling times and if not applied properly.

1. Introduction

With the aim of ensuring that EU policies are in line with the climate goals agreed by the Council and the European Parliament, the Fit for 55 [1] package was proposed, which among other things, includes an increase in the CO_2eq emissions from the buildings reduction target by 2030 compared to 2005. The measures include, among others, enhancing buildings’ energy performance when constructing new buildings and the energy renovation of existing buildings. All of this also complies with the European Commission’s long-term strategy for a low-carbon Europe by 2050 [2,3]. The European parliament adopted draft measures to increase the building renovation rate as well as to minimize buildings’ energy consumption and CO_2eq emissions. European Climate Law [4] enshrined both the 2030 and the 2050 targets into binding European law. According to the adopted measures, all new buildings should be zero-emission buildings (ZEBs) from 2028, while new buildings occupied, operated or owned by public authorities should be from 2026 [5], thus advancing even further the development set up by nearly zero energy buildings (NZEBs).

To achieve these goals, it is important to predict the energy consumption of both new buildings and renovated buildings. Most software used in the EU for energy performance simulation use simplified methods to calculate heat transfer through building elements based on steady-state boundary conditions averaged over a period of time [6,7]. This calculation can lead to unsatisfactory results when boundary conditions change rapidly and when the thermal mass of the element is significant [6,8]. As a result, the predicted energy demand may be lower or higher than the actual energy consumption, leading to undersized or oversized heating and cooling systems. Furthermore, if the boundary conditions are considered to be steady-state, the peak energy demand shifts because the thermal mass is not taken into account, leading to a thermal comfort problem for the building occupants. Nowadays, dynamic methods for simulating building energy demand are more and more widely used because they are powerful and important tools for optimizing building energy demand [8,9,10,11,12].

The dynamic methods are based on the non-stationary solution of the heat equation. This solution can be found numerically (finite volume method, finite element method or finite difference method) or analytically. Due to the long analysis period and small time steps used for the analysis, numerical methods have proven to be inadequate for building energy simulations because they rely on solving large systems of equations at each time step [13]. Analytical methods such as conduction transfer functions (CTFs) have been shown to be much better for building energy simulation problems [14] and are implemented in modern dynamic energy simulation software such as TRNSYS [15] and EnergyPlus [16]. Although the CTF methods are superior to the abovementioned simplified calculations, the aim of this research is to show that they have their limitations and lead to unsatisfactory results if not properly applied. It will be shown that in the case of heavyweight façade systems and small computational time steps, CTFs can provide low-quality results (become more and more unstable with longer time steps, particularly with massive building elements) [17]. The reason for using small time steps, which in some cases can be less than one hour, could be to analyze more accurately the impact of complex control systems on the overall energy performance of the building primarily associated with the exploitation of renewable energy sources as well as predictive model control of buildings [18,19]. Another reason is that NZEBs and ZEBs require more accurate dynamic simulations of heavy and highly insulated building elements with shorter time steps (less than 15 min) [20,21,22]. The problem also lies in the algorithms developed for yearly energy simulations with a time step of one hour. Comparisons of different methods and their main problems found in the literature are summarized in

Table 1. Comparison of different methods from the literature.

Reference	Method	Time Step Size	Main Problem
[8]	EN ISO 52016-1 hourly method	1 h	Hourly assumption leads to non-negligible inaccuracy in the calculated energy consumption
[6]	Finite element method (FEM)	Not specified	Higher computational time
[9]	GIS technology	Not specified	Too many specific data are needed and too much equipment per building
[10]	Comparison of different simulation software	Not specified	Large computer memory and processing or computation time for more accurate results of building energy demand in real time
[11]	CTF	0.5–1 h	Termination due to CTF calculation instability
[15]	CTF	0.25 and 1 h	Instability for smaller step sizes
[16]	Recurrent neural network via the Multi-Input Multi-Output strategy	24 and 168 h	Gaps between predicted and actual energy use at the peak and base load times are significant
[17]	CTF	0.02–1 h	Problems for short-term analysis for heavyweight building elements
[18]	DesignBuilder/EnergyPlus	0.25–1 h	Huge computation time, complexity of the model implementation
[19]	CTF	0.1–1 h	Using a 60 min time step could underestimate the peak load up to 70%

.

Among the most commonly used CTF methods are the Laplace method and the State-Space method. This paper compares these two methods for calculating heat losses through a heavyweight multilayer wall with a contact thermal insulation façade system by implementing the algorithms from the research work of Stephenson and Mitalas [23] for the Laplace method and Seem [24] for the State-Space method in Mathematica [25]. The software Mathematica was chosen due to its advanced mathematical and numerical tools and matrix manipulation. The results from Mathematica are then compared with the results of two commonly used software for the dynamic calculation of energy performance—TRNSYS and EnergyPlus. The first software uses the Laplace method [21,26] to calculate the CTF coefficients, while the second uses the State-Space method [27,28,29]. The finite difference method (FDM) [30] is used as a reference point for the boundary heat flux density. The aim is to show how well the Laplace method and the State-Space method perform for small time steps and heavyweight building elements and what possibilities there are for improving these methods.

2. Methodology

If we consider a multilayer wall, as shown in Figure 1, in the general case, the wall is composed of n layers and two boundary conditions. The boundary condition for x = 0 is T_int and for x = L T_ext. For the building energy performance simulation, the surface heat fluxes on the interior (q_int) and the exterior (q_ext) surface must be calculated since these two fluxes represent the energy demand for the element under consideration. Only heat flux due to conduction is considered and boundary conditions are surface temperatures on the internal and external surfaces. Convection and radiation on both surfaces are not considered, i.e., R_si and R_se are taken as 0 (m² K)/W. Solar radiation as well as all other radiation from the environment is not taken into account for the calculation. This assumption simplifies the calculation of the surface heat fluxes and the implementation of the Steady-State and Laplace method algorithms in Mathematica. The same assumption is implemented in both Mathematica’s algorithms, EnergyPlus and TRNSYS.

In most cases, heat transfer through building elements such as walls or ceilings can be considered one-dimensional. The heat equation for one-dimensional heat conduction is given by Equation (1).

$$\frac{\partial T(x,t)}{\partial t}=\frac{\lambda }{\rho \cdot c_p}\cdot \frac{{\partial }^{2}T(x,t)}{\partial x^2},$$

(1)

With the boundary and initial conditions:

$$\begin{array}{l}T(0,t)=T_{\mathrm{int}}\\ q(0,t)=q_{\mathrm{int}}\\ \\ T(L,t)=T_{ext},\\ q(L,t)=q_{ext}\\ \\ T(x,0)=0\end{array}$$

(2)

where T(x,t) is the temperature at point x at time t in °C, λ is the thermal conductivity in W/(m K), $\rho$ is the density in kg/m³, c_p is the specific heat in J/(kg K), T_int is the internal surface temperature in °C, T_ext is the external surface temperature in °C, q_int is the internal surface heat flux in W/m² and q_ext is the external surface heat flux in W/m².

In this research, Equation (1) is solved using:

• The State-Space method [24];
• The Laplace method [23];
• The finite difference method;
• The State-Space method, as implemented in EnergyPlus;
• The Laplace method, as implemented in TRNSYS.

For the solution of each method, the internal surface heat flux is calculated. The solution for all five methods is then compared to analyse the impact of each method on the solution, and also to analyse the efficiency of each algorithm used by the abovementioned methods. This analysis is performed for time steps ranging from one hour to six minutes and for lightweight and heavyweight building elements (walls). As a reference solution for the comparison (i.e., as the exact solution), the solution of the FDM is taken. Furthermore, the steady-state solution is calculated to show the difference between the transient (dynamic) and steady-state building energy demand simulation.

2.1. The Finite Difference Method

The finite difference method (FDM) is a method used for obtaining a numerical solution to Equation (1). As with most numerical solutions, the FDM uses a discrete approximation of the continuous partial differential equation. This means that the solution is known only at a finite number of points within the element. The number of these points is user-defined and generally increases the accuracy of the numerical solution. The discrete approximation of Equation (1) is as follows [31]:

$$\frac{T_i^{j+1}-T_i^j}{\mathsf{\Delta }t}=\frac{1}{\mathsf{\Delta }{x}^2}\frac{{\lambda }_1\cdot T_{i-1}^j-\left({\lambda }_{i-i}+{\lambda }_i\right)\cdot T_i^j+{\lambda }_2\cdot T_{i+1}^j}{\left({\rho }_{i-1}\cdot c_{i-1}+{\rho }_i\cdot c_i\right)/2},$$

(3)

Equation (3) represents the implicit FDM formulation for solving the heat equation. In the matrix formulation, Equation (3) can be written as follows:

$$\left(\begin{array}{ccccc}1& 0& 0& \dots & 0\\ -r_1& 1+r_1+r_2& -r_2& \dots & 0\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& \dots & -r_{i-1}& 1+r_{i-1}+r_i& -r_i\\ 0& \cdots & 0& 0& 1\end{array}\right)\cdot {\left(\begin{array}{c}{T}_1\\ T_2\\ ⋮\\ T_{i-1}\\ T_i\end{array}\right)}_{j+1}={\left(\begin{array}{c}{T}_{\mathrm{int}}\\ 0\\ ⋮\\ 0\\ T_{ext}\end{array}\right)}_j,$$

(4)

The subscript i denotes a node and j point in time. The coefficients r₁ and r₂ are calculated for each node according to Equation (5).

$$\begin{array}{l}{r}_{i-1}=\frac{\mathsf{\Delta }t}{\mathsf{\Delta }{x}^2}\cdot \frac{{\lambda }_{i-1}}{\left({\rho }_{i-1}\cdot c_{i-1}+{\rho }_i\cdot c_i\right)/2}\\ r_i=\frac{\mathsf{\Delta }t}{\mathsf{\Delta }{x}^2}\cdot \frac{{\lambda }_i}{\left({\rho }_{i-1}\cdot c_{i-1}+{\rho }_i\cdot c_i\right)/2}\end{array},$$

(5)

By solving Equation (4) for each time step j, the internal heat flux q_int can be calculated using Fourier’s law of thermal conduction as follows:

$$q_{\mathrm{int}}={\lambda }_1\cdot \frac{T_2^j-T_{\mathrm{int}}^j}{\mathsf{\Delta }x},$$

(6)

The solution calculated using Equation (6) is used as a reference (i.e., as an exact solution) for the comparison.

2.2. The State-Space Method

The State-Space method is usually used to describe a linear system with r inputs, m outputs and n state variables that are written in the general form as described in [32]:

$$\dot{x}(t)=A(t)\cdot x(t)+B(t)\cdot u(t),$$

(7)

$$y(t)=C(t)\cdot x(t)+D(t)\cdot u(t),$$

(8)

Here x is the state vector, u is the input vector, y is the output vector, t is time, and A, B, C and D are matrices of dimension n × n, n × r, m × n and m × r, respectively. Since, in the case of heat conduction, the system is time-variant (i.e., the temperature inside the element changes with time), Equations (7) and (8) can be written by discrete time-variant State-Space method as in [32]:

$$\dot{x}(k+1)=A(k)\cdot x(k)+B(k)\cdot u(k),$$

(9)

$$y(k)=C(k)\cdot x(k)+D(k)\cdot u(k),$$

(10)

Here the matrices A, B, C and D have the same dimensions as in the continuous time case—Equations (7) and (8). In the case of heat conduction, the state vector x represents the temperatures of each node of the element, the input vector u represents the boundary conditions (i.e., the boundary temperatures shown in Figure 1) and the output vector y represents the surface heat fluxes q_i and q_e shown in Figure 1.

Equations (9) and (10) can be used to solve Equation (1) by enforcing the finite difference grid over the analysed building element:

$$\left(\begin{array}{c}\partial T_2/\partial t\\ \partial T_3/\partial t\\ ⋮\\ \partial T_{n-2}/\partial t\\ \partial T_{n-1}/\partial t\end{array}\right)=\left(\begin{array}{ccccc}-r_2-r_3& r_3& 0& \dots & 0\\ r_3& -r_3-r_4& r_4& \dots & 0\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& \dots & r_{n-3}& -r_{n-3}-r_{n-2}& r_{n-2}\\ 0& \cdots & 0& r_{n-2}& -r_{n-2}-r_{n-1}\end{array}\right)\cdot \left(\begin{array}{c}{T}_2\\ T_3\\ ⋮\\ T_{n-2}\\ T_{n-1}\end{array}\right)+\left(\begin{array}{ccc}{r}_1& \dots & 0\\ ⋮& \ddots & ⋮\\ 0& \cdots & r_n\end{array}\right)\cdot \left(\begin{array}{c}{T}_{\mathrm{int}}\\ ⋮\\ T_{ext}\end{array}\right)$$

(11)

The coefficients r₁ and r₂, which represent the left node “i − 1” and the node “i” for each row of the matrix, are calculated for each node using Equation (12).

$$\begin{array}{l}{r}_{i-1}=\frac{1}{\mathsf{\Delta }{x}^2}\cdot \frac{{\lambda }_{i-1}}{\left({\rho }_{i-1}\cdot c_{i-1}+{\rho }_i\cdot c_i\right)/2}\\ r_i=\frac{1}{\mathsf{\Delta }{x}^2}\cdot \frac{{\lambda }_i}{\left({\rho }_{i-1}\cdot c_{i-1}+{\rho }_i\cdot c_i\right)/2}\end{array},$$

(12)

where index “i” denotes the node in the FDM grid (i = 2 to n) and “i − 1” first previous node in the FDM grid.

The interior heat flux can be written in the matrix form as follows:

$$\left(q_{\mathrm{int}}\right)=\left(\begin{array}{c}-{\lambda }_1/\mathsf{\Delta }x\\ ⋮\\ 0\end{array}\right)\cdot \left(\begin{array}{c}{T}_2\\ ⋮\\ T_{i-1}\end{array}\right)+\left[\begin{array}{c}{\lambda }_1/\mathsf{\Delta }x\\ 0\end{array}\right]\cdot \left(\begin{array}{c}{T}_{\mathrm{int}}\\ T_{ext}\end{array}\right),$$

(13)

If Equations (11) and (13) are compared to Equations (9) and (10), then matrices A, B, C and D are equal to:

$$\begin{array}{l}A=\left(\begin{array}{ccccc}-r_2-r_3& r_3& 0& \dots & 0\\ r_3& -r_3-r_4& r_4& \dots & 0\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& \dots & r_{i-3}& -r_{i-3}-r_{i-2}& r_{i-2}\\ 0& \cdots & 0& r_{i-2}& -r_{i-2}-r_{i-1}\end{array}\right)\\ B=\left(\begin{array}{ccc}{r}_1& \dots & 0\\ ⋮& \ddots & ⋮\\ 0& \cdots & r_i\end{array}\right)\\ C=\left(\begin{array}{c}-{\lambda }_1/\mathsf{\Delta }x\\ ⋮\\ 0\end{array}\right)\\ D=\left[\begin{array}{c}{\lambda }_1/\mathsf{\Delta }x\\ 0\end{array}\right]\end{array},$$

(14)

The matrices (14) are the input values for the algorithm [24] to calculate the CTF coefficients Z_j, Y_j and φ_j. The algorithm [24] is implemented in Mathematica [25], and all operations with matrices were performed using this powerful software. From the CTF coefficients, the interior heat flux is calculated as follows:

$$q_{\mathrm{int}}(k\cdot \mathsf{\Delta }t)={\displaystyle \sum _{j=0}^{N_Z}{Z}_j\cdot T_{\mathrm{int},k\cdot \mathsf{\Delta }t-j\cdot \mathsf{\Delta }t}}+{\displaystyle \sum _{j=0}^{N_Y}{Y}_j\cdot T_{ext,k\cdot \mathsf{\Delta }t-j\cdot \mathsf{\Delta }t}}-{\displaystyle \sum _{j=1}^{N_{\mathsf{\Phi }}}{\mathsf{\Phi }}_j\cdot q_{\mathrm{int},k\cdot \mathsf{\Delta }t-j\cdot \mathsf{\Delta }t}},$$

(15)

where Z_j is the interior CTF coefficient, Y_j is the cross CTF coefficient and φ_j is the flux CTF coefficient. N_Z, N_Y and N_φ represent the number of Z_j, Y_j and φ_j coefficients, respectively.

It is common practice to use Equation (16) to verify that the CTF coefficients obtained by the State-Space method give the correct steady-state heat transfer, with X_k being the exterior CTF coefficient.

$$\frac{{\displaystyle \sum _{\mathrm{k}=1}^{{\mathrm{N}}_{\mathrm{a}}}{\mathrm{X}}_{\mathrm{k}}}}{1-{\displaystyle \sum _{\mathrm{k}=1}^{{\mathrm{N}}_{\mathrm{k}}}{\mathsf{\Phi }}_{\mathrm{k}}}}=\frac{{\displaystyle \sum _{\mathrm{k}=1}^{{\mathrm{N}}_{\mathrm{b}}}{\mathrm{Y}}_{\mathrm{k}}}}{1-{\displaystyle \sum _{\mathrm{k}=1}^{{\mathrm{N}}_{\mathrm{k}}}{\mathsf{\Phi }}_{\mathrm{k}}}}=\frac{{\displaystyle \sum _{\mathrm{k}=1}^{{\mathrm{N}}_{\mathrm{b}}}{\mathrm{Z}}_{\mathrm{k}}}}{1-{\displaystyle \sum _{\mathrm{k}=1}^{{\mathrm{N}}_{\mathrm{k}}}{\mathsf{\Phi }}_{\mathrm{k}}}}\approx \mathrm{U},$$

(16)

2.3. The Laplace Method

The Laplace method uses Laplace transforms to solve Equation (1) [33]. In matrix form, this solution can be written as follows:

$$\begin{array}{l}\left(\begin{array}{c}{\overline{\mathrm{T}}}_{\mathrm{int}}\\ {\overline{\mathrm{q}}}_{\mathrm{int}}\end{array}\right)={\displaystyle \prod _{\mathrm{k}=1}^{{\mathrm{N}}_{\mathrm{L}}}\left(\begin{array}{cc}\cosh \left(\sqrt{\mathrm{s}/{\alpha }_{\mathrm{k}}}\cdot {\mathrm{L}}_{\mathrm{k}}\right)& \frac{\sinh \left(\sqrt{\mathrm{s}/{\alpha }_{\mathrm{k}}}\cdot {\mathrm{L}}_{\mathrm{k}}\right)}{{\lambda }_{\mathrm{k}}\cdot \sqrt{\mathrm{s}/{\alpha }_{\mathrm{k}}}}\\ {\lambda }_{\mathrm{k}}\cdot \sqrt{\mathrm{s}/{\alpha }_{\mathrm{k}}}\cdot \sinh \left(\sqrt{\mathrm{s}/{\alpha }_{\mathrm{k}}}\cdot {\mathrm{L}}_{\mathrm{k}}\right)& \cosh \left(\sqrt{\mathrm{s}/{\alpha }_{\mathrm{k}}}\cdot {\mathrm{L}}_{\mathrm{k}}\right)\end{array}\right)\cdot \left(\begin{array}{c}{\overline{\mathrm{T}}}_{\mathrm{e}\mathrm{x}\mathrm{t}}\\ {\overline{\mathrm{q}}}_{\mathrm{e}\mathrm{x}\mathrm{t}}\end{array}\right)},\\ \hspace{1em}\hspace{1em}\hspace{1em}=\left(\begin{array}{cc}\mathrm{A}\mathrm{A}(\mathrm{s})& \mathrm{B}\mathrm{B}(\mathrm{s})\\ \mathrm{C}\mathrm{C}(\mathrm{s})& \mathrm{D}\mathrm{D}(\mathrm{s})\end{array}\right)\cdot \left(\begin{array}{c}{\overline{\mathrm{T}}}_{\mathrm{e}\mathrm{x}\mathrm{t}}\\ {\overline{\mathrm{q}}}_{\mathrm{e}\mathrm{x}\mathrm{t}}\end{array}\right)\end{array}$$

(17)

Here AA(s), BB(s), CC(s) and DD(s) are the elements of the wall transmission matrix H that relates to the surface temperatures and surface heat flux, with N_L being the number of layers, s being the frequency in the Laplace domain, and L_k and α being the thickness and the thermal diffusivity of each layer, respectively.

For the purpose of building performance simulation, and because DD(s) = AA(s), Equation (17) is transformed into Equation (18) [34].

$$\left(\begin{array}{c}{\overline{\mathrm{q}}}_{\mathrm{int}}\\ {\overline{\mathrm{q}}}_{\mathrm{e}\mathrm{x}\mathrm{t}}\end{array}\right)=\frac{1}{\mathrm{B}\mathrm{B}(\mathrm{s})}\left(\begin{array}{cc}\mathrm{A}\mathrm{A}(\mathrm{s})& -1\\ 1& -\mathrm{A}\mathrm{A}(\mathrm{s})\end{array}\right)\cdot \left(\begin{array}{c}{\overline{\mathrm{T}}}_{\mathrm{int}}\\ {\overline{\mathrm{T}}}_{\mathrm{e}\mathrm{x}\mathrm{t}}\end{array}\right),$$

(18)

In their work, [23] introduced a procedure that uses Z-Transforms [35] to convert the solution from the frequency domain (s) back to the time domain (t). This procedure is based on the partial fraction expansion [36,37] and finding roots of the denominator BB(s) on some predefined interval. The root-finding algorithm is one of the main drawbacks of the Laplace method, as it can lead to an excessive computation time and potential calculation errors [38]. Although the Laplace method is an analytical method, except for the root-finding procedure, which is mostly conducted numerically, the State-Space method has been shown to be superior for shorter time steps and heavyweight building elements [39,40]. Some researchers have tried to overcome the problem of root-finding such as [41], which proved to be efficient. They discovered that the roots of the heat transfer function BB(s) are separated by the roots of the transfer function AA(s). This new discovery ensured that no roots of the heat transfer function BB(s) are overlooked when performing the root-finding procedure. Nevertheless, they did not overcome the problem of using numerical methods for finding the actual roots, which was performed by the classical bisection method. In this paper, the root-finding procedure is performed using Mathematica’s FindRoot function and Brent’s method [42], in contrast to TRNSYS, which uses the bisection method.

From the CTF coefficients (b_j, c_j, d_j), the interior heat flux is calculated as follows:

$$q_{\mathrm{int}}(k\cdot \mathsf{\Delta }t)={\displaystyle \sum _{j=0}^{N_b}{b}_j\cdot T_{ext,k\cdot \mathsf{\Delta }t-j\cdot \mathsf{\Delta }t}}-{\displaystyle \sum _{j=0}^{N_c}{c}_j\cdot T_{\mathrm{int},k\cdot \mathsf{\Delta }t-j\cdot \mathsf{\Delta }t}}-{\displaystyle \sum _{j=1}^{N_d}{d}_j\cdot q_{\mathrm{int},k\cdot \mathsf{\Delta }t-j\cdot \mathsf{\Delta }t}},$$

(19)

where b_j is the interior CTF coefficient, c_j is the cross CTF coefficient and d_j is the flux CTF coefficient. N_b, N_c and N_d represent the number of b_j, c_j and d_j coefficients, respectively.

It is common practice to use Equation (20) to verify that the CTF coefficients obtained by the Laplace method give the correct steady-state heat transfer.

$$\frac{{\displaystyle \sum _{\mathrm{k}=1}^{{\mathrm{N}}_{\mathrm{a}}}{\mathrm{a}}_{\mathrm{k}}}}{{\displaystyle \sum _{\mathrm{k}=1}^{{\mathrm{N}}_{\mathrm{k}}}{\mathrm{d}}_{\mathrm{k}}}}=\frac{{\displaystyle \sum _{\mathrm{k}=1}^{{\mathrm{N}}_{\mathrm{b}}}{\mathrm{b}}_{\mathrm{k}}}}{{\displaystyle \sum _{\mathrm{k}=1}^{{\mathrm{N}}_{\mathrm{k}}}{\mathrm{d}}_{\mathrm{k}}}}=\frac{{\displaystyle \sum _{\mathrm{k}=1}^{{\mathrm{N}}_{\mathrm{b}}}{\mathrm{c}}_{\mathrm{k}}}}{{\displaystyle \sum _{\mathrm{k}=1}^{{\mathrm{N}}_{\mathrm{k}}}{\mathrm{d}}_{\mathrm{k}}}}\approx \mathrm{U},$$

(20)

2.4. EnergyPlus

Version of the software used in this research is version 9.6. The building elements were modelled in a single zone environment and only heat conduction through the wall was considered, i.e., convection and radiation on both internal and external surface were set to zero in EnergyPlus settings. That means that the temperature on of the internal and external surface was the same as the air temperature. EnergyPlus uses Steady-State method for the approximation of the solution for the heat equation as shown below. The same weather file used in Mathematica was also used in EnergyPlus.

The following equation shows the basic form of a CTF solution for the internal surface heat flux [43]:

$$q_{\mathrm{int}}(k\cdot \mathsf{\Delta }t)=-{\displaystyle \sum _{j=0}^{N_Z}{Z}_j\cdot T_{\mathrm{int},k\cdot \mathsf{\Delta }t-j\cdot \mathsf{\Delta }t}}+{\displaystyle \sum _{j=0}^{N_Y}{Y}_j\cdot T_{ext,k\cdot \mathsf{\Delta }t-j\cdot \mathsf{\Delta }t}}+{\displaystyle \sum _{j=1}^{N_{\mathsf{\Phi }}}{\mathsf{\Phi }}_j\cdot q_{\mathrm{int},k\cdot \mathsf{\Delta }t-j\cdot \mathsf{\Delta }t}},$$

(21)

where Z_j is the interior CTF coefficient, Y_j is the cross CTF coefficient and φ_j is the flux CTF coefficient. N_Z, N_Y and N_φ represent the number of Z_j, Y_j and φ_j coefficients, respectively.

2.5. Trnsys

Version of the software used in this research is 17. For consistency with both EnergyPlus and algorithms implemented in Mathematica for Steady-State and Laplace method algorithms, only conduction through the wall was modelled. Convection and radiation on the internal and external surface were not considered. The temperature of the internal surface was the same as the air temperature. The building element was also modelled in a single zone environment. The same weather file that was used in Mathematica was used in TRNSYS.

The instantaneous heat flux entering or leaving the zone for an exterior wall can be modelled according to Equation (19) as described in [44].

3. Results

In this section, the results of all the abovementioned methods are compared to check the efficiency of each method for different time steps and two types of building elements (i.e., heavyweight and lightweight) shown in

Table 2. Characteristics of the heavyweight building wall element.

Description	d	λ	cp	ρ
	(cm)	(W/(m K))	(J/(kg K))	(kg/m3)
Interior stucco	2	1.0	1000	1800
Heavyweight concrete	20	2.6	1000	2500
Mineral wool	16	0.035	1030	100
Exterior stucco	1	0.9	1000	1800
U-value: 0.21370 W/(m2 K)

and

Table 3. Characteristics of the lightweight building wall element.

Description	d	λ	cp	ρ
	(cm)	(W/(m K))	(J/(kg K))	(kg/m3)
Gypsum sheathing board	2.5	0.25	900	900
Mineral wool	5	0.032	1030	10
Oriented strand board (OSB)	1.5	0.13	1700	650
Mineral wool	16	0.032	1030	10
Gypsum sheathing board	1.3	0.25	900	900
Polystyrene	10	0.042	1260	30
U-value: 0.10857 W/(m2 K)

. The time steps used for the analysis were 1 h, 0.5 h, 0.25 h and 0.1 h. The thermal transmittance (U-value) of each building element was also calculated to determine the steady-state surface heat flux as:

$${\mathrm{q}}_{\mathrm{int}}=\mathrm{U}\cdot \left({\mathrm{T}}_{\mathrm{int}}-{\mathrm{T}}_{\mathrm{e}\mathrm{x}\mathrm{t}}\right),$$

(22)

The U-value in this research refers to surface-to-surface heat conduction that does not take into account the values of the surface heat transfer resistances due to the convection and radiation, i.e., R_si = 0 (m² K)/W and R_se = 0 (m² K)/W. U-values calculated this way cannot be used as the design values. If the design values are to be calculated from the material properties and thicknesses shown in Table 2 and Table 3, then the procedure and values of the surface heat transfer coefficients (h_si and h_se) given by ISO 6946 standard are to be used accordingly.

The boundary conditions used for the analysis were 20 °C for the interior surface temperature (T_int), and the values for the exterior surface temperatures (T_ext) were taken from [45] for the City of Zagreb (Figure 2). The data taken from [45] provide climate data with a 30 min time step, and to be able to use it for shorter time steps, the data are transformed into a continuous function using Mathematica’s algorithms. The calculation period was one month, with a warm-up period of six days to neutralize the effects of the initial conditions. Since only heat conduction through the elements was considered, the internal and external surface heat transfer coefficients were both taken as 0 W/(m² K).

To check the stability of the results, Equations (16) and (20) are used (

Table 4. Steady-state heat transfer for heavyweight wall element.

Time Step	State-Space	EnergyPlus	Laplace	TRNSYS
	${\displaystyle \sum \mathbf{X}}/{\displaystyle \sum \mathbf{\Phi }}$		${\displaystyle \sum \mathbf{a}}/{\displaystyle \sum \mathbf{d}}$
1 h	0.21379	0.21371	0.21370	0.21371
0.5 h	0.21378	0.21357	0.21370	0.21356
0.25 h	0.21379	0.16264	0.21370	0.18437
0.1 h	0.21372	5.19508	0.19196	843.73377
	${\displaystyle \sum \mathrm{Y}}/{\displaystyle \sum \mathsf{\Phi }}$		${\displaystyle \sum \mathrm{b}}/{\displaystyle \sum \mathrm{d}}$
1 h	0.21379	0.21370	0.21370	0.21370
0.5 h	0.21379	0.21370	0.21370	0.21370
0.25 h	0.21379	0.21376	0.21370	0.21385
0.1 h	0.21379	0.20323	0.21250	0.15343
	${\displaystyle \sum \mathrm{Z}}/{\displaystyle \sum \mathsf{\Phi }}$		${\displaystyle \sum \mathrm{c}}/{\displaystyle \sum \mathrm{d}}$
1 h	0.21379	0.21371	0.21370	0.21367
0.5 h	0.21378	0.21638	0.21370	0.21326
0.25 h	0.21379	0.25976	0.21370	0.32998
0.1 h	0.21373	14.49152	0.19549	364.89093

and

Table 5. Steady-state heat transfer for lightweight wall element.

Time Step	State-Space	EnergyPlus	Laplace	TRNSYS
	${\displaystyle \sum \mathbf{X}}/{\displaystyle \sum {\Phi }}$		${\displaystyle \sum \mathbf{a}}/{\displaystyle \sum \mathbf{d}}$
1 h	0.10803	0.10735	0.10857	0.10857
0.5 h	0.10803	0.10734	0.10857	0.10856
0.25 h	0.10803	0.10730	0.10857	0.10851
0.1 h	0.10803	0.10156	0.10857	0.11196
	${\displaystyle \sum \mathrm{Y}}/{\displaystyle \sum \mathsf{\Phi }}$		${\displaystyle \sum \mathrm{b}}/{\displaystyle \sum \mathrm{d}}$
1 h	0.10803	0.10735	0.10857	0.10857
0.5 h	0.10803	0.10736	0.10857	0.10857
0.25 h	0.10803	0.10736	0.10857	0.10857
0.1 h	0.10803	0.10723	0.10857	0.10796
	${\displaystyle \sum \mathrm{Z}}/{\displaystyle \sum \mathsf{\Phi }}$		${\displaystyle \sum \mathrm{c}}/{\displaystyle \sum \mathrm{d}}$
1 h	0.10803	0.10736	0.10857	0.10857
0.5 h	0.10803	0.10736	0.10857	0.10834
0.25 h	0.10803	0.11020	0.10857	0.10744
0.1 h	0.10803	0.14306	0.10857	0.33648

).

As a result of the analysis, the interior heat flux is calculated using Equations (6), (15), (19), (21) and (22) for time steps 1, 0.5, 0.25 and 0.1 h (Figure 3, Figure 4, Figure 5 and Figure 6). In addition, the total energy transferred through the interior surface of the element is calculated by integrating the heat flux function from t_min = 144 h to t_max = 864 h (Figure 6b) the results for EnergyPlus and TRNSYS are not shown because they are numerically unstable, as can be seen in Figure 7d.

Table 6. Algorithm computation time for each method.

Method	Computation Time (s)
Time Step	1 h	0.5 h	0.25 h	0.1 h	0.05 h
FDM	1.291	2.360	4.289	10.183	20.058
State-Space	1.296	2.010	3.454	7.905	15.393
Laplace	3.065	3.581	4.618	8.282	14.852

and Figure 8 show the computational time required to calculate the interior heat flux function in the period from t_min to t_max for the FDM, the State-Space method and the Laplace method as implemented in Mathematica.

4. Comparison and Discussion

Two different types of wall elements were modelled—heavyweight and lightweight. For each wall element, the internal heat flux was calculated using six different procedures to show the effect of the wall’s thermal mass on the solution for different time steps used for the dynamic thermal simulation for all six procedures.

The solution of the FDM was taken as the reference point for the comparison. The number of subdivisions for the FDM solution was 100 with the step size of 0.390 cm for the heavyweight and 0.363 cm for the lightweight wall element. The solution of the FDM is calculated as described in Section 2.1. The same subdivision used for the FDM was used for the State-Space method.

In Table 4 and Table 5, the stability of all procedures can be seen. For the heavyweight wall element, all methods show stable solutions for time steps of 1 h and 0.5 h (Figure 3a and Figure 4a). For time steps of 0.25 h and 0.1 h, the State-Space method and the Laplace method show stable solutions, while EnergyPlus and TRNSYS show a high deviation from the solution obtained with the FDM (Figure 5a and Figure 6a). For a time step of 0.25 h, this difference is as high as 31% for EnergyPlus and 78% for TRNSYS. In the case of a time step of 0.1 h, both EnergyPlus and TRNSYS give unreliable results. Although the algorithms of the State-Space method and the Laplace method implemented in the Mathematica for the purpose of this research are also implemented in EnergyPlus and TRNSYS, the heat flux calculated using EnergyPlus and TRNSYS is unstable for smaller time steps (Figure 5a and Figure 6a). This is due to the simplified procedures used for the calculation of matrix exponential in the case of the State-Space method implemented in EnergyPlus, and the root-finding algorithm for the Laplace method implemented in TRNSYS. In contrast, Mathematica uses a robust built-in algorithm for finding the matrix exponential and finding the roots of the denominator BB(s). In the case of the State-Space method, the number of subdivisions also affects the solution; however, not in terms of the stability of the CTF coefficients, but rather in terms of the overall accuracy of the calculated heat flux. As the number of subdivisions increases, the CTF coefficients for the State-Space method converge to a finite value.

For the lightweight wall element, all procedures show stable solutions for time steps of 1 h, 0.5 h and 0.25 h (Figure 3b, Figure 4b and Figure 5b). For time steps of 0.1 h, the State-Space method and the Laplace method show stable solutions, while EnergyPlus and TRNSYS show a high deviation from the solution obtained with the FDM (Figure 5b and Figure 6b). For a time step of 0.1 h, the difference is 46% for EnergyPlus and 303% for TRNSYS.

The difference in total energy transferred through the interior surface can be seen in Figure 7. For the lightweight wall element, the State-Space method and the Laplace method behave similarly for all time steps. There is no difference between the total energy transferred through the inner surface between the FDM and the State-Space method. For the Laplace method, this difference is negligible (0.39%). For the heavyweight wall element, the State-Space method is also stable for all time steps, and the difference between the total energy transmitted through the inner surface calculated using the FDM and the State-Space method is 0–0.04%. The Laplace method gives stable results for the time steps of 1 h, 0.5 h and 0.25 h, and the difference between the solution of the FDM and the Laplace method is 0.02%. For a time step of 0.1 h, the difference is 12.03%.

For the lightweight and heavyweight wall elements, EnergyPlus calculation produces more stable results than TRNSYS (Figure 7). For the lightweight wall element, the difference between the total energy transferred through the inner surface between the FDM and EnergyPlus is 0.65–3.11% for time steps of 1 h, 0.5 h and 0.25 h, and for the time step of 1 h, that difference is 45.64%. On the other hand, the difference between FDM and TRNSYS is 0.04–0.31% for the 1 h and 0.5 h time steps, and 31.01% for the 0.25 h time step. For the 0.1 h time step, the results are unreliable as can be seen in Figure 7d.

If the steady-state solution is considered, it can be seen that although the amplitude of the steady-state heat flux function is greater than that of the FDM, the difference between the total energy transferred through the inner surface between the FDM and the steady-state method is only 2.23% for the lightweight wall element and 2.19% for the heavyweight wall element (Figure 3, Figure 4, Figure 5, Figure 6 and Figure 7). However, when the steady-state method is used to calculate the heat flux, the thermal mass of the element is not taken into account, which could lead to potential problems when designing mechanical systems and thermal comfort for the building’s occupants. This is especially important since one of the future directions of building energy management lies in the implementation of model predictive control (MPC).

The computation times for the State-Space method, Laplace method and FDM are shown in Figure 8 and Table 6. For larger time steps, the State-Space method and the FDM perform similarly, while the Laplace method is much slower due to the interval for finding the roots of the denominator BB(s), which is 10⁻⁸ for both the lightweight and heavyweight wall elements. For smaller time steps, the computation time for the State-Space method and the Laplace method is similar. In contrast, the computation time of the FDM becomes higher due to the number of linear systems that must be solved for each time step.

The State-Space method and the Laplace method implemented in the conventional software for dynamic energy simulations software could be improved by implementing more advanced algorithms for the calculation of the matrix exponential for the State-Space method and the root-finding procedure for the Laplace method, as these are the main drawbacks of the CTF methods. As is shown in the results, by using more advanced algorithms, both the State-Space method and the Laplace method were stable for smaller time steps in comparison to the results obtained from EnergyPlus and TRNSYS. This improvement for the calculation of the matrix exponential and the root-finding procedure comes at the cost of the computational time. However, if the goals set by the EU for the NZEBs and ZEBs are to be achieved, it is essential to improve the accuracy of the dynamic calculation for smaller time steps.

5. Conclusions

Dynamic methods are based on the non-stationary solution of the heat equation. This solution can be found numerically or analytically. Because of the long analysis period and the small time steps used for the analysis, numerical methods have proven to be insufficient for building energy simulations. Analytical methods such as conduction transfer functions (CTFs) have proven to be much better for building energy simulation problems and are implemented in modern dynamic energy simulation software such as TRNSYS and EnergyPlus. Among the most commonly used CTF methods are the Laplace method and the State-Space method. This research compares these two methods for the calculation of heat losses in heavyweight wall elements by implementing the algorithms from the research paper of Stephenson and Mitalas for the Laplace method and Seem for the State-Space method. These results are then compared with the results from TRNSYS and EnergyPlus.

TRNSYS uses the Laplace method to calculate the CTF coefficients and the bisection method to calculate the roots of the function BB(s), while EnergyPlus uses the State-Space method to calculate the CTF coefficients and the simplified methods such as the algorithm by Moler and Van Loan to calculate the matrix exponential. The goal of this paper is to show how well the Laplace method and the State-Space method perform for small time steps and heavyweight building elements and to show what opportunities exist to improve these methods by implementing more efficient algorithms for root-finding and the matrix exponential procedure. The finite difference method (FDM) is used as a reference point for the boundary heat flux density.

For the time steps of 1 h and 0.5 h, all methods (State-Space, Laplace, EnergyPlus and TRNSYS) were found to be adequate for both the heavyweight and lightweight wall elements. For the heavyweight wall element and a time step of 0.25 h, the difference between the total energy transferred through the inner surface was about 31% for EnergyPlus and 78% for TRNSYS compared to the solution obtained with the FDM. For the lightweight wall element, the results were stable for the time step of 0.25 h, but for the time step of 0.1 h, the difference between EnergyPlus and FDM was 45.64% while the difference between TRNSYS and FDM was 303%.

The State-Space method and the Laplace method implemented in Mathematica were stable for all time steps for both the heavyweight and lightweight wall elements due to the Mathematica’s powerful numerical and matrix manipulation algorithms. The maximum difference between the FDM and Mathematica was for the Laplace method and a time step of 0.1 h, where the difference was 12.03%. Therefore, one should choose the time step size accordingly, and avoid the smaller time steps if possible.

For the calculation of the total energy consumption, the steady-state method proved to be satisfactory. The difference between the FDM and the steady-state method was around 2.23% for the lightweight wall element and 2.19% for the heavyweight wall element.

The conclusions of this paper could be of interest for the implementation in NZEBs and ZEBs where thermal mass influences are relatively higher as well as for better power demand management strategies, which require accurate transient simulation of heavyweight and highly insulated building elements with short time steps (lower than or equal to 0.1 h, i.e., 6 min). On the other hand, the limitation of the research is that it involves a complex technique and requires high computational costs. This paper shows that the State-Space method and the Laplace method implemented in EnergyPlus and TRNSYS could be improved by implementing more advanced algorithms for the calculation of the matrix exponential for the State-Space method and the root-finding procedure for the Laplace method since these are the main drawbacks of the CTF methods. Further research will address the integration of more efficient algorithms with equal or smaller deviations from theoretical values.