Thermal - electrical analogy in dynamic simulations of buildings: comparison of four numerical solution methods THERMAL - ELECTRICAL ANALOGY IN DYNAMIC SIMULATIONS OF BUILDINGS: COMPARISON OF FOUR NUMERICAL SOLUTION METHODS

: The lumped capacitance method is widely used in dynamic modelling of buildings. Models differ in complexity, solution methods and ability to simulate transient behaviour of described objects. The paper presents a mathematical description of a simple 1R1C thermal network model of a building zone. Four numerical methods were applied to solve differential equation describing its dynamics. For validation purposes two test cases (600 and 900) from the BESTEST procedure were used. In both cases detailed results were given. Better ability of the simulation model to reproduce transient behaviour of the modelled buildings was noticed in case of the lightweight object (case 600). Annual heating and cooling demand was within the reference range for heavyweight one (case 900). The kind of the computation method had no significant effect on simulation results.


INTRODUCTION
Determination of the amount of energy for space heating and cooling is crucial in auditing, certification and thermal modernisation of buildings [1]. For this purpose steady-state, quasi steady-state and dynamic methods [2] can be used. The latter ones enable transient analyses in hourly (or less) time step. They are also commonly applied in numerous simulation tools because of the rapid development of the computer technology. Dynamic methods include weighting factor and heat balance methods [3][4][5]. The second group contains response factor method, conduction transfer functions, finite difference method and lumped capacitance method [6][7].
Because of its simplicity and low computational demand the lumped capacitance method was commonly applied in researches [8][9][10][11][12][13] and standards, as VDI 6007 [14][15] or PN-EN ISO 13790 [16]. It is based on the thermal-electrical analogy [17][18]. In this approach the voltage source corresponds to the temperature, current -heat flux, resistor -thermal resistance and capacitor -thermal capacity. From this a dynamic model of a building in the form of an electrical scheme can be built. From the circuit theory there can be written nodal equations for each node of that network [19]. These are differential equations for nodes with connected capacitor (dynamic element) or algebraic equations for non-capacitive nodes.
As analytical solutions are usually difficult to obtain, numerical methods are applied in computer simulation tools. It should be done carefully, because there is a risk of instabilities between subsequent calculation time steps because of the non-continuous nature of boundary conditions given by weather, occupants or HVAC systems. Also significant changes in input values (outdoor temperature, ventilation rate, etc.) may occur between subsequent time steps [20][21][22][23]. These factors influence calculated instantaneous indoor temperature and thermal power required to maintain set point temperature. Thus, accurate modelling of transient states is important from the point of view of design procedures. This paper presents a simple RC lumped capacitance model of a building zone, under varying external and internal conditions, solved by four numerical methods. They were applied in a typical spreadsheet, not requiring specialized software.
Selected waveforms of the indoor temperature were shown and peak heating and cooling power and annual energy for space heating and cooling were calculated. Obtained results were discussed and finally concluding remarks were given.

Introduction
As aforementioned, a number of different schemes of RC networks have been presented in literature. To avoid complicated mathematical description a simple model with a single capacitor seemed to be a reasonable solution. It has been found that a number of researchers has presented such models in the last decades obtaining satisfactory accuracy [24][25][26][27][28][29][30][31][32][33][34][35][36][37][38]. Therefore, a single node 1R1C model (one capacitor and one resistor) thermal network model of a building zone ( Fig. 1) has been chosen. Obviously, this model makes it impossible to accurately account for a number of phenomena, but in spite of that it was applied in several studies [39][40][41][42][43] which confirmed its practical usability in different cases. The total heat transfer between the external and internal environment includes thermal transmittance through external partitions (walls, windows, roof and by floor to the ground) and by ventilation: The heat flux ϕg includes solar gains, internal gains and thermal power (heating or cooling) delivered directly to the internal air: The single capacitor includes the thermal capacity of the considered building.

Mathematical description
Assuming directions of heat fluxes presented in Fig. 1, from the Kirchhoff's Current Law (KCL) for the node Ti, it can be written that: where: and: dt dT Inserting Eqs. (4) and (5) into (3) yields: and finally: ( ) Eq. (7) is a first order ordinary differential equation in the general form of: with the initial condition: Assuming ϕg = 0 and Ti(0) = 0 Eq. (7) can be easily solved analytically [19]: The expression: is called the thermal time constant and indicates a time required for a given body (wall, floor, whole building, etc.) to respond to a change in its ambient temperature.

Introduction
Dynamic simulations of buildings are typically performed in hourly, or sub-hourly, time steps for the whole year. Boundary conditions, given by weather data, occupancy and control schedules are given in the same time resolution. Within that period, at consecutive time steps, all required boundary conditions are known and then differential equations are solved. Then output results are used to obtain solution at the next time step and so on until the end time is reached.
In order to clearly determine the current moment at which calculations are carried out, the convention for counting the simulation time should be established.

ISO 8601 standard
The ISO 8601 standard [44] recommends the use of time scales applying the 24-hour time keeping system for the identification of time points within a calendar day.
Hour is represented by two digits from 00 to 24. But the value of 24 is allowed only to indicate the end of a calendar day within a time interval. Minute is represented by two digits from 00 to 59. Second is represented by two digits from 00 to 60, but the value of 60 is allowed only to indicate a positive leap second or an instant within that second.
Beginning and end of day is represented by 00:00:00 and 24:00:00, respectively. But the end of one calendar day 24:00 coincides with 00:00 at the start of the next calendar day. Thus, the time 24:00 on 2 January is the same as 00:00 on 3 January. Hence, this method of time representation should be used in computer simulation programs with caution.

Simulation tools
Application of abovementioned recommendations to computer codes requires caution so as not to cause errors during transitions between successive hours.
Usually the default time step in simulation tools to calculate buildings thermal loads is one hour regardless of the fact that sub-hourly values are also possible. It is so, because weather data files are in hourly formats [45][46][47].
For example, in DOE-2, EnergyPlus and Helios programs a day is from hour 1 to hour 24. Hour 1 means the time interval from 00:00:01 to 1:00:00 (see Fig. 2). In DeST a day is from hour 0 to hour 23. Hour 0 represents the time interval from 00:00:00 to 00:59:59 [48][49]. The first convention was adopted in the paper. Omitting seconds for convenience it means that hour 1 covers the time interval from 00:01 to 01:00.

Introduction
To solve Eq. (7) finite difference methods can be applied. They use and approximation of the derivative from that equation by the arithmetic difference: forward, backward or central [50][51].

Forward Euler
Forward Euler method (explicit) uses forward difference approximation at the n-th time step: Typically the length of the time step is assumed to be constant for the whole simulation period: and this condition is followed in the paper. Hence, the general scheme of this method can be written as: Rearranging Eq. 14 we obtain: Looking at Eq. (15) it can be easily noticed that in an explicit formula the right side of that equation has only known variables. Taking the internal temperature Ti as the unknown variable y and applying Eq. (15) into Eq. (6) we get for the n-th calculation time step: The thermal capacity was assumed to be constant and therefore it wasn't assigned with the subscript "n".
From Eq. (16) the unknown internal temperature at the n+1 time step can be derived on the base of the value from the previous moment: . (17) This method is very easy to apply in computing procedures. Its disadvantage is the instability, which can be eliminated by the suitable choice of time step.
To obtain stability condition Eq. (8) should be rewritten in the form: Solution of Eqn. (18) is stable if: Comparing Eq. (19) with Eq. (7) it can be easily noticed that: and:

Backward Euler
Backward Euler method (implicit) uses backward difference approximation. It is similar to Eq. (12) Consequently, following the previous section, the general solution scheme for this method has the form: The unknown value yn+1 appears on both sides of the Eq. (23). Applying Eq. (23) to Eq. (7) we obtain: .  (24) This method is computationally more expensive per time step than the forward method, but it is numerically more stable.
In this case stability condition for Eq. (18) has the form:

Trapezoid method
This method is similar to the previous one, but more accurate and fast. For λ<0 it is absolutely stable. In general form applied to multiple dimensional partial differential equations it is known as the Crank-Nicholson method. It uses the approximation of the derivative at the midpoint by the expression: The general solution scheme for this method has the following form: The method is implicit because the unknown yn+1 appears on both sides of the Eq. (27).

Heun's method
This method is also known as improved Euler method. Instead of the derivative value calculated at the beginning of the interval, as it was in the Euler method, the derivative is also calculated at the end of the interval. The first result is called a predictor: Then: Heun's method uses forward Euler's method as predictor and trapezoid method as corrector. Thanks to the introduced modification this method is more accurate than the classic Euler method.
Applying Eq. (7) and (8)  Having given k1 and k2 variables they are inserted into Eq. (31) and the resulting temperature is obtained.
Stability condition for this method is described by the equation:

BESTEST procedure
To show an application of the aforementioned methods in thermal computations an exemplary test building from the BESTEST procedure [57] was used.
It was developed by the National Renewable Energy Laboratory in collaboration with the International Energy Agency.
It is a validation method based on comparative testing of building simulation programs including several test cases evaluating the influence of different physical processes on the simulation results. The basis for comparison is a range of results from a number of reference programs, such as EnergyPlus, BLAST, DOE2, COMFIE, ESP-r and others. A number of test cases of BESTEST have been incorporated into ANSI/ASHRAE Standard 140 [58].

Test cases
The BESTEST method covers a number of different tests to check buildings and their systems. To compare calculation algorithms shown in Section 4 there were chosen cases 600 and 900, i.e. base tests for low and heavy mass building, respectively.
In these tests the ANSI 140 standard uses an exemplary test building (Fig. 3) [59][60][61][62][63]. It also describes the thermal and physical properties of building elements, windows properties, schedules of internal gains, ventilation and thermostat control and heat transfer coefficients for all tests. The model of the outdoor climate is given in the form of 8760 hourly values of climatic data for Denver in the USA (39.8° N, 104.9° W, 1609 m above sea level). In the test 600 the building is of lightweight construction. Constant infiltration of 0.5 air change per hour is assumed. Internal load: 200 W continuous, 60% radiative, 40% convective, 100% sensible. It uses mechanical heating and cooling air system, 100% convective, 100% efficient, with no duct losses and no capacity limitation and no latent heat extraction. Temperature is controlled by non-proportional-type dual set-point thermostat with deadband. Internal setpoint temperatures: Ti,H,set = 20°C and Ti,C,set = 27°C. The test 900 uses the same building model as previously, except heavier materials used for construction of the wall and floor. Everything else with the building remained the same.
From the abovementioned data there were calculated values (Tab. 1) of elements of the RC model of the building (Fig. 1) The procedure to calculate hourly heating/cooling power ϕHC,nd,ac to maintain the required indoor air temperature (Fig. 4) was taken from [16]. Firstly, the calculation method (see section 4) to obtain internal temperature is chosen. Next, it is checked if heating or cooling is required (case 3 from For this value Ti is computed again, the resulting value is named as Ti,10 and then unrestricted heating or cooling power is calculated from the equation: ϕHC,nd,un = ϕHC,nd,10 (Ti,set -Ti,0)/(Ti,10 -Ti,0). (36) Now it is checked whether available cooling or heating power is sufficient (case 2 or case 4). If ϕC,max < ϕHC,nd,un < ϕH,max, then ϕHC = ϕHC,nd,un, Ti = Ti,set and the calculation is completed. If this condition is not satisfied, then (case 1 or case 5) the internal temperature is calculated for ϕHC,nd,un = ϕH,max, when ϕHC,nd,un > 0 or ϕHC,nd,un = ϕC,max, when = ϕHC,nd,un < 0.
In this case the set point temperature is not reached and calculations for the current time-step end.

Introduction
Calculations were performed in two ways. In the first one the analytical solution given by Eq. (10) was considered.
In the second part calculation models, described in the section 4, were applied in a spreadsheet along with the internal air temperature control strategy. Solar heat gains were determined according to [16]. Hourly solar irradiance incident on vertical surfaces was computed applying the Perez model [65].

Analytical solution
To compare accuracy of the presented methods with analytical solution given by Eq. (10) the test building from Fig. 3 was used. Assuming values of Htot and C from Table 1 and external air temperature step change of ∆Te = 50 K variation of internal air temperature was calculated using presented numerical methods. Results for test case 600 and 900 are shown in Fig. 5 and Fig. 6, respectively. In both cases Ti waveforms are the same. They differ only in time scale resulting from the thermal time constants of both buildings.
All methods provided good fitting of the calculated values with the analytical solution. But it was rather simple example. In case of varying conditions the situation is more complicated. Table 2 shows minimum and maximum reference values of annual heating and cooling energy and peak power from BESTEST and calculation results for all methods in the case 600.

BESTEST tests
Heating demand and peak cooling power were underestimated in all methods and amounted from 3751 kWh (3-rd method) to 4028 kWh (1) and from 4241 W (2) to 4710 W (1), respectively. So, their discrepancy was relatively small comparing to the computed average values.
Annual cooling energy was within the required range only in forward Euler method. Peak heating power was calculated properly by all methods. Based on the results presented in Table 2 and Table 3 the usefulness of the analysed calculation models for energy assessment of buildings can be determined. But to evaluate their ability to simulate transient states a deeper analysis is needed. It can be performed on the base of the more detailed charts shown in Fig. 7 and Fig. 8.
Regardless of the type of a building all methods reproduced the hourly heating and cooling load in the same way. The shape and variability of presented waveforms is the same. Maximum and minimum hourly values of thermal power differed from 116 W (7-th hour) to 1466 W (11) and from 45 W (17) to 223 W (12) in the case 600 and 900, respectively. Lower discrepancies between individual methods occurred for the heavy building case. Smaller differences, in relation to the reference values, were noticed for a light building.
Better compliance of computed hourly power values with the reference range was visible during low variability of the external conditions (hours 1-7 and 20-24) not influencing thermal balance of the building.
To illustrate this problem global solar irradiance (Isol) and external air temperature (Te) during 4-th January were shown in Fig. 9. Significant and rapid change in solar gains and external air temperature caused a decrease in the demand for heating power within a short period of time. Smaller discrepancies were noticed in case of the thermally light object, which react faster to changes in external and internal conditions. Differences between presented methods and reference waveforms are caused by the simplicity of the RC model. In this model solar irradiance incident on the external walls, calculated using the method from [16], is included on the internal side of a zone.

CONCLUSIONS
The lumped capacitance 1R1C thermal network model of a building zone was presented and numerically solved by four methods. In spite of its simplicity the model showed abilities to simulate thermal dynamics of simple buildings.
Presented simulations confirmed that numerical methods selected to solve presented problem affected final results depending on the type of the building. But literature review has revealed that not only calculation method, but the layout and physical background of the thermal network is very important, as well. It was confirmed by numerous experiences associated with the use of more complex models with one capacitor [66][67][68][69][70] indicating that its modification may result in better calculation accuracy while maintaining its simplicity and low computational demand.