Thermal—Airflow Coupling in Hourly Energy Simulation of a Building with Natural Stack Ventilation

Natural ventilation dominates in Polish residential buildings. It is a simple and low-cost system but its performance is affected by varying environmental conditions. Hence, setting up constant ventilation airflow results in errors when calculating heating and cooling energy. In this paper, an attempt to integrate the buoyancy effect in natural ventilation of a residential building at hourly resolution with the hourly simulation method of EN ISO 13790 to obtain energy use for space heating and cooling is presented. The ping-pong coupling algorithm was proposed and applied. Hourly variation of ventilation airflow rate was from −26.8 m3/h (flow from outdoor to the interior of the building) to 87.2 m3/h with 55 m3/h on average. The lack of a cooling system resulted in overheating during summer and indicated the necessity of its application or use of other techniques to reduce solar gains. Application of the cooling system resulted in an hourly ventilation rate from−38.0 m3/h to 87.2 m3/h. Detailed simulation in EnergyPlus and statistical analysis proved the applicability of the proposed method in stack-induced ventilation assessment. The coefficient of determination R2 = 0.936, mean squared error MAE = 5.72 m3/h and root mean square error RMSE = 7.86 m3/h.


Introduction
Regardless of the building's energy standard, an important element in its thermal balance is heat for heating (in the cool season) or cooling (in the warm season) of the ventilation air [1][2][3][4]. It is so because of hygienic reasons: continuous operation of ventilation, providing fresh outdoor air and removing polluted indoor air, is necessary during the presence of people.
Due to significant energy consumption by buildings and rising energy costs, numerous efficiency-related initiatives have been launched in many countries recently [5][6][7]. Additionally, new building codes directed to reduce transmission losses through external partitions have been introduced [8][9][10]. Consequently, the importance of proper assessment of ventilation heat loss has increased [11,12].
The ventilation systems are divided into [13] natural ventilation (NV), mechanical ventilation (MV) and hybrid ventilation (HV). The first solution, because of low investment costs and simplicity [14], is preferred in residential buildings.
According to the results of the National Census of Population and Housing from 2011 [15], residential buildings form the most significant part of the Polish buildings stock. Of them [16] only 16% are the objects built after the year 2002. In addition, only 5.5% of multifamily buildings retrofitted between 2010 and 2016 had the ventilation system modernised [17]. It follows from this that, among various ventilation systems, natural ventilation plays the most important role in the Polish residential sector. Therefore, the problem of assessment of ventilation heat loss in total energy consumed for heating and cooling of buildings is of great practical importance.
In general, the calculation, measurement and hybrid methods are used in the energy evaluation of buildings [18][19][20][21]. The first type is the most commonly used because of its advantages [22]. Additionally, if properly calibrated based on in-situ measurements, simulation methods can produce very reliable results [23].
As ventilation airflow in a building depends on ambient and geometric conditions, in the first group of studies the computational fluid dynamics (CFD) simulation tools were involved to assess indoor thermal comfort, pollution transport, air distribution or input boundary conditions for energy simulations in naturally ventilated objects.
Malkawi et al. [24] presented an assessment of cooling energy consumption in a residential building with mixed ventilation. The FloVENT and FlowDesigner programs were used in CFD simulations of natural ventilation and Airflow Network, the built-in module of EnergyPlus, for mixed-mode ventilation. The authors concluded that Airflow Network could not deal with complex airflow issues which resulted in differences in flow rate predictions in relation to CFD tools. On the other hand, about 30-40 CFD simulations for a one-year study on mixed-mode natural ventilation are required.
A similar combination of tools was used in [25] to obtain boundary conditions on the building's facade pressure (CFD), predict the natural ventilation's potential and assess the effect of the selected windows opening positions (Rhinoceros 3D) on the energy demand of the three-story building (EnergyPlus).
In [26] authors linked 3D parametric modeling and CFD to assess wind-driven natural ventilation in a residential test building. The Rhinoceros tool was used for the geometric model and OpenFOAM for CFD simulations.
CFD application in naturally ventilated buildings was also presented for offices [27] and industrial facilities [28]. These programs can easily cope with the thermal-airflow coupled computations. Such coupling also occurs in the analysed case as buoyancyinduced ventilation airflow depends on the difference between ambient air and indoor air temperatures. At the same time, the indoor temperature depends on indoor thermal conditions influenced by the amount of ventilation air entering a building's zone. This way an airflow-thermal coupling is established. The main four types of algorithms used to solve such couplings are distinguished in literature: the sequential method, ping-pong method, onion method and fully integrated method [29][30][31][32][33].
However, as CFD modelling is rather laborious, time-consuming, and requires the use of sophisticated commercial tools, it is preferred in the design studies rather than in energy assessment of buildings. For this purpose, building simulation tools, such as TRNSYS [34] or EnergyPlus [35][36][37][38] but also Autodesk Green Building Studio [39], IDA ICE [23], Matlab, DesignBuilder, Revit and Python [40], involving the aforementioned coupling methods, are preferred.
On the other hand, when the problem of fast energy calculations for certification, auditing or design purposes arises, then the related standards shall be used [41][42][43]. In such a case EN 13790 [44] or its recently introduced successor EN ISO 52016 [45] are recommended.
In [46] authors studied the use of the natural ventilation in a residential house in terms of a comfortable indoor environment and energy performance in the most efficient way in the early building design stage. CFD simulation was used to choose the most suitable architectural solution in terms of natural ventilation performance and indoor thermal comfort at the design stage of a residential building. Then, energy effects were evaluated using the monthly quasi-steady-state calculation method of EN ISO 13790.
In the study [47] on the energy efficiency of a residential multifamily building in Poland, annual heating energy was calculated for 180 building variants, including three types of ventilation (exhaust, exhaust with night reduction and mechanical ventilation with heat recovery). A monthly method from EN ISO 52016-1, similar to that of EN ISO 13790, was used.
This assumption is satisfactory in practice in pre-design studies in energy auditing or in the fast evaluation of retrofit strategies where standard calculations are often required. In step-ahead studies, regarding a more detailed assessment of thermal dynamics, the inclusion of time variation of a ventilation flow is also needed. It depends on environmental conditions, especially air temperature and wind velocity [57,58]. For these reasons, monthly methods are not suitable here and dynamic methods have been used recently to model varying airflow in buildings with natural ventilation: COMIS [59], TRNSYS [60], Energy-Plus [61] and Matlab [62]. These tools are of high quality but they are too complicated or expensive for use in less demanding applications. For these reasons, intermediate solutions are sought that take account of current standards and are quick to implement.
Jędrzejuk and Rucińska [63] proposed an interesting simple and efficient solution that could be applied here. In their analysis of artificial cooling in a single family building, they used the 5R1C model of a building from EN ISO 13790.
This model has several advantages, such as easy modification for various applications [64][65][66] and low computational requirements [67][68][69]. It can be also easily applied in a spreadsheet [70,71] not requiring specialised commercial simulation tools. Despite its simplicity, it produces reliable results [72]. For these reasons, it was used in studies where varying ventilation flow was considered, such as broiler house [73], residential building [74], greenhouse [75] or building-façade with integrated solar thermal collectors [76].
Because of its flexibility, the 5R1C model was also adopted in simulations related to ventilation. Narowski et al. [77] separated infiltration and ventilation airflows using two separate thermal conductances. The resulting 6R1C model was solved in a spreadsheet but the calculation procedure was more complicated than in the case of the generic 5R1C model.
Michalak [78] proposed the 4R1C model with equivalent ventilation heat flux. It was calculated from the difference between outdoor and indoor temperatures. To overcome circular references during computations in the developed method he proposed the internal air temperature at the given time step to be taken from the previous one. The model was successfully validated against EN 15265, BESTEST and EnergyPlus simulation with varying, ventilation-induced, airflow.
This review shows the lack of applications of simple hourly dynamic models used in energy performance evaluation of buildings in simulations with natural stack ventilation systems. Mainly, sophisticated CFD or building performance simulation tools are used for these solutions. Therefore, the novelty of the paper comes directly from these considerations. Its aim is to describe an application of the 5R1C dynamic model in the simulation of heating and cooling demand in hourly resolution including the effect of buoyancy-driven natural ventilation. The usability of the proposed solution is confirmed by comparison with a more advanced tool.
The next section briefly presents a calculation of an airflow rate in buoyancy-driven ventilation, the 5R1C resistance-capacitance model of a building from EN ISO 13790 and a mathematical model to obtain stack-induced ventilation airflow. Then, the test building, weather data and simulation assumptions are given. After them, the results are presented, and concluding remarks are given.

Modelling of Natural Ventilation
Two types of natural ventilation are met in buildings. Wind-driven ventilation is caused by the different pressures created by the wind on a building's envelope and by openings on the perimeter forming a flow path for ventilation air. Buoyancy-driven ventilation is based on differences between the density of warmer indoor and colder outdoor air, resulting in a flow of indoor air from a building's interior through higher located openings. At the same time, cooler and denser outdoor air flows into a building through inlet openings [79][80][81][82]. The physical nature of natural ventilation makes it very difficult to model [36]. However, numerous experimental and theoretical studies conducted in recent decades have led to convenient mathematical relationships used to compute ventilation airflow in various conditions and technical solutions met in residential buildings.
As Poland lies in the heating-dominant climate, when the difference between indoor and outdoor temperature is larger than during the remaining periods of the year, the buoyancy effect dominates over the wind effect. Therefore, in Polish residential buildings, both single and multifamily, stack ventilation dominates [83]. Several experiments carried out in Poland [84,85] confirmed that during winds with low velocity the air change rate (ACH) resulting from the stack effect slightly exceeded the required standard minimum value. For these reasons and to simplify presented considerations, only the buoyancy effect was taken into account.
In the case of buoyancy-driven ventilation in a building with two openings: top A top and bottom A bottom the ventilation airflow rate is given by the relationship [86,87]: where A eff is the effective opening area, given by: and C d is the discharge coefficient (a typical value is 0.6), g is gravitational acceleration, h is the height between openings, and T i and T e are internal and external temperatures, respectively.

The 5R1C Model
The EN ISO 13790 standard introduced the thermal network model of a building zone consisting of five thermal resistances and one capacitance ( Figure 1).
Two types of natural ventilation are met in buildings. Wind-driven ventilation is caused by the different pressures created by the wind on a building's envelope and by openings on the perimeter forming a flow path for ventilation air. Buoyancy-driven ventilation is based on differences between the density of warmer indoor and colder outdoor air, resulting in a flow of indoor air from a building's interior through higher located openings. At the same time, cooler and denser outdoor air flows into a building through inlet openings [79][80][81][82].
The physical nature of natural ventilation makes it very difficult to model [36]. However, numerous experimental and theoretical studies conducted in recent decades have led to convenient mathematical relationships used to compute ventilation airflow in various conditions and technical solutions met in residential buildings.
As Poland lies in the heating-dominant climate, when the difference between indoor and outdoor temperature is larger than during the remaining periods of the year, the buoyancy effect dominates over the wind effect. Therefore, in Polish residential buildings, both single and multifamily, stack ventilation dominates [83]. Several experiments carried out in Poland [84,85] confirmed that during winds with low velocity the air change rate (ACH) resulting from the stack effect slightly exceeded the required standard minimum value. For these reasons and to simplify presented considerations, only the buoyancy effect was taken into account.
In the case of buoyancy-driven ventilation in a building with two openings: top Atop and bottom Abottom the ventilation airflow rate is given by the relationship [86,87]: where Aeff is the effective opening area, given by: and Cd is the discharge coefficient (a typical value is 0.6), g is gravitational acceleration, h is the height between openings, and Ti and Te are internal and external temperatures, respectively.

The 5R1C Model
The EN ISO 13790 standard introduced the thermal network model of a building zone consisting of five thermal resistances and one capacitance ( Figure 1).  The external envelope of a building is represented by two kinds of partitions. Thermally "light" elements (doors, windows, curtain walls, glazed walls, etc.) are described by the H tr,w thermal transmission coefficient. Thermally "heavy" elements (walls, ceilings) are included in the H tr,op thermal transmission coefficient. It is then split into the external (H tr,em ) and the internal (H tr,ms ) parts: with: The single capacitor, C m , represents the thermal mass of the building. The external environment is represented by two variables. The first one, external air temperature (T e ) is connected with H tr,w and H tr,em coefficients. The supply air temperature (T sup ) is connected through the heat transfer by ventilation (H ve ) with the internal air temperature (T i ). In the case of natural ventilation: The indoor air node (T i ) and the central node (T s ) are linked through the coupling conductance H tr,is given by Equation: The heat fluxes from internal heat sources (ϕ int ) and solar irradiance (ϕ sol ) are split into three parts: ϕ ia , ϕ st , and ϕ m , respectively, connected to the indoor air, the central node, and the thermal mass temperature nodes. Heating or cooling power (ϕ HC ) is supplied to (heating) or extracted from (cooling) the indoor air node.

The Calculation Algorithm
In this paper, the ping-pong coupling method (called also loose coupling, quasidynamic or decoupled approach) was used. In this approach, the thermal and the airflow models are run in sequence. At a given time step the output from the airflow model is passed to the thermal model. Then, the thermal model returns its output to the input of the airflow model at a subsequent time step.
In the considered model the ventilation heat transfer coefficient is computed from the relationship: with the volumetric heat capacity of air assumed to be constant at ρ a c a = 1200 J/m 3 K according to EN ISO 13790. Inserting Equation (1) into Equation (7) we get: Taking into account the aforementioned considerations, we get the calculation algorithm of hourly heating and cooling energy demand presented in Figure 2. At the first time step (n = 1) there should be assumed indoor air temperature from the previous, initial step. As it is unknown, it can be assumed as the set internal air temperature for heating, Tint,H,set or for cooling, Tint,C,set if calculations start in the heating or cooling period, respectively. As simulations are usually performed in hourly time step At the first time step (n = 1) there should be assumed indoor air temperature from the previous, initial step. As it is unknown, it can be assumed as the set internal air temperature for heating, T int,H,set or for cooling, T int,C,set if calculations start in the heating or cooling period, respectively. As simulations are usually performed in hourly time step for the whole year, the step number, n, changes sequentially from 1 to 8760.
As the ventilation rate is computed by taking the indoor air temperature from the previous step it may be stated here that no convergence is reached at the end of each calculation step. It is not so because the required thermal power is calculated to obtain a set internal air temperature. Additionally, it should be remembered that indoor air temperature fluctuations are relatively small when compared to ambient air. It is so, especially during the heating period when the indoor-outdoor difference is the largest and the resulting stack effect is the most significant.
The block "Solve the circuit and compute T i " contains the procedure to obtain heating or cooling power (ϕ HC ) required to maintain a certain set-point indoor air temperature. This procedure starts from the calculation of internal (ϕ int ) and solar (ϕ sol ) heat gains. Then they are split into three parts, as follows: The solution is based on a Crank-Nicholson scheme. The temperatures are the average over one hour except for thermal mass temperature T m,t and T m,t−1 which are instantaneous values at time t and t − 1, respectively.
T m temperature at a given time step, t, is given by: where: H tr,2 = H tr,1 + H tr,w (15) and: At each time step, the average values of nodes temperatures are given by: and the internal air temperature: Control of T i can be set to maintain its value depending on the user's requirements. The procedure to calculate the required hourly thermal power, ϕ HC , is given in the standard and has been presented by several authors recently [66,73,75,88]. It can be easily implemented in a spreadsheet which makes it very flexible and available for a wide audience.

Modelling of Stack Effect with EnergyPlus
As no measurements were taken to provide a basis for comparison there was chosen a detailed simulation in EnergyPlus. Both infiltration and ventilation can be modelled but, depending on the user's choice, they can be disabled when needed. This way only infiltration or ventilation can be taken into account or even a situation when there is no air exchange inside a building may be also analysed.
This tool offers a wide range of possible options. However, for convenience, there are presented those that can be used for direct comparison. Obviously, for comparative purposes, among various available solutions and algorithms of interest are those expressed in the same form as Equation (1).
The first model is of general form and can be used to evaluate infiltration or ventilation flow rate. It is the "Design Flow Rate" method in which airflow rate is given by: If wind impact is not considered then C = 0 and D = 0. The same can be done with the A factor. However, the proper setting of B is rather problematic. In [89] there were only given general indications taken from other sources.
The first model to evaluate infiltration flow rate due to stack effect is the effective leakage area model given by the relationship: In the considered case C w = 0 should be assumed. For a one-story house, C s = 0.000145 (L/s) 2 /(cm 4 ·K) is suggested.
The second model is infiltration by flow coefficient expressed as: When ignoring the wind effect it is similar to Equation (20). In general, the pressure exponent is set between 0.6 and 0.7. Usually n = 0.67 is used. In this model, for a one-story house, C s = 0.069 (Pa/K) n or C s = 0.054 (Pa/K) n is recommended with and without flue, respectively.
When ventilation is modelled the "Wind and Stack with Open Area" option can be used. Setting "Opening effectiveness" parameter in the program to zero, wind-induced flow is omitted and we obtain total volumetric airflow rate due to stack effect given by: In practice, the most difficult is the estimation of ∆H NPL , the height from the midpoint of the lower opening to the neutral pressure level. The reader is referred to Chapter 16 of the ASHRAE Handbook of Fundamentals [90] for guidance.
From the above, it can be easily deduced that Equation (23) is the most similar to Equation (1). In such a case, F schedule = 1, A opening = A eff and ∆H NPL = h/2 should be set.

Model Evaluation
When no measurement data are available, the detailed simulation in EnergyPlus can be considered as a reference for the evaluation of the presented model. Then, its performance can be described by several statistical measures. Assumingx i is the predicted (by the model) value of a given quantity (x), x i is its actual (reference) value, x i is the average value of x and n is the total number of analysed samples, then it can be written that [91][92][93][94] mean absolute error is given by: This metric is used to evaluate the average deviation between the predicted value and the true value.
The root mean square error: measures how much, on average, the predicted value deviates from the actual (reference) value. Hence, a small RMSE is desired. It is closely connected with the mean square error: providing the average of the summation value of the square error of the prediction compared to actual (reference) values. The operation of squaring enhances the importance of larger errors. Hence, when comparing MAE and RMSE, the RMSE is more useful when large errors are especially undesirable. The last considered statistic is the coefficient of determination: measures the amount of variance of the predicted variables given by the model in the case of the most optimal prediction model R 2 = 1.

Case Building
To perform necessary simulations, a single-story residential building without a cellar was chosen, located in Nowy Sącz, the city located in southern Poland in the third climatic zone according to the PN-EN 12831 standard.
According to the typical meteorological year (TMY) for Nowy Sącz [95] mean annual temperature was 8.46 • C and varied from 0.51 • C in January to 19.65 • C in July. Hourly air temperature (Figure 3) was from −12.0 • C on 8 January at 6:00 to 32.2 • C on 5 July at 12:00. Wind velocity varied from 0 to 10 m/s on 24 February. For 7757 h (88.6% of the year) wind velocities were lower than 3 m/s which means that the impact of wind on natural ventilation is insignificant [85].
The building is placed on a rectangular plan of 10 m × 7 m and is inhabited by four persons. Longer walls (front and back) are oriented in the east-west axis. It has a usable area of 58. According to the typical meteorological year (TMY) for Nowy Sącz [95] mean annual temperature was 8.46°C and varied from 0.51 °C in January to 19.65 °C in July. Hourly air temperature (Figure 3) was from −12.0°C on 8 January at 6:00 to 32.2 °C on 5 July at 12:00. Wind velocity varied from 0 to 10 m/s on 24 February. For 7757 h (88.6% of the year) wind velocities were lower than 3 m/s which means that the impact of wind on natural ventilation is insignificant [85].        To calculate the values of elements of the resistance-capacitance circuit from Figure 1 the physical properties of the materials obtained from manufacturers were used. The design heat transmission coefficients were determined following PN-EN ISO 6946. Thermal capacities of the building's partitions were calculated using the detailed method from PN-EN ISO 13786 for a calculation period of 24 h. The final results are given in Table 1. Heat transfer by ventilation was computed using the procedure given in Section 2.3 assuming volumetric heat capacity of air of 1200 J/(m 3 K) recommended in EN ISO 13790. The model was then simulated using a spreadsheet.

Simulation Cases
For comparative purposes, two simulations were performed. In the first case, it was assumed that no cooling system is in the building. The model of buoyancy-induced stack ventilation, described in Section 2.1, was used. In the second case, the cooling system was added. The last one was a detailed simulation in EnergyPlus tool with no cooling, as in the first case.

Case 1
In the first case, annual energy for heating amounted to 5493 kWh. It varied from 6 kWh in June to 1045 kWh in January. Hourly heating power varied ( Figure 5) from 0 to 3.02 kW on 8 January at 6:00.

Simulation Cases
For comparative purposes, two simulations were performed. In the first case, it was assumed that no cooling system is in the building. The model of buoyancy-induced stack ventilation, described in Section 2.1, was used. In the second case, the cooling system was added. The last one was a detailed simulation in EnergyPlus tool with no cooling, as in the first case.

Case 1
In the first case, annual energy for heating amounted to 5493 kWh. It varied from 6 kWh in June to 1045 kWh in January. Hourly heating power varied ( Figure 5) from 0 to 3.02 kW on 8 January at 6:00. Ventilation heat flux ( Figure 6) was from −921.0 W on 8 January at 5:00 to 21.4 W on 6 September at 12:00. It means that, except for short periods, the heat was extracted from the interior of the building by ventilation air during the year. It has a negative impact on thermal balance of the building during the heating season, increasing heating demand. In hot periods this effect means that the excess heat is removed, thus enhancing the building's cooling, especially in case of large solar gains.  Ventilation heat flux ( Figure 6) was from −921.0 W on 8 January at 5:00 to 21.4 W on 6 September at 12:00. It means that, except for short periods, the heat was extracted from the interior of the building by ventilation air during the year. It has a negative impact on thermal balance of the building during the heating season, increasing heating demand. In hot periods this effect means that the excess heat is removed, thus enhancing the building's cooling, especially in case of large solar gains.

Simulation Cases
For comparative purposes, two simulations were performed. In the first case, it was assumed that no cooling system is in the building. The model of buoyancy-induced stack ventilation, described in Section 2.1, was used. In the second case, the cooling system was added. The last one was a detailed simulation in EnergyPlus tool with no cooling, as in the first case.

Case 1
In the first case, annual energy for heating amounted to 5493 kWh. It varied from 6 kWh in June to 1045 kWh in January. Hourly heating power varied ( Figure 5) from 0 to 3.02 kW on 8 January at 6:00. Ventilation heat flux ( Figure 6) was from −921.0 W on 8 January at 5:00 to 21.4 W on 6 September at 12:00. It means that, except for short periods, the heat was extracted from the interior of the building by ventilation air during the year. It has a negative impact on thermal balance of the building during the heating season, increasing heating demand. In hot periods this effect means that the excess heat is removed, thus enhancing the building's cooling, especially in case of large solar gains.   The negative value of ventilation airflow means that air is flowing from the building's interior to the outdoor environment. The opposite situation may take place during hot summer months, which can be clearly stated when analysing indoor air temperature ( Figure 8). Hence, either passive or active cooling is required. The negative value of ventilation airflow means that air is flowing from the building's interior to the outdoor environment. The opposite situation may take place during hot summer months, which can be clearly stated when analysing indoor air temperature ( Figure 8). Hence, either passive or active cooling is required. The negative value of ventilation airflow means that air is flowing from the building's interior to the outdoor environment. The opposite situation may take place during hot summer months, which can be clearly stated when analysing indoor air temperature ( Figure 8). Hence, either passive or active cooling is required.

Case 2
In this case, an additional cooling system was assumed. Annual energy for heating and cooling amounted to 5496 kWh and 143 kWh, respectively. It varied from 7 kWh in June to 1045 kWh in January and from 2 kWh in May to 84 kWh in July. Hourly peak heating and cooling power ( Figure 9) amounted to 3.03 kW on 8 January at 6:00 and 0.9 kW on 6 July at 14:00.

Case 2
In this case, an additional cooling system was assumed. Annual energy for heating and cooling amounted to 5496 kWh and 143 kWh, respectively. It varied from 7 kWh in June to 1045 kWh in January and from 2 kWh in May to 84 kWh in July. Hourly peak heating and cooling power ( Figure 9) amounted to 3.03 kW on 8 January at 6:00 and 0.9 kW on 6 July at 14:00. The negative value of ventilation airflow means that air is flowing from the building's interior to the outdoor environment. The opposite situation may take place during hot summer months, which can be clearly stated when analysing indoor air temperature ( Figure 8). Hence, either passive or active cooling is required.

Case 2
In this case, an additional cooling system was assumed. Annual energy for heating and cooling amounted to 5496 kWh and 143 kWh, respectively. It varied from 7 kWh in June to 1045 kWh in January and from 2 kWh in May to 84 kWh in July. Hourly peak heating and cooling power ( Figure 9) amounted to 3.03 kW on 8 January at 6:00 and 0.9 kW on 6 July at 14:00.  Although the extreme values are the same as in the previous case, the waveform shape during the hot period when cooling was active has changed. It was directly connected with a higher indoor-outdoor temperature difference in relation to the previous case.
Assuming the annual average value of the ventilation airflow rate of 53.2 m 3 /h annual energy for heating and cooling of 6326 kWh and 374 kWh was obtained, respec- Although the extreme values are the same as in the previous case, the waveform shape during the hot period when cooling was active has changed. It was directly connected with a higher indoor-outdoor temperature difference in relation to the previous case.
Assuming the annual average value of the ventilation airflow rate of 53.2 m 3 /h annual energy for heating and cooling of 6326 kWh and 374 kWh was obtained, respectively. It varied from 8 kWh in June to 1183 kWh in December and from 5 kWh in May to 221 kWh in July.
In this consideration, the assumed discharge coefficient of 0.6 was assumed. However, in [85] from measurements, the value of 0.83 was obtained. This parameter is the most troublesome for estimation. Hence, an uncertainty connected with its proper calculations affects the final results.

Simulation in EnergyPlus
As aforementioned, simulations in EnergyPlus reflected assumptions set in the first case. Analysis of the air infiltration rate was given in this section to give a better view of the possibilities of the analysed model.
The annual average value of the ventilation airflow rate was 56.9 m 3 /h when not considering flow direction and 55.5 m 3 /h if summing both positive and negative values. However, regardless of the direction of flow, contaminated indoor air is replaced by fresh outdoor air. Hence, for design purposes, the indoor-outdoor temperature difference in Equation (1) is computed as the absolute value.
In this study, there were separated situations with opposite flow directions ( Figure 11). Calculated coefficient of determination R 2 = 0.936 indicated a very strong correlation. Its value means that 96.7% of variations in the reference ventilation airflow from EnergyPlus can be explained by the proposed model.
Although the extreme values are the same as in the previous case, the waveform shape during the hot period when cooling was active has changed. It was directly connected with a higher indoor-outdoor temperature difference in relation to the previous case.
Assuming the annual average value of the ventilation airflow rate of 53.2 m 3 /h annual energy for heating and cooling of 6326 kWh and 374 kWh was obtained, respectively. It varied from 8 kWh in June to 1183 kWh in December and from 5 kWh in May to 221 kWh in July.
In this consideration, the assumed discharge coefficient of 0.6 was assumed. However, in [85] from measurements, the value of 0.83 was obtained. This parameter is the most troublesome for estimation. Hence, an uncertainty connected with its proper calculations affects the final results.

3.4.Simulation in EnergyPlus
As aforementioned, simulations in EnergyPlus reflected assumptions set in the first case. Analysis of the air infiltration rate was given in this section to give a better view of the possibilities of the analysed model.
The annual average value of the ventilation airflow rate was 56.9 m 3 /h when not considering flow direction and 55.5 m 3 /h if summing both positive and negative values. However, regardless of the direction of flow, contaminated indoor air is replaced by fresh outdoor air. Hence, for design purposes, the indoor-outdoor temperature difference in Equation (1) is computed as the absolute value.
In this study, there were separated situations with opposite flow directions ( Figure  11). Calculated coefficient of determination R 2 = 0.936 indicated a very strong correlation. Its value means that 96.7% of variations in the reference ventilation airflow from Ener-gyPlus can be explained by the proposed model. Figure 11. Hourly ventilation airflow from EnergyPlus and the presented model. Figure 11. Hourly ventilation airflow from EnergyPlus and the presented model.
The average deviation between the value from the analysed model and the reference value from EnergyPlus measured by MAE = 5.72 m 3 /h. The average squared difference between the reference and predicted values given by MSE = 61.75 m 6 /h 2 . Finally, the average distance between the predicted and reference values given by RMSE = 7.86 m 3 /h.
As MAE and RMSE, when compared to the annual average ventilation airflow, are quite similar, it can be deduced that there are no significant large errors in the predictions from the proposed model. The histogram of hourly differences between predictions from the 5R1C model and detailed simulation ( Figure 12)  between the reference and predicted values given by MSE = 61.75 m 6 /h 2 . Finally, the average distance between the predicted and reference values given by RMSE = 7.86 m 3 /h. As MAE and RMSE, when compared to the annual average ventilation airflow, are quite similar, it can be deduced that there are no significant large errors in the predictions from the proposed model. The histogram of hourly differences between predictions from the 5R1C model and detailed simulation ( Figure 12) proves this outcome.

Conclusions
In this study, the thermal network model of a building zone from the EN ISO 13790 standard was linked with the mathematical model of buoyancy-driven natural ventilation to calculate the hourly heating and cooling energy demand of a residential building.
To couple airflow and thermal models, the so-called ping-pong coupling algorithm was applied. This way, calculations were performed by avoiding circular references and using the generic method provided by EN ISO 13790 applied in a spreadsheet. The chosen thermal model of a building showed its applicability in combination with the proposed algorithm. The method is simple and does not require the use of sophisticated building simulation tools. Presented statistical analysis proved its quality sufficient for simple calculations in terms of simple energy analyses. It also enables predictions of the airflow direction and ventilation heat flux which is the first step to assessing the necessity to apply, for example, night cooling or other ventilation strategies to minimise overheating risk or preheating to reduce ventilation loss.
It should be emphasized that the simplifications associated with this type of description of natural ventilation will certainly deviate from the results obtained from CFD tools or measurements. The simplifications associated with this type of description of natural ventilation will certainly deviate from the results obtained from CFD tools or measurements. Nevertheless, simulation in an hourly step with the choice of Cd according to the results of accurate studies allows the buoyancy phenomenon to be taken into account more accurately than, for example, an averaged ventilation flow over the year, independent of indoor and outdoor conditions.
In the next step, other coupling algorithms could be studied and comparisons with the more detailed tools can be used for comparison. It also seems reasonable to include the impact of wind velocity and direction into the calculation procedure to more accurately include the impact of environmental conditions on the performance of natural ventilation.

Conclusions
In this study, the thermal network model of a building zone from the EN ISO 13790 standard was linked with the mathematical model of buoyancy-driven natural ventilation to calculate the hourly heating and cooling energy demand of a residential building.
To couple airflow and thermal models, the so-called ping-pong coupling algorithm was applied. This way, calculations were performed by avoiding circular references and using the generic method provided by EN ISO 13790 applied in a spreadsheet. The chosen thermal model of a building showed its applicability in combination with the proposed algorithm. The method is simple and does not require the use of sophisticated building simulation tools. Presented statistical analysis proved its quality sufficient for simple calculations in terms of simple energy analyses. It also enables predictions of the airflow direction and ventilation heat flux which is the first step to assessing the necessity to apply, for example, night cooling or other ventilation strategies to minimise overheating risk or preheating to reduce ventilation loss.
It should be emphasized that the simplifications associated with this type of description of natural ventilation will certainly deviate from the results obtained from CFD tools or measurements. The simplifications associated with this type of description of natural ventilation will certainly deviate from the results obtained from CFD tools or measurements. Nevertheless, simulation in an hourly step with the choice of C d according to the results of accurate studies allows the buoyancy phenomenon to be taken into account more accurately than, for example, an averaged ventilation flow over the year, independent of indoor and outdoor conditions.
In the next step, other coupling algorithms could be studied and comparisons with the more detailed tools can be used for comparison. It also seems reasonable to include the impact of wind velocity and direction into the calculation procedure to more accurately include the impact of environmental conditions on the performance of natural ventilation.

Conflicts of Interest:
The author declares no conflict of interest. h is heat transfer coefficient between the air node, T i , and the surface node, T s , with a fixed value of h is = 3.45 W/m 2 K; n pressure exponent,q design design ventilation airflow rate, m 3 /s s shelter factor,w local wind speed, m/s ρ a air density, kg/m 3 ϕ int heat flow rate due to internal heat sources, W ϕ sol heat flow rate due to solar heat sources, W ϕ ia heat flow rate to the internal air node, W ϕ st heat flow rate to the central node, W ϕ m heat flow rate to the mass node, W ϕ ve heat flow rate by ventilation, W ϕ HC heating or cooling power supplied to or extracted from the indoor air node, W ∆H NPL Height from the midpoint of lower opening to the neutral pressure level, m