Thermal comfort: standard compliant modelling of natural ventilation in adaptive time step solvers

Natural ventilation can help maintain thermal comfort during summer periods. In simulation models for thermal comfort analysis, occupant behaviour related to natural ventilation is often included through definition of ventilation rule sets. For example, windows are opened when room temperature exceeds a given limit while outside temperatures remain below a threshold. Such rules are included in building codes, such as DIN 4108-2 [1]. However, these rules are defined with respect to fixed step time integration methods, often using hourly step sizes. Straight-forward implementation of such rules in error controlled adaptive time step integration solvers (including drivers for Modelica models) fails. The article discusses strategies and model variants that work with such modern solvers and discusses new model parameters such as the delay time or dead-band size. The impact of these parameter choices is illustrated with sensitivity studies using the thermal room model THERAKLES [3-5].


Introduction
The importance of calculating summer heat protection has increased in recent years. This may be due to architecture with ever larger glass facades, smaller storage masses or higher internal loads. Due to the increasingly complex requirements, the old analytical methods, which have been mostly used until now, are reaching their limits. Numerical methods have been known for a long time, but have only recently been used more widely. This is mainly due to the fact that these procedures are now also taken into account in the standards. In Germany, for example, since 2013, the relevant standard DIN 4108-2 [1] contains a chapter with requirements and boundary conditions for the use of numerical calculation methods. However, no calculation method is prescribed or validation or certification of such a method is required. Only a reference is made to ISO 13791 [2].
In numerical calculation methods, two basic time integration types can be distinguished. There are methods with constant time steps and those with variable time step size adjustment. The former are easier to implement and therefore account for the largest share of the procedures used. Methods with dynamic time-step control have some advantages, although they are more complex and challenging to implement. Table 1 shows a summary of the major differences between constant and variable step methods.

Dynamic step size adjustment based on local truncation error analysis
Modern building simulation models such as in THERAKLES [3][4][5], IDA-ICE [6], NANDRAD [7] and most Modelica [8] engines use state-of-the-art adaptive step size time integration methods (DASSL [9], CVODE [10], etc.). There are different types of methods and mathematical formulations, yet with the same basic idea: • construct a solution at the next time level with some estimated time step size using history information only (a predictor equation, using for example the current solution and time derivatives of solution variables) • compute a solution at the next time level using the system equations (implicit solution, usually involving solving of a non-linear equation system; also called the corrector equation, which different method order than the predictor equation) • evaluate the difference between predictor and corrector equation and estimate the local truncation error (in case of several equations, a vector norm of the error estimates is obtained) • compare the error estimate with the acceptable limit; if limit is exceeded, reduce time step size accordingly and try again, otherwise proceed and potentially enlarge time step when error estimate is small enough  Step sizes used by an adaptive step size integrator. Note that integration steps need not be aligned with output time points. Instead, interpolation (linear/higher order) is used to obtain output values at designated time grid points (hourly values etc.) With dynamic time-step control, special challenges arise with step-shaped control strategies with strong dependence on a state variable in order to prevent unstable behaviour (oscillations, very small time steps). Possible solutions to this problem will be illustrated using the example of ventilation control according to the German standard DIN 4108-2.

Ventilation model DIN 4108-2
The German standard DIN 4108-2 [1] describes two methods for determining summer heat protection. In addition to a simple manual calculation method (solar input characteristic value), the requirements and boundary conditions for the thermal simulation are described in Chapter 8.4. For ventilation, a distinction is made between a basic air exchange and the increased air exchanges during day and night. The basic air exchange rate is constant and only depends on the type of building.
The increased daytime air exchange rate allows the air exchange rate to be increased to a value of up to 3 1/h under the following conditions: • During occupancy hours • Room air temperature > 23°C This shall mimic either user behaviour of opening a window or of a mechanical ventilation system that increases flow rate. For the increased night air exchange rate, the air exchange rate can be increased to a maximum of 5 1/h if: • Possibility of night-time ventilation must be provided (system, user) • Room air temperature > Heating setpoint temperature • Indoor air temperature < outdoor air temperature The following equations summarize the zonal energy balance and ventilation model, respectively:. The dependence of the air exchange rate on the room air temperature creates a strong dependency and can lead to numerical problems if the ventilation is changed abruptly (see discussion below). The standard itself says nothing about the type of adjustment (step change, or gradual change), or duration of an individual change of air change rate. Given the availability and predominant use of fixed step simulation tools such as EnergyPlus, TRNSYS, TAS etc. at the time of writing of the standard, it can be speculated, that authors only considered a step-2 E3S Web of Conferences 1 0 (2020) 72, 6004 NSB 2020 wise change of air change rate in fixed step integrators. Interestingly, though, no comment is made in the standard on the minimum/maximum length of fixed steps used in such simulators.

Naïve direct control model implementation
In a first attempt, the defined control strategy was implemented directly into the ventilation model within the simulation program THERAKLES. This single-zone simulation model uses an implicit time integration method with modified Newton iteration. The control equations are thus evaluated in every Newton iteration. Whenever room air temperatures are iteratively updated, the ventilation rates may need to be updated, also. This in turn leads to changes of the air change related energy source in the zone. The direct non-linear feedback between energy density (and thus temperature) and the energy source due to air change is at the heart of all the numerical problems discussed below.
To test the implementation variants, a simple test room was designed where increased air flow was enabled during daytime. The base air flow rate is 0.5 1/h, increased air flow rate is 3.0 1/h during daytime. The other parameters of the test room are as follows: • 2 outer walls (U-value 0.29W/m 2 K, 15m 2 ) to north and south with each one 12m 2 window (U-value 2.8W/m 2 K) • 2 outer walls (U-value 0.29W/m 2 K, 30m 2 ) to east and west with each one 25m 2 window (U-value 2.8W/m 2 K) • All windows with automatic shading system made of blinds with shading of 60% • Ground floor 50m 2 with constant ground temperature of 10°C with U-value 0.5W/m 2 K • Flat roof 50m 2 with U-value 0.2W/m 2 K • Residential building according DIN 4108-2 • Constant inner load 4.1667Wh/m 2 d The parameters of the test room are not important in the following because the focus of the paper is on the numerical modelling of ventilation control.
When running the simulation, we observe that already at the first time the room air temperature exceeds the 23 °C limit, the simulation becomes extremely slow. A closer look reveals, that the simulation progresses with tiny time steps of about 10 -5 s, and the room air temperature oscillates between 22.998 °C and 23.001 °C. According to the control scheme defined in the standard, the ventilation rate toggles between 0.5 and 3 1/h from step to step. Since the outside air temperature is well below 15 °C at this simulation time, the corresponding change of energy sink in the room air energy balance is significant. Clearly, such rapid changes cause the error test to fail for all but the tiniest time steps.

Mathematical analysis
The formulated control scheme has three distinct characteristics: • The criterion switching between daytime/occupancy and night control method (i.e. potentially different increased air flow rates) causes a weak discontinuity to arise at the switching time points. This is also called a time-event.
• The criterion to toggle increased ventilation rate on/off at a given room air temperature limit results in a weak discontinuity, when the corresponding energy source/sink term changes. This is also called a state-event.
• The criterion to toggle increased ventilation rate on/off when the ambient temperature reaches/exceeds the room air temperature, results in a weak discontinuity only in the second derivative of the energy density (see Fig. 2c for an illustration).
With weak discontinuity we mean that the state variable itself (energy density/temperature) remains continuous, yet the first time derivative of the conserved quantity changes abruptly from one value to another (see Fig. 2a and 2b). Error controlled time integrators typically handle discontinuities by determining the exact time point at which the discontinuity occurs. For state-events, this involves an expensive root-finding procedure. If one were to reproduce this problem with Modelica and an integrator engine such as DASSL, the solver would stall due to excessive root finding work.
If only weak discontinuities arise in a model formulation, these expensive root-finding methods can 3 E3S Web of Conferences 1 0 (2020) 72, 6004 NSB 2020 usually be avoided. Integrators such as CVODE use the aforementioned error testing procedure to detect such weak discontinuities. For example, if the first derivative of a variable were only slowly changing, the predictor and corrector solutions would match very closely. However, if a trigger point (in our case the temperature limit) is reached in the corrector step, and the time derivative changes abruptly, the corrector solution differs significantly from the predictor and an error test failure is triggered (see Fig. 3). Obviously, a flip-flop type control algorithm, such as the discussed ventilation model, will trigger such error test failures at each step where the gradient changes abruptly -resulting in the observed temperate oscillations and toggled air change rates, and ultimately lead to tiny integration steps.

Time delay approach
Considering the source of the numerical problems, it may be advisable to reduce the frequency of ventilation mode changes. Motivated also by observed occupant behaviour, a time delay could be formulated. After a window was opened for increased ventilation, it is likely to stay open for at least a few minutes, before it is closed again. Hence, a minimum time delay of at least 5 min could be enforced in the simulation, before the air flow rate may be changed again.

Hysteresis/dead band controllers
Another alternative to avoid the constant on/off switching of increased air flow rate may be the use of dead-band control algorithms. Hereby, increased air flow rates are enabled when the indoor air temperature limit of 23 °C is exceeded (likely resulting in a few convergence and error test failures). However, the room air temperature must than first fall below another threshold, say 21 °C, before the base ventilation rate is used again. Actually, this mimics pretty much the behaviour of human occupants, which will usually wait a little and only close the window when they feel too cold.
We still expect error test failures whenever a switch occurs, but if this only happens every 5 minutes or after temperature has dropped somewhat, the simulation will eventually finish. We can further avoid convergence errors by taking the adjustment of the ventilation rate out of the actual Newton iteration. With respect to the time delay variant, we expect that with increasing delay times, the number of control mode switches reduces, as does the error fail count. Thus, simulations with longer delay times should be faster, though potentially less accurate.

Step-wise state adjustment
The implementation of time delays or dead-band controllers typically requires the use of an additional transient equation in the equation system, to monitor the state of the ventilation system. However, we find it more convenient to adjust the ventilation state discretely after a completed integration step. So, instead of evaluating the ventilation rate control rule inside the Newton iteration, and potentially causing convergence errors, we only modify the physical behaviour in between steps.

Interpolation/step smoothing approach
Another common procedure to avoid discontinuous signals in dynamic equations involves ramping or smoothing transition curves. In fact, it would be possible to define a temperature interval, say between 23 and 24 °C, where the air change rate increases linearly from the base ventilation rate to the increased ventilation rate (or in a polynomial or otherwise smoothened shape). While this certainly helps with the numerical problem, it stands to reason what physical interpretation such a ramping implies.
Given the simplicity of the well-mixed zonal air node model, and the equally simple way to treat ventilation (regardless of actual temperature differences, stack pressures, room and window geometry, and air mixing conditions, etc.), it is probably sufficient to simply consider the overall impact of ventilation on the energy balance. Hence, it is probably not meaningful to interpret an interpolated air change rate of say 1.7 1/h as "partially opened window" (though there are people who claim to have seen half opened windows).
It is interesting to note, that the overall energy loss of the zonal energy balance computed with this interpolated approach is very similar to the energy loss in the naïve implementation (if that would ever finish the calculation) or time delay variants. Effectively, the switching of air change rates at a high frequency looks like the output of a pulse-width-modulation (PWM) controller, which generates also an on average reduced signal (see Fig. 4).

Variant analysis
We will now test the described algorithms with a typical application case: climatic data according DIN 4108-2 region C (is equal to German TRY 2011 Nr. 12 -Mannheim) [11], a room with large windows, moderate thermal storage capacity and potential to be ventilated properly described in chapter 2.1. Fortunately, all model variants finished calculated rather quickly, with surprisingly small differences in results and little difference in simulation time.

Detailed look at a warm day
First, we will evaluate the time delay variants. A timer variable was added to the solver that was updated whenever the air flow rate was toggled between base and increased rate. Then, a minimum delay time had to pass, before the air change rate could be changed again.  As expected, longer delay times mean less frequent switching. In all variants the increased ventilation is toggled on at the same time instant, yet for the smaller time delays it is turned off at least once again during the morning hours. At this particular day the hourly time delay variant only increases ventilation rate once. However, in other days some switching is also observed with the hourly time delay variant.
Next, the dead-band control variant is tested. Hereby, we select the dead-band such, that the requirements of the building code are fulfilled, i.e. the increased ventilation is turned off below 23 °C. It is turned on, when the temperature exceeded 23 °C by the selected dead-band size. Fig. 6 shows the results obtained for 0.2 K and 1.0 K dead-band width, respectively. It would also be possible to shift the dead-band below 23 °C (or centred around 23 °C), and thus compute less overheating hours, but we selected the more conservative dead-band location in our tests.  Note, that the variant with the small dead-band requires many air change rate adjustments, and requires much more computational work.
Finally, we look at the variant with linear interpolation. Also for this variant, we can tune the control model by selecting different interpolation interval lengths. Fig. 7 shows the results obtained with a 1.0 K interpolation interval. Starting from 23 °C up to 24 °C the air change rate is linearly changed from base to increased rate. In this variant, we see a gradual increase of air change rate, which only depends on the room temperature.

Figures 8 show a comparison of all different room air
temperatures, calculated with the described control models. Fortunately, the temperature curves look more or less the same, with the exception of the actual switching periods. This period is shown in detail in Fig. 9, where the effect of the individual control strategies on the room air temperature can be distinguished. Particularly interesting is the 0.2 K dead-band variant, where switching occurs every few seconds. Naturally, we expect the longest simulation time for this variant.

Annual analysis and standard compliance
With respect to the anticipated application as standard compliant, simulation-based evaluation of thermal comfort in summer, the most critical question is: how much is the actual compliance criterion, the over-heating hours, influenced by the choice of control algorithm. Surprisingly little, as can be seen in Fig. 10 and Table 2. The simulation duration, however, differs quite a bit among the different variants, see Table 2. The numbers look similar for other room constructions and climatic locations. So we can take the analysed example as somewhat representative.  Given the very similar calculation results, it really doesn't matter which approach is taken when implementing such a control model in an adaptive time step integration solver. While the linear interpolation variant is surely the easiest one to implement (and also works in Modelicabased models, where step-completed model adjustments are not easily possible), we favour the delay time variant.
In general, the described ways of handling discrete control signals in modern error controlled solvers, can be applied similarly to other models.