On the Procedure of Draught Rate Assessment in Indoor Spaces

Featured Application: In case the ventilation system design leads to low-frequency high-amplitude air velocity pulsations, strong time variation of the draught rate (DR) index may exist and the thermal environment categorization requires for space-distributed long-term computations and measurements. Abstract: The objective of the paper is to demonstrate the importance of the unsteady Computational Fluid Dynamics (CFD) simulations and long-term measurements for the reliable assessment of thermal comfort indoors, for proper categorization of the indoor thermal environment and for identifying the reasons for complaints due to draught discomfort. Numerical simulations and experimental measurements were applied in combination to study ventilation in a ﬁeld laboratory, a university classroom with a controlled indoor environment. Strong unsteadiness of the airﬂow was registered both in the unsteady RANS results and the real-scale long-term velocity data measured with thermo anemometer. Low-frequency high-amplitude velocity ﬂuctuations observed lead to substantial time variation of the draught rate. In case of categorization of a thermal environment, the point measurements or steady-state RANS computations would lead to wrong conclusions as well as they cannot be used for identiﬁcation of the reasons for people’s complaints due to draught discomfort if strong unsteadiness of the airﬂow exists. It is demonstrated that the length of the time interval for draught rate (DR) assessment may not be universal if low-frequency high-amplitude pulsations are present in the room airﬂow.


Introduction
The indoor environmental quality (IEQ) depends on both indoor air quality (IAQ) and human comfort (thermal, visual, acoustic). A poor IEQ could precondition different discomforts and health problems that result in short or and long-term absence from work and decreased productivity [1][2][3].
Olesen et al. [4] reported that in more than 40% of enclosures people complain from the comfort or health-related issues. The situation is worsening in schools, where the occupation density is 1.8-2.4 m 2 /person, compared to 10 m 2 /person in offices [5]. Improving the working environment and The air is supplied to the room by the four-way type MMU-AP0092H fan coils that could be used either in an isothermal regime or for cooling/heating purposes. Four fan coils are installed in the room in the staggered arrangement to provide various conditions for ventilation tests. The size of a fan coil is 1.0 × 1.0 × 0.25 m. It has four peripheral supply sections with the size of lin × win = 0.5 × 0.07 m equipped with turning blades (louvers) and one central suction section with the size of 0.5 × 0.5 m (Figure 1b-d). The supply sections are marked as cardinal points: north, east, south, and west ( Figure 1b). The louver at each supply section has five preset angles of inclination, θ, to the ceiling, namely 25°, 30°, 35°, 40°, and 45°. Figure 1d demonstrates the louver inclined at the extreme position of 45°. In the present study, two values of θ were considered, θ = 25° and θ = 45°, and it is one of the varied parameters that specifies the list of cases considered, see Table 1. For each fan coil the fan speed could vary in the range up to 300 rpm, the maximum fan speed was set for the present experiments.

Description of Measurement Technique and Equipment
The experimental program comprised two sets of measurements. In each experiment simultaneous measurements of air velocity and air temperature at four points were performed. An 8-channel thermo anemometer composed of a multichannel power supply module SENSOR HT-430 and four thermo anemometer transducers SENSOR HT-428, each equipped with a SENSOR HT-412 (omnidirectional probe sensitive to both air velocity and air temperature) were used. The HT-412 probes meet all requirements of ISO 7726:1998 [38] standard for class C (comfort)-category desirable, i.e., accuracy better than ± (0.02 m/s + 7% of readings) and response time (90%) better than 0.2 s. The The air is supplied to the room by the four-way type MMU-AP0092H fan coils that could be used either in an isothermal regime or for cooling/heating purposes. Four fan coils are installed in the room in the staggered arrangement to provide various conditions for ventilation tests. The size of a fan coil is 1.0 × 1.0 × 0.25 m. It has four peripheral supply sections with the size of l in × w in = 0.5 × 0.07 m equipped with turning blades (louvers) and one central suction section with the size of 0.5 × 0.5 m (Figure 1b-d). The supply sections are marked as cardinal points: north, east, south, and west ( Figure 1b). The louver at each supply section has five preset angles of inclination, θ, to the ceiling, namely 25 • , 30 • , 35 • , 40 • , and 45 • . Figure 1d demonstrates the louver inclined at the extreme position of 45 • . In the present study, two values of θ were considered, θ = 25 • and θ = 45 • , and it is one of the varied parameters that specifies the list of cases considered, see Table 1. For each fan coil the fan speed could vary in the range up to 300 rpm, the maximum fan speed was set for the present experiments.

Description of Measurement Technique and Equipment
The experimental program comprised two sets of measurements. In each experiment simultaneous measurements of air velocity and air temperature at four points were performed. An 8-channel thermo anemometer composed of a multichannel power supply module SENSOR HT-430 and four thermo anemometer transducers SENSOR HT-428, each equipped with a SENSOR HT-412 (omnidirectional probe sensitive to both air velocity and air temperature) were used. The HT-412 probes meet all requirements of ISO 7726:1998 [38] standard for class C (comfort)-category desirable, i.e., accuracy Appl. Sci. 2020, 10, 5036 5 of 20 better than ± (0.02 m/s + 7% of readings) and response time (90%) better than 0.2 s. The thermo anemometers were calibrated for two measurement ranges for better accuracy of the polynomial fit. The first range was between 0.03 m/s and 2.65 m/s, which corresponds to room airflow velocities. The second range was between 2 m/s and 5 m/s, which corresponds to the velocities of fan coil jets. The mean absolute error (MAE) and the normalized mean bias error (NMBE) were calculated for each sensor and both measurement ranges. The maximum sensor MAE was 0.01 m/s for measurement range 1, and 0.03 m/s for measurement range 2. It roughly corresponds to 1% normalized mean absolute error for both ranges. There was a small positive bias for the sensors in range 1 with normalized mean bias error between 0.02% and 0.05%. For range 2, the bias was practically zero with NMBE in the range ± 0.002%. The 8 voltage signals (4 proportional to air velocity and 4 proportional to air temperature) generated by the HT-412 probes were simultaneously read and stored on the hard drive of a PC by a data acquisition system manufactured by DATAQ Instruments. This system was composed of a DI-2108-P high-speed DAQ board controlled by the WinDaq recording and playback software.
The first measurement set was devoted to the characterization of airflow near the fan coil exhausts. These measurements aimed to evaluate the level of uniformity in the supply flow rates and to get data for inlet boundary conditions to be set in the CFD simulations. All fan coils were set to isothermal recirculation mode with the louver inclination of θ = 45 • (case #3 in Table 1, with longer time sample). For each fan coil, the measurements of air velocity and air temperature at the middle points of four supply sections were performed simultaneously. The sensors were placed below the middle point of each opening at the distance of 0.017 m from the surface (Figure 1d). The total measuring period for each fan coil was 35 min (2100 s) long, with the sampling rate of 5 Hz. The results of this measurement set are discussed in Section 3.1.
In addition to the local measurements, at the same fan coil operation conditions the measurements of the volume flow rate Q, at each supply section of each fan coil were performed by an ALNOR Balometer. Table 2 summarizes the Q values. Table 2. Air flow rate, Q, m 3 /h, through the supply openings of each fan coil.

Fan coil Number
East Section South Section West Section North Section   #1  212  146  190  196  #2  174  158  203  191  #3  170  163  197  168  #4  164  189  193  173 The second set of measurements was devoted to the characterization of the airflow in the room occupied zone. Air velocity and air temperature data were collected simultaneously at four points placed one above another along a vertical line at coordinates z = 0.1 m, 0.6 m, 1.1 m, 1.7 m. According to ISO 7730:2005 [11], z = 0.1 m corresponds to the ankle level of a person, z = 0.6 m corresponds to the level of body center of gravity for a seated person, z = 1.1 m is the head level of a seated person and the body center of gravity level of a standing person, and z = 1.7 m is the head level of a standing person. The line was placed at x = 1.9 m, y = 7.4 m being the intersection of two vertical planes: the plane passing through the middle of the north supply section of fan coil #4 and the plane passing through the middle of the east supply section of fan coil #3. Monitoring points are located along this vertical line since there the supply jets from three fan coils meet -north jet of fan coil #4, the east jet of fan coil #3, and south jet of fan coil #2. The reason for selecting this vertical line is to compare velocity fluctuations parameters in the above-specified jets with the velocity fluctuation parameters in the monitoring points. Based on this a conclusion may be drawn about the reasons of velocity fluctuations at the monitoring points. Measurement point positions will be illustrated below in Section 4.1.
The measurements in the occupied zone were performed for three cases presented in Table 1. Various time intervals used for measurement data collection are given in Table 1. Both the louver inclination angle value and the fan coil operation mode varied from case to case. Isothermal recirculation mode with no outdoor air supply and no cooling/heating was provided for cases #2 and #3, while case #1 was non-isothermal. To avoid direct sunlight influence, all the measurements were done in the afternoon. However, it was not possible to keep the room temperature constant during a long measurement period. Thus, in case #1 during the whole period of 3600 s, the fan coils have been working in the cooling mode to keep an appropriate moderate temperature level in the room.

Computational Aspects
Two sets of CFD simulations of turbulent airflow in the room were performed. The aim of the first set was to reproduce the experimental conditions as much as possible. It allowed collecting the detailed data on the airflow in the room and draught rate (DR) indices. It also allowed validating the numerical data using the experimental data. The turbulent airflow in the first set of computations was modelled using the unsteady Reynolds-Averaged Navies-Stokes (URANS) approach with a fine computational grid; all geometry and boundary conditions were set as close to the experimental conditions as possible. Particular computational aspects for these CFD simulations are described in detail in Section 2.3.1.
The main goal of the second set of the CFD simulations was to demonstrate the possibilities of the up-to-date engineering CFD tools (the steady-state RANS) with a coarse computational grid to predict the DR indices. Keeping, in general, the room geometry, the problem was set with non-isothermal conditions and took into account the effects the room equipment, such as desks and chairs, as well as the occupants. Details of the problem formulation for the second set of computations are given in Section 2.3.2.
As it was mentioned in the Introduction, the RANS approach applied in the current study for turbulence modelling depends on the particular turbulence model used, and the uncertainty due to the turbulence model influence could be high. However, for multi-jet configuration considered in the study, it is not possible to use less-empirical scale-resolving simulation due to enormous computational resources required. From the literature data [39,40], it is evident that a k-ε turbulence model and the baseline SST k-ω model accurately predict jet flows and qualitatively reproduce separation zones behavior that gave a basis for the turbulence model choice in the current computations.

Unsteady RANS (URANS) Simulations for Computational Fluid Dynamics (CFD) Validation
Air was assumed as an incompressible fluid with constant physical properties (ρ = 1.2 kg/m 3 , µ = 1.8 × 10 −5 kg/m·s). The boundary conditions at the supply sections were based on the supply flow rate values, measured in the experiments, and the effective cross-sectional area of the supply section. The mean value of volume flow rate entering the room through one supply section (averaged in time and overall 16 supply sections, Table 2) was equal to 190 m 3 /h (mass flow rate of 0.0647 kg/s). The area of the supply section was decreased as compared with the l in value given in Section 2.1. The cross-section of the supply channel has smaller length then l in (Figure 1d); the effective opening length l in, effective = 0.177 m was set in the URANS computations that gave the mean inlet velocity value of U bulk = 4.26 m/s (it corresponds to the local velocity values detected in the first set of measurements).
The uniform distribution of U bulk was set at each supply section of each fan coil (being an inlet for the computational domain). The Reynolds number calculated using the supply section width was equal to Re = ρw in U bulk /µ = 2 × 10 4 . Two cases with different airflow angle values, θ = 25 • and θ = 45 • , were considered (they correspond to the fixed positions of the louver during the measurements, cases #2 and #3, Table 1). At the outlet boundaries (i.e., at each square suction section of each fan coil), the pressure outlet boundary condition was set assuming uniform pressure distribution. A no-slip boundary condition was set at the solid walls.
A quasi-structured grid consisted of hexahedral cells was built with the grid generator ANSYS ICEM CFD 18.2. The total cell number was 33 million cells; the cells with the smallest size were located in the corners of the inlet and outlet sections: the dimensions of these cells were 3 × 3 × 3 mm. The cells adjacent to the solid walls had the same height of 20 mm (the only exception was the solid wall at the surface of the fan coils). Depending on the local flow features, this height corresponds to the normalized distance from the center of the first near-wall cell to the wall, y + , in the range from the values less than one up to values exceeding 60. The sensitivity analyses of different computational parameters influence on the CFD results were performed previously and are discussed in [37]. In particular, the results of the mesh sensitivity analysis were presented in detail in [37].
The turbulence modelling of the airflow in the room was performed using the URANS approach coupled with the standard k-ε turbulence model in combination with the enhanced wall treatment. To specify the inlet boundary conditions for the turbulence characteristics the ratio of the turbulent to molecular viscosity of ν turb /ν = 10 and the turbulent intensity of 5% were set at each supply section.
The numerical simulations were performed using the CFD package ANSYS Fluent (version 18.2) based on the finite volume method with the cell-centered variable arrangement. The iterative SIMPLEC method was used to link the continuity and momentum equations. The second-order schemes for the spatial and time discretization were chosen. The value of the time step was 0.5 s. The time step value, ∆t, was varied in [37], from ∆t = 0.004 s to ∆t = 0.5 s, and it was found that the lower frequencies and the time-averaged flow patterns were the same in the ∆t range considered. The computations were performed using the resources of Peter the Great St.Petersburg Polytechnic University Supercomputer Center (scc.spbstu.ru). In total, up to 504 cores were used for the parallel computations.

RANS Simulations
The steady-state RANS simulation was set for the modified geometry model, as compared with the geometry shown in Figure 1a. Namely, the ventilation ducts were removed, and the form of columns and ceiling beams was simplified. The fan coil geometry model was also modified: the fan coils were placed directly at the ceiling, and the height of each fan coil was increased to 0.3 m. The supply opening dimensions were also changed, l in,RANS × w in,RANS = 0.4 × 0.1 m, so that the opening area was 3.2 times larger than in URANS problem formulation. Opposite to the URANS problem formulation (and the experimental conditions) with empty room interior, the room for the RANS simulation contained 8 desks, 24 chairs, and 24 persons of simplified shape in seating positions, and one person standing upright. The detailed geometry model for this case was presented in [41].
The RANS simulation was performed for non-isothermal conditions with heat transfer problem solving. At each surface the constant values of heat flux, q, were set ( Table 3). The mass flow inlet boundary condition was set at each supply opening, with the mass flow rate of 0.227 kg/s. This increased mass flow rate value allowed to keep almost the same value of the bulk velocity as in URANS computations, U bulk = 5 m/s. Two cases with different air flow angle values were considered again, with the flow angle corresponded to θ = 25 • and θ = 45 • . It was assumed that cooled air is supplied to the room, and the constant temperature value of 17.7 • C was set at the inlets. A uniform velocity distribution, U bulk_out = −0.69 m/s, was set at suction (outlet) fan coil sections. The model included also four addition outlets located at the sidewall y = 5.4 m (Figure 1b,c). They were assumed to be open, and the pressure outlet conditions were set there. The grid used in the calculations was built in Gambit 2.4.6 grid generator. The grid consisted of cubic elements and the edge length of each cell was equal to 0.1 m. Thus, a grid suitable for engineering usage was generated, with a small number of cells (approximately 200,000). The chosen coarse grid allowed obtaining a steady-state solution due to high numerical dissipation. The problem was solved in ANSYS Fluent 6.3 using RANS approach connected with k-ω SST turbulence model. The Boussinesq approximation for buoyancy forces was used (according to the thermal balance evaluation the bulk temperature in the room was set to 24.5 • C).

Airflow near the Fan Coil
This section presents the near fan coil velocity experimental data obtained for case #3 (Table 1) with θ = 45 • . Figure 2 is an example of the velocity behavior: the data at four points located near supply sections of fan coil #4 are presented. In addition to the instantaneous velocity data (red curves), the evolution curves of the velocity time-averaged for 60 s, <V> 60 s , are plotted (green curves). Table 4 presents the integral information on the mean velocity near the supply sections of fan coil #4 collected from the whole measurement interval of 2100 s, as well as the lowest and the highest velocity values over the measurement sample. Table 5 presents the corresponding information on the turbulence intensity: where N is the number of measurements of instantaneous local air velocity V i , the summation is from i = 1 to i = N. The values of Tu (Table 5) are averaged over three intervals: 10 s, 30 s, and 60 s (the variables <Tu> 10 s , <Tu> 30 s , and <Tu> 60 s , correspondingly).

Airflow near the Fan Coil
This section presents the near fan coil velocity experimental data obtained for case #3 (Table 1) with θ = 45°. Figure 2 is an example of the velocity behavior: the data at four points located near supply sections of fan coil #4 are presented. In addition to the instantaneous velocity data (red curves), the evolution curves of the velocity time-averaged for 60 s, <V>60 s, are plotted (green curves). Table 4 presents the integral information on the mean velocity near the supply sections of fan coil #4 collected from the whole measurement interval of 2100 s, as well as the lowest and the highest velocity values over the measurement sample. Table 5 presents the corresponding information on the turbulence intensity: Tu = 100% [∑(Vi -<V>) 2 / (N-1)] 0.5 / <V>, where N is the number of measurements of instantaneous local air velocity Vi, the summation is from i = 1 to i = N. The values of Tu (Table 5) are averaged over three intervals: 10 s, 30 s, and 60 s (the variables <Tu>10 s, <Tu>30 s, and <Tu>60 s, correspondingly).  Velocity time variation records presented in Figure 2 and the data in Table 4 clearly show that the near fan coil velocity characteristics vary much from one supply opening to another. The range of variation of air velocity at the east supply section is the smallest one (0.26 m/s) while the largest one is observed at the north supply section (1.82 m/s). Velocity variation range over the supply sections of fan coil #4 is 1.87 m/s, which is comparable with its range of variation at the north supply   Velocity time variation records presented in Figure 2 and the data in Table 4 clearly show that the near fan coil velocity characteristics vary much from one supply opening to another. The range of variation of air velocity at the east supply section is the smallest one (0.26 m/s) while the largest one is observed at the north supply section (1.82 m/s). Velocity variation range over the supply sections of fan coil #4 is 1.87 m/s, which is comparable with its range of variation at the north supply section. The range of variation of mean air velocity at the supply sections of fan coil #4 is 0.69 m/s. The smallest mean air velocity appears at the north supply section (3.42 m/s) and the largest one appears at the west supply section (4.11 m/s). It is evident that the pattern of the jets generated by the supply sections of fan coil #4 is different, which preconditions dissimilar interaction of the supply air jets with the room air along the direction of the jets spread. Table 5 demonstrates that the turbulence intensity also varies significantly from one supply opening to another. Following the velocity variation data, the lowest mean <Tu> 60 s = 0.83% is at the east opening, where the range of variation of velocity is the smallest (0.27 m/s), while the highest mean <Tu> 60 s = 6.00% is at the north opening, where the range of variation of velocity is the largest (1.82 m/s). The <Tu> variation at each supply section is also pronounced. For example, at the north supply section the maximum value of <Tu> 60 s is 7.05%, while the minimum value of <Tu> 60 s is 5.13%. However, with an increase in the length of the interval for Tu assessment the range of variation of Tu at each supply section becomes narrower. Minimum Tu increases and maximum Tu decreases while the mean value increases slightly.
Various time variation patterns of both air velocity and Tu at the inspection points are due to different configuration, shape, and area of the supply channels between the fan and the supply port at the face plate of the fan coil. The measurements have shown that the other three fan coils also demonstrate noticeable variation in the supply opening velocity characteristics. This variation is in accordance with the flow rate variation visible from the data measured independently ( Table 2). Both the flow rate and velocity measurements at the supply openings provided a basis for the CFD problems boundary conditions setting. It could be concluded, however, that the uniform velocity and turbulence intensity distribution adopted in the computations could be treated as the first approach only.

Airflow in the Room
The velocity characteristics in the occupied zone were measured during a one-hour period. However, not the whole sample could be used for analysis, as the case was non-isothermal during the whole period of measurements. Fan coils operated in a cooling mode to reduce the initial temperature and to keep the room air temperature at an approximately constant moderate level.  [23.9, 25.9] • C. Therefore, the assessment of Tu and DR at the inspected locations was done for the temperature quasi-steady-state phase only, i.e., for the period starting from the instant of 1020 s (marked with the blue dashed line in Figure 3).
The time variation of air velocity at the inspected locations during the whole period of observation is presented in Figure 4. Besides the instantaneous velocity data (red curves), Figure 4 presents the evolution of the time-averaged air velocity with the averaging period of 1 min, <V> 60 s (green curves), and the averaging period of 3 min, <V> 180 s (black curves). The blue dashed vertical line marks the start of the temperature quasi steady-state phase of the records.
The temperature quasi-steady-state period could be divided into two sub-periods according to two different phases of the velocity behavior. At z = 1.  (Figure 4b) the plateau is relatively short, and a gradual decrease of the mean air velocity starts at the moment t = 1680 s, while at z = 1.1 m the plateau is longer again, and velocity reduction starts at t = 2400 s. The general conclusion from the measurement results is that in addition to the high-frequency turbulence, low frequency pulsations of the velocity field in the room exist, with a period of ten to twenty minutes or even more. only.

Airflow in the Room
The velocity characteristics in the occupied zone were measured during a one-hour period. However, not the whole sample could be used for analysis, as the case was non-isothermal during the whole period of measurements. Fan coils operated in a cooling mode to reduce the initial temperature and to keep the room air temperature at an approximately constant moderate level. Figure 3 presents the air temperature time variation at the four inspected points. The temperature records at all heights can be divided into 2 phases: the transient period and the quasi-steady-state period. The transient phase takes about 1000 s; the air temperature at the inspected locations drops from 32.44 ± 0.1 °C to 25.16 ± 0.2 °C. During the quasi-steady-state period the air temperature at point z = 0.1 m varies in the interval [24.7, 26.5] °C; at z = 0.6 m the interval is [24.6, 26.4] °C; at z = 1.1 m it is [24.3, 26.2] °C, and at z = 1.7 m it is [23.9, 25.9] °C. Therefore, the assessment of Tu and DR at the inspected locations was done for the temperature quasi-steady-state phase only, i.e., for the period starting from the instant of 1020 s (marked with the blue dashed line in Figure 3).  The time variation of air velocity at the inspected locations during the whole period of observation is presented in Figure 4. Besides the instantaneous velocity data (red curves), Figure 4 presents the evolution of the time-averaged air velocity with the averaging period of 1 min, <V>60 s (green curves), and the averaging period of 3 min, <V>180 s (black curves). The blue dashed vertical line marks the start of the temperature quasi steady-state phase of the records.  (Figure 4b) the plateau is relatively short, and a gradual decrease of the mean air velocity starts at the moment t = 1680 s, while at z = 1.1 m the plateau is longer again, and velocity reduction starts at t = 2400 s. The general conclusion from the measurement results is that in addition to the high-frequency turbulence, low frequency pulsations of the velocity field in the room exist, with a period of ten to twenty minutes or even more.

Global Unsteady Airflow Pattern
The results from the numerical simulation of the turbulent airflow in the room are shown in Figure 5. The vertical section, chosen for post-processing, crosses fan coil #2 (the right part of the plot) and fan coil #4 (the left part of the plot), and demonstrates the flow structure only partially. Even at this particular plot, the interaction between two and even three jets is visible. Another source of unsteadiness is the jet impingement to the sidewalls. Detailed analysis of the airflow pattern evolution through velocity field animation at various sections allows concluding that all jets undergo quasi-periodic low-frequency oscillations. Strong three-dimensional unsteadiness is also detected in the low-velocity occupied zone. Note that the used URANS approach resolves only very large eddies evolution, and all other scales (pulsations at higher frequencies) are modelled.
For the flow angle θ = 45° (case #3, Table 1) the unsteadiness is even more pronounced as the jets are detached from the ceiling and spread freely.

Comparison of CFD and Experimental Data
The CFD velocity magnitude data were collected at the four points for which the experimental data are available. For both the cases, the samples of 7500 s correspond to the self-oscillations with the invariant mean value, i.e., the statistically developed regime. Figure 6 compares the evolution of computed and measured velocity magnitude for case #2, θ = 25°. The entire measured sample of one hour is given in the plot; the computed sample is presented partially. As the URANS technique resolves only low-frequency pulsations, the time-averaged measured velocity data with the averaging period of 1 min, <V>60 s, are used for comparison. This type of measurement data processing filters the high frequencies that are useless for the URANS data validation. The low-frequency pulsations kept in the experimental data after processing are reproduced qualitatively in the URANS computations, though there is some quantitative difference. Thus, at the highest point, at z = 1.7 m (Figure 6d) the period of the oscillations is approximately 300 s in the numerical solution and approximately 200 s in the experimental data. The large amplitude of the measured pulsations at this point is attributed to the non-isothermal effects; if the temperature quasi-steady-state period in the time interval [1200, 2700] is considered, the measured amplitude is lower than the computed one. The amplitude of the pulsations obtained in CFD simulation is also higher than the measured amplitude at the lowest point, at z = 0.1 m (Figure 6a); the satisfactory agreement between the computed and measured frequencies is demonstrated. On the contrary, at the remaining two points higher amplitude of pulsations is obtained in the measurements. The main conclusion is that both the experimental and computational data confirm the low-frequency pulsations in the occupied zone that is very important for the thermal comfort evaluation. The vertical section, chosen for post-processing, crosses fan coil #2 (the right part of the plot) and fan coil #4 (the left part of the plot), and demonstrates the flow structure only partially. Even at this particular plot, the interaction between two and even three jets is visible. Another source of unsteadiness is the jet impingement to the sidewalls. Detailed analysis of the airflow pattern evolution through velocity field animation at various sections allows concluding that all jets undergo quasi-periodic low-frequency oscillations. Strong three-dimensional unsteadiness is also detected in the low-velocity occupied zone. Note that the used URANS approach resolves only very large eddies evolution, and all other scales (pulsations at higher frequencies) are modelled.
For the flow angle θ = 45 • (case #3, Table 1) the unsteadiness is even more pronounced as the jets are detached from the ceiling and spread freely.

Comparison of CFD and Experimental Data
The CFD velocity magnitude data were collected at the four points for which the experimental data are available. For both the cases, the samples of 7500 s correspond to the self-oscillations with the invariant mean value, i.e., the statistically developed regime. Figure 6 compares the evolution of computed and measured velocity magnitude for case #2, θ = 25 • . The entire measured sample of one hour is given in the plot; the computed sample is presented partially. As the URANS technique resolves only low-frequency pulsations, the time-averaged measured velocity data with the averaging period of 1 min, <V> 60 s , are used for comparison. This type of measurement data processing filters the high frequencies that are useless for the URANS data validation. The low-frequency pulsations kept in the experimental data after processing are reproduced qualitatively in the URANS computations, though there is some quantitative difference. Thus, at the highest point, at z = 1.7 m (Figure 6d) the period of the oscillations is approximately 300 s in the numerical solution and approximately 200 s in the experimental data. The large amplitude of the measured pulsations at this point is attributed to the non-isothermal effects; if the temperature quasi-steady-state period in the time interval [1200, 2700] is considered, the measured amplitude is lower than the computed one. The amplitude of the pulsations obtained in CFD simulation is also higher than the measured amplitude at the lowest point, at z = 0.1 m (Figure 6a); the satisfactory agreement between the computed and measured frequencies is demonstrated. On the contrary, at the remaining two points higher amplitude of pulsations is obtained in the measurements. The main conclusion is that both the experimental and computational data confirm the low-frequency pulsations in the occupied zone that is very important for the thermal comfort evaluation. The computed and measured mean velocity data are compared in Figure 7. The computed data, extracted from the time-averaged computational solution, are presented as solid lines covering the whole room height (both the jet zone and the occupied zone). The symbols illustrate all the experimental data available for the occupied zone. The comparison performed for the two cases demonstrates perfect agreement between the measured and computed mean velocity data at the three lower points, z ≤ 1.1 m. The non-isothermal effects influence the values of <V>3600 s, resulting in some increase in the values as compared with <V>360 s (Figure 7a). The largest difference between the measured and computed data is seen at the highest point, z = 1.7 m. This difference could occur due to the point being not far from the jet, especially for case #3 (Figure 7b), and the non-isothermal effects could be more pronounced there. The computed and measured mean velocity data are compared in Figure 7. The computed data, extracted from the time-averaged computational solution, are presented as solid lines covering the whole room height (both the jet zone and the occupied zone). The symbols illustrate all the experimental data available for the occupied zone. The computed and measured mean velocity data are compared in Figure 7. The computed data, extracted from the time-averaged computational solution, are presented as solid lines covering the whole room height (both the jet zone and the occupied zone). The symbols illustrate all the experimental data available for the occupied zone. The comparison performed for the two cases demonstrates perfect agreement between the measured and computed mean velocity data at the three lower points, z ≤ 1.1 m. The non-isothermal effects influence the values of <V>3600 s, resulting in some increase in the values as compared with <V>360 s (Figure 7a). The largest difference between the measured and computed data is seen at the highest point, z = 1.7 m. This difference could occur due to the point being not far from the jet, especially for case #3 (Figure 7b), and the non-isothermal effects could be more pronounced there. The comparison performed for the two cases demonstrates perfect agreement between the measured and computed mean velocity data at the three lower points, z ≤ 1.1 m. The non-isothermal effects influence the values of <V> 3600 s , resulting in some increase in the values as compared with <V> 360 s (Figure 7a). The largest difference between the measured and computed data is seen at the highest point, z = 1.7 m. This difference could occur due to the point being not far from the jet, especially for case #3 (Figure 7b), and the non-isothermal effects could be more pronounced there.

Draught Rate Evaluation Discussion
According to [25], DR values could be calculated if the mean velocity magnitude, <V>, turbulent intensity, Tu, and ambient temperature, T a , are known: 62 (0.37<V>Tu + 3.14). (1) In the present study, three data sets for DR evaluation are available. The first data set is from the RANS solution: the steady-state mean velocity components and Tu coupled with the temperature field. The second data set is from the URANS solution: instantaneous 3D velocity and turbulence kinetic energy fields ready for statistical analysis; URANS study is an isothermal case, and to apply Equation (1) an appropriate constant T a value should be set. The third data set comprises results obtained with local thermo anemometer measurements: the measured values of <V>, Tu, and T a at four points in the occupied zone.

RANS-Based Draught Rate (DR) Evaluation
The calculation of DR based on the results of RANS CFD simulations is a straightforward procedure since all needed quantities (velocity, turbulence intensity, and air temperature) are provided by the solution. Steady-state distributions of DR are given in Figure 8  According to [25], DR values could be calculated if the mean velocity magnitude, <V>, turbulent intensity, Tu, and ambient temperature, Ta, are known: DR = (34 − Ta) (<V> − 0.05) 0.62 (0.37<V>Tu + 3.14). (1) In the present study, three data sets for DR evaluation are available. The first data set is from the RANS solution: the steady-state mean velocity components and Tu coupled with the temperature field. The second data set is from the URANS solution: instantaneous 3D velocity and turbulence kinetic energy fields ready for statistical analysis; URANS study is an isothermal case, and to apply Equation (1) an appropriate constant Ta value should be set. The third data set comprises results obtained with local thermo anemometer measurements: the measured values of <V>, Tu, and Ta at four points in the occupied zone.

RANS-Based Draught Rate (DR) Evaluation
The calculation of DR based on the results of RANS CFD simulations is a straightforward procedure since all needed quantities (velocity, turbulence intensity, and air temperature) are provided by the solution. Steady-state distributions of DR are given in Figure 8   According to the plots, the RANS-obtained DR values correspond, in general, to a small percentage of people, dissatisfied due to draught. However, the room classification differs from one level to another, and there is strong spatial non-uniformity in the distributions due to localization of the steady-state jets. At θ = 25° when the jets are attached to the ceiling (Figure 8a)  According to the plots, the RANS-obtained DR values correspond, in general, to a small percentage of people, dissatisfied due to draught. However, the room classification differs from one level to another, and there is strong spatial non-uniformity in the distributions due to localization of the steady-state jets. At θ = 25 • when the jets are attached to the ceiling (Figure 8a)  There are 16 high-DR spots at z = 1.1 m that correspond to the jet penetration locations, where 25% ≤ DR ≤ 30% there. These spots are more pronounced at z = 1.7 m, with 30% ≤ DR ≤ 35% inside the jets.
The conclusion from the engineering RANS data analysis is that the ventilation scheme used is designed well and provide satisfying DR values. However, the RANS computations provide information on mean components of velocity that define the mean velocity magnitude, or simply "velocity". To evaluate the thermal comfort indices (and DR in particular), the time-averaged velocity magnitude, <V>, called the mean "speed" distributions, must be used [42]. It is known that in case of large velocity fluctuations that is detected for the test classroom, both in the experiments and in the unsteady computations, the "speed" values could be much higher than the "velocity" values [40].

URANS-Based DR Evaluation
URANS simulations are the first step on the way of the velocity fluctuations accounting. After time-averaging, URANS data provide the 3D fields of mean velocity magnitude, <V>, and turbulent kinetic energy, k. Tu could be computed from k using the following relation: Tu = (2<k>/3) 0.5 /<V>. Equation (1) could be applied to the isothermal URANS data after some temperature level is set. The DR distributions computed with the indoor air temperature taken as 27 • C are presented in Figure 9. The DR fields are plotted at the same horizontal cross-sections: z = 0.1 m, z = 0.6 m, z = 1.1 m, and z = 1.7 m. The integral DR value area-averaged over four z-planes under consideration is about 10%, but thermal comfort evaluation based on the averaged DR value is incomplete. First, it does not take into account the DR stratification visible in the occupied zone. Second, it does not take into account strong DR non-uniformity at each section, both spatial and temporal. The latter was not considered in the RANS simulation, but in the URANS fields, per the unsteady behavior of velocity distribution, discussed in Section 4.1, the DR fields are substantially unsteady. There are 16 high-DR spots at z = 1.1 m that correspond to the jet penetration locations, where 25% ≤ DR ≤ 30% there. These spots are more pronounced at z = 1.7 m, with 30% ≤ DR ≤ 35% inside the jets. The conclusion from the engineering RANS data analysis is that the ventilation scheme used is designed well and provide satisfying DR values. However, the RANS computations provide information on mean components of velocity that define the mean velocity magnitude, or simply "velocity". To evaluate the thermal comfort indices (and DR in particular), the time-averaged velocity magnitude, <V>, called the mean "speed" distributions, must be used [42]. It is known that in case of large velocity fluctuations that is detected for the test classroom, both in the experiments and in the unsteady computations, the "speed" values could be much higher than the "velocity" values [40].

URANS-Based DR Evaluation
URANS simulations are the first step on the way of the velocity fluctuations accounting. After time-averaging, URANS data provide the 3D fields of mean velocity magnitude, <V>, and turbulent kinetic energy, k. Tu could be computed from k using the following relation: Tu = (2<k> / 3) 0.5 / <V>. Equation (1) could be applied to the isothermal URANS data after some temperature level is set. The DR distributions computed with the indoor air temperature taken as 27 °C are presented in Figure 9. The DR fields are plotted at the same horizontal cross-sections: z = 0.1 m, z = 0.6 m, z = 1.1 m, and z = 1.7 m. The integral DR value area-averaged over four z-planes under consideration is about 10%, but thermal comfort evaluation based on the averaged DR value is incomplete. First, it does not take into account the DR stratification visible in the occupied zone. Second, it does not take into account strong DR non-uniformity at each section, both spatial and temporal. The latter was not considered in the RANS simulation, but in the URANS fields, per the unsteady behavior of velocity distribution, discussed in Section 4.1, the DR fields are substantially unsteady. Due to pronounced unsteady effects, the instantaneous distributions given in Figure 9 differ from the RANS-obtained fields presented in Figure 8a. In general, the URANS solution reports stronger draft discomfort than the RANS data. Thus, the spots with local values of DR exceeding 20% Due to pronounced unsteady effects, the instantaneous distributions given in Figure 9 differ from the RANS-obtained fields presented in Figure 8a. In general, the URANS solution reports stronger draft discomfort than the RANS data. Thus, the spots with local values of DR exceeding 20% at two z-coordinate levels: z = 1.1 m and z = 1.7 m, change their locations from one instant to another, and the resulting high-DR zones cover a larger area as compared with the RANS data (category III in EN 16798-1:2019 [28]). The DR values exceed 10% at z = 0.1 m and z = 0.6 m (category II in EN 16798-1:2019 [28] instead of category I, resulted from the RANS data).

DR Evaluation Based on the Measurement Data
The measured DR time evolution curves plotted in Figure 10 represent data at four measurement points, processed with three different averaging procedures. The DR values are calculated using three moving intervals: the period of one minute, <DR> 60 s ; the period of three minutes, <DR> 180 s , that corresponds to the practical procedure for evaluation of DR in ventilated spaces presented in [36]; and the period of fifteen minutes, <DR> 900 s , which is in coherence with the experimental protocol, reported in [25]. The data processing started at the instant of t = 1020 s (the blue dashed lines in Figure 4) that corresponds to the beginning of the temperature quasi-steady-state phase. At all z levels considered, the time-averaged DR varies in time. The DR evolution reproduces to some extent the unsteady velocity behavior, presented previously in Figures 4 and 6: characteristic periods of 200-300 s are registered for the two quantities. It is evident that the DR assessment depends much on the averaging period. Remarkably that for the longest averaging period of 900 s, the time-averaged DR values also change in time significantly, see the red curve, especially at z = 0.6 m ( Figure 10b) and z = 1.7 m (Figure 10d).
Appl. Sci. 2020, 10, x FOR PEER REVIEW 15 of 20 at two z-coordinate levels: z = 1.1 m and z = 1.7 m, change their locations from one instant to another, and the resulting high-DR zones cover a larger area as compared with the RANS data (category III in EN 16798-1:2019 [28]). The DR values exceed 10% at z = 0.1 m and z = 0.6 m (category II in EN 16798-1:2019 [28] instead of category I, resulted from the RANS data).

DR Evaluation Based on the Measurement Data
The measured DR time evolution curves plotted in Figure 10 represent data at four measurement points, processed with three different averaging procedures. The DR values are calculated using three moving intervals: the period of one minute, <DR>60 s; the period of three minutes, <DR>180 s, that corresponds to the practical procedure for evaluation of DR in ventilated spaces presented in [36]; and the period of fifteen minutes, <DR>900 s, which is in coherence with the experimental protocol, reported in [25]. The data processing started at the instant of t = 1020 s (the blue dashed lines in Figure  4) that corresponds to the beginning of the temperature quasi-steady-state phase. At all z levels considered, the time-averaged DR varies in time. The DR evolution reproduces to some extent the unsteady velocity behavior, presented previously in Figures 4 and 6: characteristic periods of 200-300 seconds are registered for the two quantities. It is evident that the DR assessment depends much on the averaging period. Remarkably that for the longest averaging period of 900 s, the time-averaged DR values also change in time significantly, see the red curve, especially at z = 0.6 m ( Figure 10b) and z = 1.7 m (Figure 10d). In addition to the time evolution curves, Figure 11 presents the integral range of the DR time variation: for three intervals of averaging (60 s, 180 s, 900 s) the minimum, maximum, and median values over the thermal quasi-steady-state period are given. With the increase in the averaging interval, the range of DR variation becomes narrower. The minimum value grows up, the maximum value decreases, while the mean and the median values remain almost the same for the whole temperature quasi-steady-state period. At each level, the difference between the maximum DR and minimum DR is the lowest for the 900 s long interval; this difference varies from one point to another: at z = 1.7 m this difference is 23.6%, at z = 1.1 m it is 9.1%, at z = 0.6 m it is 19.8%, and at z = 0.1 m it is 6.8%. For the three intervals of assessment, the maximum difference between the DR median values is observed at z = 1.7 m, where the median value of <DR>60 s is 2.0% lower than that of <DR>900 s. In addition to the time evolution curves, Figure 11 presents the integral range of the DR time variation: for three intervals of averaging (60 s, 180 s, 900 s) the minimum, maximum, and median values over the thermal quasi-steady-state period are given. With the increase in the averaging interval, the range of DR variation becomes narrower. The minimum value grows up, the maximum value decreases, while the mean and the median values remain almost the same for the whole temperature quasi-steady-state period. At each level, the difference between the maximum DR and minimum DR is the lowest for the 900 s long interval; this difference varies from one point to another: at z = 1.7 m this difference is 23.6%, at z = 1.1 m it is 9.1%, at z = 0.6 m it is 19.8%, and at z = 0.1 m it is 6.8%. For the three intervals of assessment, the maximum difference between the DR median values is observed at z = 1.7 m, where the median value of <DR> 60 s is 2.0% lower than that of <DR> 900 s . The conclusion is that it is necessary to pay much attention to the measurement period for a reliable DR estimation. If the median values from a set of numerous measurements are used, the DR assessment interval for the room under consideration could be 180 s long since median values for <DR>180 s and <DR>900 s are practically equal. However, for the practical categorization of the room, it is insufficient, as a single measurement of 180 s long could be performed during a high-DR or low-DR period and the input used for categorization would be far from the median. The recommendation for the averaging period choice is that the range of variation of DR, i.e., the difference between the maximum and the minimum <DR> values, for the sufficient averaging period, must be lower than 10% of the median value. Note that among the four points under consideration, at two points this condition is not fulfilled even for the longest period of 900 s (see Figure 11c,d).

Categorization of the Classroom Environment Based on Measured DR
Information for the categorization for DR at the monitored locations is presented in Figure 12. At the monitored location, based, for example, on the measured data at z = 1.1 m, the probability to categorize the thermal environment with respect to DR, using a 60 s long interval for assessment, is as follows: 19.6% as Category I, 23.0% as Category II, and 45.4% as Category III. A probability of 12.0% exists that the categorization evaluated at this point does not meet the requirements of the standard (categorization is attributed to be "Fail"). If 900 s long intervals are used for measurements at the same point, there is no probability of the DR categorization to be assessed as "Fail". This example shows both the importance of the long-term measurements for assessment of the thermal environment with respect to DR and the way of how to cope with the draught complaints.   The conclusion is that it is necessary to pay much attention to the measurement period for a reliable DR estimation. If the median values from a set of numerous measurements are used, the DR assessment interval for the room under consideration could be 180 s long since median values for <DR> 180 s and <DR> 900 s are practically equal. However, for the practical categorization of the room, it is insufficient, as a single measurement of 180 s long could be performed during a high-DR or low-DR period and the input used for categorization would be far from the median. The recommendation for the averaging period choice is that the range of variation of DR, i.e., the difference between the maximum and the minimum <DR> values, for the sufficient averaging period, must be lower than 10% of the median value. Note that among the four points under consideration, at two points this condition is not fulfilled even for the longest period of 900 s (see Figure 11c,d).

Categorization of the Classroom Environment Based on Measured DR
Information for the categorization for DR at the monitored locations is presented in Figure 12. At the monitored location, based, for example, on the measured data at z = 1.1 m, the probability to categorize the thermal environment with respect to DR, using a 60 s long interval for assessment, is as follows: 19.6% as Category I, 23.0% as Category II, and 45.4% as Category III. A probability of 12.0% exists that the categorization evaluated at this point does not meet the requirements of the standard (categorization is attributed to be "Fail"). If 900 s long intervals are used for measurements at the same point, there is no probability of the DR categorization to be assessed as "Fail". This example shows both the importance of the long-term measurements for assessment of the thermal environment with respect to DR and the way of how to cope with the draught complaints.

Conclusions
The paper presents the results of a combined numerical and experimental study of room air movement in a test university classroom with a controlled indoor environment. Both the 3D field results computed with the unsteady RANS technique and the real-scale long-term point velocity data measured with the thermo anemometer report about strong unsteadiness of the airflow. Lowfrequency high-amplitude velocity fluctuations emerging due to airflow instability lead to strong time variation of the draught rate that corresponds to possible uncertainty in the percentage of dissatisfied persons due to unwanted local cooling/heating related to draught.
For large indoor spaces with numerous interacting air supply jets, it is recommended to assess DR based on time records of both air velocity and air temperature not shorter than half an hour. The categorization procedure must be done for a sufficiently long interval. The recommendation/guidance of such kind must be given in the standards since not all people dealing with categorization and assessment of the indoor thermal environment are professionals in the field of fluid mechanics. However, as it is not possible to specify a universal length of the time interval for the DR assessment procedure in case of a low-frequency unsteadiness occurs, a preliminary study of the DR sensitivity to the averaging period should be done during a categorization procedure.
In case of categorization of the thermal environment, the time-point measurements or steadystate RANS computations would lead to wrong conclusions as well as they cannot be used for identification of the reasons for complaints due to draught discomfort, if strong unsteadiness of the airflow is expected. In general, the presented analyses demonstrated the importance of the combination of the unsteady computer simulations and long-term measurements for the correct assessment of both DR and thermal comfort, as well as for the proper categorization of the indoor thermal environment.
The general recommendation is to perform measurements and simulations under design conditions, with respect to both occupancy and mode of operation. If geometry details and occupants are modelled in an engineering study, the grid has to be coarse to reduce the computational cost. Even in this case, the CFD problem formulation must be unsteady, and the turbulence model and the quality of the mesh must be sufficient to reproduce the main low-frequency pulsations. A study of

Conclusions
The paper presents the results of a combined numerical and experimental study of room air movement in a test university classroom with a controlled indoor environment. Both the 3D field results computed with the unsteady RANS technique and the real-scale long-term point velocity data measured with the thermo anemometer report about strong unsteadiness of the airflow. Low-frequency high-amplitude velocity fluctuations emerging due to airflow instability lead to strong time variation of the draught rate that corresponds to possible uncertainty in the percentage of dissatisfied persons due to unwanted local cooling/heating related to draught.
For large indoor spaces with numerous interacting air supply jets, it is recommended to assess DR based on time records of both air velocity and air temperature not shorter than half an hour. The categorization procedure must be done for a sufficiently long interval. The recommendation/ guidance of such kind must be given in the standards since not all people dealing with categorization and assessment of the indoor thermal environment are professionals in the field of fluid mechanics. However, as it is not possible to specify a universal length of the time interval for the DR assessment procedure in case of a low-frequency unsteadiness occurs, a preliminary study of the DR sensitivity to the averaging period should be done during a categorization procedure.
In case of categorization of the thermal environment, the time-point measurements or steady-state RANS computations would lead to wrong conclusions as well as they cannot be used for identification of the reasons for complaints due to draught discomfort, if strong unsteadiness of the airflow is expected. In general, the presented analyses demonstrated the importance of the combination of the unsteady computer simulations and long-term measurements for the correct assessment of both DR and thermal comfort, as well as for the proper categorization of the indoor thermal environment.
The general recommendation is to perform measurements and simulations under design conditions, with respect to both occupancy and mode of operation. If geometry details and occupants are modelled in an engineering study, the grid has to be coarse to reduce the computational cost. Even in this case, the CFD problem formulation must be unsteady, and the turbulence model and the quality of the mesh must be sufficient to reproduce the main low-frequency pulsations. A study of the DR variation effects in case of an indoor space is filled with the furniture and occupants is planned as future work.
To summarize, the DR assessment must be done for a sufficiently long interval allowing identification of the velocity field variations and the DR variation in the interval for assessment should be as small as possible. However, even if these conditions are fulfilled, the indoor space categorization could be uncertain, e.g., the indoor space could be related to different categories for successive relatively long periods of time. It is necessary to investigate what happens with the occupant's complaints if categorization in an indoor space varies in time. Thus, additional research is needed on how people perceive the thermal environment under the low-frequency unsteadiness. Analysis of the possible complaints about the thermal environment that cannot be discovered by local short-time measurements is also a field for future work.