Unsteady RANS Simulation of Air Distribution in a Ventilated Classroom with Numerous Jets

. The study is devoted to the Unsteady Reynolds-Averaged Navies-Stokes (URANS) simulation of ventilation in an isothermal room with numerous jets supplied from ceiling diffusers. The computations of the airflow under the test conditions considered were carried out in the classroom of the Technical University of Sofia with no occupants. The room floor has a simple rectangular form, but several columns, beams, window sills, and four radiators are located inside the room that makes the geometry more complex. Air is supplied to the room through four ceiling fan coils, the Reynolds number is 2×10 4 . Calculations were carried out using the ANSYS Fluent 18.2 software with the standard k -  turbulence model chosen. Computational meshes of up to 33 million hexahedral cells clustered to the inlet and outlet sections were used. The main aim of the study presented is to analyze and discuss the complicated 3D flow structure in the room and to give foundation for future measurements of air velocity field in the


Introduction
The Indoor Air Quality (IAQ) control is of key importance in the design of Heating, Ventilation and Air Conditioning (HVAC) systems. The requirements for the HVAC systems design are based on keeping a certain level of air speed and temperature in the ventilated space, as well as contaminant source control and efficient contaminant removal. There are several recognized standards for ventilation system design and acceptable IAQ, such as American Society of Heating, Refrigerating and Air-Conditioning Engineers (ASHRAE) standard [1] or European standard [2]. To satisfy the IAQ standard, it is necessary, in particular, to define the required level of ventilation in air changes per hour or the outside air supply rate that is sufficient to prevent formation of large stagnant zones, to keep the temperature values within the prescribed range, and to limit the air pollution level.
To comply with these requirements from integral point of view, simplified lumped-parameters approaches are applied successfully for design and optimization of HVAC systems (see, e.g., [3]). These approaches use simplified analytical and empirical models and provide information on mass-averaged parameters in the ventilated room. Along with these methods, Computational Fluid Dynamics (CFD) simulation has become widely used in HVAC industry. Started in late nineties, the CFD techniques application has been continuing to rise during the last two decades. Opposite to simplified approaches, CFD modeling of turbulent indoor airflow allows to produce detailed threedimensional (3D) field information on air quality, e.g., size and location of stagnation and high-velocity zones, peculiarities of temperature distributions with respect to gradient values and location, as well as spatial and temporal information on contaminant transport.
In engineering practice, CFD models based on the Reynolds-Averaged Navier-Stokes (RANS) equations solution are widely used due to moderate computational resources required. RANS approach has been applied to solve complex ventilation problems such as design of HVAC systems in airplane cabins [4,5], indoor swimming pools [6,7], a solar house [8], and an ice rink arena [9]. RANS-based CFD studies aimed at development of life support systems on board of the International Space Station (ISS), especially in case of off-nominal operation and other contingency scenarios, are presented, e.g., in [10][11][12]. RANS approach solves a set of transport equations obtained from the Navier-Stokes equations by means of the averaging procedure that results in the unknown Reynolds stress tensor emergence. To close the equations, the Reynolds stresses must be modeled, and the commonly used method is to apply the Boussinesq turbulent viscosity approximation [13]. To define the turbulent viscosity, a semi-empirical turbulence model should be involved. RANS results depend strongly on the particular turbulence model used, and the uncertainty due to the turbulence model influence could be very high, especially when a fully developed turbulent flow, e.g., a turbulent jet, combines with a moderate Reynolds number flow in one problem, that is typical for ventilation tasks [14]. Consequently, RANS solutions demand comprehensive validation. The validation of RANS results could be performed either using more accurate vortex-resolving approaches (Large Eddy Simulation, LES, or Direct Numerical Simulation, DNS) or well-documented benchmark experimental data obtained in different test room configurations.
Vortex-resolving LES approach solves the filtered Navier-Stokes equations. In LES, large scales of motion are resolved directly, while small scales must be modeled with an appropriate subgrid-scale (SGS) model. LES application gives information on instantaneous components of velocity, Vx, Vy and Vz (sure, the term "instantaneous" covers the resolved frequencies only). The overview and outlook of LES could be found in [15]. LES, being an eddy-resolving technique, requires for large computational resources. However, it gives more accurate prediction of turbulent flows, and with the development of the parallel computations becomes more and more popular not only in fundamental studies.
To validate the CFD results, it is necessary to compare the numerical data with the experimental benchmark data. Extended description of most representative mixing ventilation benchmark experiments with the correspondent reference list is given in [16]. Examples of well-documented experimental data are as follows. The most popular isothermal benchmark test is the 2D test by Nielsen et al performed in 1978 [17]. The test data on mean velocity and fluctuations measured with the Laser-Doppler Anemometry (LDA) are available in [17] and on the website http://www.cfd-benchmarks.com/. These data have been used in more than 50 papers for CFD validation during the last two decades. Most of the contributions presented RANS simulation results, and among others, PhD theses by Bennetsen [18] and Voight [19] provided the most complete sets of RANS data under conditions of the test [17]. Recent contributions [20,21] presented accurate RANS and wall-modeled LES (WMLES) data on mixing ventilation in the test room [17].
Another example of isothermal ventilation test flow for more complicated geometry is the configuration considered by Posner et al. in 2003 [22] where LDA and particle image velocimetry (PIV) measurements were performed in a room with a large flow obstruction. Experimental data for this test is limited -the velocity profiles are given only in two sections of the room near the jet zone. Detailed computational study compared the CFD and experimental data was performed in [23,24] using both RANS [23] and LES [24] approaches.
Several tests deal with a single 3D jet supplied from a diffuser to the confined space. Mocikat et al. proposed one of the tests in 2003 [25]. The test room includes a huge obstacle that interrupt the jet propagation. The velocity measurements were obtained using LDA technique with an accuracy of about 1%. RANS simulation of airflow under the test conditions was performed in [26], and pronounced unsteady effects were reported.
In mixing ventilation systems diffusers are often located at side walls, above the inhabitants, and it is necessary to predict correctly not only the jet behavior, but the interaction between the jet zone, characterized by high velocities, and the low-velocity zone of induced secondary flows (occupied zone). Recent experiments by Hurnik et al. in 2015 [16, 27] are the only example of a 3D test for a room with a sidewall jet with pronounced division into the jet and occupied zones. 3D jet supplied from a rectangular sidewall diffuser was studied using the LDA measurements (two air velocity components were measured with the uncertainty lower than 5%), while velocity fields in the occupied zone were measured with the omnidirectional low velocity thermal anemometers (LVTA) techniques. The full data set on velocity measurement results is provided in a special file attached to [27]. WMLES data of the side jet spreading in the confined space was compared with the velocity profiles in the jet zone [28] and in the occupied zone in [29], with special discussion of the mean velocity and mean speed difference.
The current paper presents results of the Unsteady Reynolds-Averaged Navies-Stokes (URANS) simulation of ventilation in an isothermal room with numerous jets supplied from ceiling diffusers. The main goal of the study is to analyze the complicated 3D flow structure in the room with jet interaction. The study will give foundation for future measurements of air velocity field in the room.

Problem formulation 2.1 Geometry model
The study covers airflow in the classroom at the Technical University of Sofia shown in Figure 1. The room has the following dimensions: the width of W = 5.4 m (along the x-axis), the length of L = 11.8 m (along the y-axis), and the height of H = 3.2 m (along the z-axis).

Fig. 1. Classroom at the Technical University of Sofia
The computational domain adopted for numerical simulation of the indoor ventilation problem is shown in Figure 2. The origin of the coordinate systems is located in the floor corner of the room ( Figure 2). The geometry model includes the window sill along one of longitudinal wall, four radiators located under the window sill (r1-r4 in Figure 2), four columns (I, II, III and IV in Figure 2), two ceiling beams placed in the corners of the room along the longitudinal direction, four ceiling diffusers (#1 -#4 in Figure 2) and corresponding ventilation ducts. All the elements are treated as rectangles and/or parallelepipeds, so that, if necessary, the shape of the elements included into the geometry model is simplified slightly.
The dimensions of the elements in the room are as follows. The window sill, with a width of 0.4 m and a height of 0.1 m, is stretched along the side wall x = 0, the width and the height of its base is equal to 0.1×0.9 m. All radiators have the same dimensions of 0.2×1.5×0.5 m, the distance from the side-wall y = 0 to the radiator r1 is 0.7 m, to r2 -3.6 m, to r3 -6.7 m, and to r4 -9.6 m; the distance from the floor to the bottom of each radiator is 0. The room is equipped with four conventional ceiling fan coils (or diffusers, marked as #1 to #4 in Figure 2). The size of each fan coil is 1×1×0.25 m. The distance from the wall x = 0 to fan coils #2 and #4 is equal to 1.4 m, the distance to fan coils #1 and #3 is equal to 2.4 m; the distance from the side-wall y = 0 to the fan coil #1 is 1 m, #2 -3.9 m, #3 -6.9 m, and #4 -9.8 m.
As shown in Figure 3, each fan coil has four peripheral supply sections equipped with turning vanes and one central suction section. The supply sections are inlet openings for the computational domain shown in Figure 2, and the suction section is an outlet opening for the computational domain. The size of the suction section is 0.5×0.5 m. The size of each supply section is equal to 0.5×0.07 m; each supply section is placed at the distance of 0.08 m from the corresponding edge of the suction section. In total, the room includes four outlet openings and sixteen inlet openings.
Air is supplied to the fan coils by means of ventilation ducts shown in Figure 2. The airflow in the ducts is not considered in the current simulation, so that the ducts are just the obstacles for the airflow in the computational domain. The vertical duct with the width of 0.5 m and the length of 0.3 m adjoins the wall y = 11.8 m (it is located at the distance of wvent = 2.6 m from the wall x = 0). The fan coils are attached to the horizontal ventilation ducts with the height of 0.1 m extending along the ceiling.

Airflow parameters
Air was assumed as an incompressible fluid with constant physical properties taken at the temperature 15С (ρ = 1.225 kg/m 3 , μ = 1.810 -5 kg/ms). The uniform velocity distribution over the inlet section area was assumed, and the velocity value of Ubulk = 4.26 m/s was set at each inlet section. The supply mass flow rate of G = 0.183 kg/s was adopted (G = ρ lin win Ubulk), the correspondent volume flow rate is Q = 537.8 m 3 /h. It was assumed that the turning vanes provide the flow angle of 45°. The Reynolds number calculated using the inlet width is equal to Re = ρwin Ubulk /μ = 2×10 4 (if the hydraulic diameter is used as the length scale, Dh = 0.1228 m, the correspondent Reynolds number is ReDh = 3.6×10 4 ). At the outlet boundaries, the pressure outlet boundary condition was set assuming the uniform pressure distribution. The no-slip boundary condition was set at the solid walls.

Turbulence modeling
Turbulent airflow in the room was computed using the Reynolds-Averaged Navies-Stokes approach. According to the literature data, for RANS simulation of ventilation problems, one of the k-ε semi-empirical turbulence models is preferable. The standard k-ε turbulence model coupled with the enhanced wall treatment option was adopted for the current computations. In this model, the turbulent viscosity is defined as νturb = Cµk 2 /ε, where k is turbulent kinetic energy, ε is its rate of dissipation, and Cµ = 0.09 is the model constant. It was not possible to get a steady-state solution, so that the unsteady RANS (URANS) approach was used.
To specify boundary conditions for turbulence characteristics the ratio of the turbulent to molecular viscosity of turb/ = 10 and the turbulent intensity of 5% were set at each inlet section.

Computational meshes
The computational meshes were created with the ANSYS ICEM CFD 18.2 mesh generator. First, three quasi-structured meshes of the same sub-map topology consisted of hexahedral cells were generated. The total cell number for the coarse mesh was 5.9 million cells, for the baseline mesh -10 million cells, and for the refined mesh -33 million cells. All meshes were nonuniform, with clustering to the inlet and outlet sections and to the walls. Figure 3b illustrates the refined mesh at the fan coil surface. The mesh covers the inlet and outlet sections, and strong clustering is visible in the figure.
The cells of the smallest size are located in the corners of the inlet and outlet sections: for the coarse and baseline meshes these cells have dimensions of 10×10×10 mm, while for the refined mesh the dimensions are 3×3×3 mm (see Table 1).
Second, two polyhedral meshes consisted of 4.5 and 7 million mesh elements were generated. The computational domain at a distance from the diffusers was filled with almost equal polyhedral cells, the volume of each cells was equal to Volmax; mesh refinement was applied in the jet zones, and beyond a small transition layer the refined cell volume was equal to Volmin. The characteristics of the mesh cells are shown in Table 1. For the hexahedral meshes, the cells adjacent to the solid walls have the same height of 20 mm (the only exception is 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 lower than one up to values exceeding 40. Figure 4 gives an instantaneous y + -distribution over several walls: zones with peak y + -values correspond to the jet impingement regions at the side walls. Note that the y + -values at the floor of the room are also relatively high, about 20. Polyhedral meshes ensured the y + -values from the same range.
The computations were carried out using the resources of Peter the Great St.Petersburg Polytechnic University supercomputer center (scc.spbstu.ru). The computational resources used included 18 nodes of the Polytechnic RSC Tornado cluster. Each node has two CPUs Intel(R) Xeon(R) E5-2697v3 at 2.60 GHz, 14 cores each CPU. In total, up to 504 cores were used for the parallel computations.

Solver settings
Numerical solutions were obtained with the CFD package ANSYS Fluent 18.2 based on the finite volume method with the cell-centered variable arrangement. The second order scheme provided spatial discretization for convective terms. The second order pressure interpolation was used. Depending on the case, to link the continuity and momentum equations the iterative SIMPLEC or PISO methods were used. The non-iterative time advancement scheme (NITA) based on the fractional step method was chosen if possible to get a stable solution. The second order implicit time integration was used.

Time step sensitivity analysis
The time step value, Δt, was varied: the values of Δt = 0.004 s, Δt = 0.05 s, and Δt = 0.5 s were used. For the coarse and baseline quasi-structured meshes, the value Δt = 0.004 s allowed to provide the Courant number values less than one over the entire computational domain. For these calculations, the NITA scheme was used. It was not possible to get solutions with the NITA method for the cases with Δt = 0.05 s and Δt = 0.5 s, so that an iterative algorithm was applied.
To accumulate representative statistics, the total length of the time sample for each case was up to 15 000 seconds. It was verified that the sample collected is sufficient to collect time-averaged data after the statistically developed regime is reached. Figure 5 illustrates the sensitivity of the unsteady solution to the time step variation. Time evolution of z-velocity is given in Figure 5 Figure 5 where the data for the baseline mesh are presented). However, as the RANS approach is used, the main result of the computations is in the mean quantities. For all three values of time step used, time-averaged flow characteristics are the same. Figure 6 illustrates the sensitivity of the time-averaged solution to the number of mesh cells variation; the case with quasi-structured meshes is shown. Mean velocity magnitude profiles along ten horizontal lines placed at two cross-sections are plotted in the figure. Positions of the lines with respect to the room height are given in Table 2; the range from 0.5 m to 2.75 m covers the entire height of the room. The profiles at the section x = 2.1 m (Figure 6c) illustrate the propagation of four jets issued in the longitudinal direction. The first pair of jets located at 3.9 m ≤ y ≤ 4.9 m is issued from the supply sections of diffuser #2, while the second pair of jets located at 9.8 m ≤ y ≤ 10.8 m is issued from the supply sections of diffuser #4. The profiles at the section y = 4.4 m (Figure 6d) illustrate the propagation of two jets issued in the transversal direction from the supply sections of diffuser #2.  It is evident from the figure that there is strong mesh dependence of the time-averaged solution in the occupied zone, at least at z < 1 m. The reason for that is in the insufficient spatial resolution provided by the coarse and baseline meshes in this region due to lack of computational cells there. Hence, prediction of jet propagation is sensitive to the spatial resolution at some distance from the initial jet region, and the jet decay is over predicted on the coarse and baseline meshes. The clustering of the quasi-structured grid used for the current simulation was mostly near the inlet and outlet sections of the diffusers, and did not cover enough the whole jet spreading at the angle of 45. On the stage of the refined mesh generation, better resolution in the occupied zone was ensured, and more cells were used for the mixing layers reproduction.

Mesh sensitivity analysis
For the polyhedral meshes used, less mesh sensitivity is detected. It is illustrated in Figure 7 where velocity profiles computed with two polyhedral meshes and the finest quasi-structured mesh are compared. The velocity profiles are almost the same for all cross-sections of the room; some differences could be observed at lines e1 and e2 that are located near the room floor.
Solution computed with the mesh of 33 million hexahedral cells could be treated as little sensitive to the mesh size.

Flow structure discussion
In accordance with the conclusions from the mesh sensitivity analysis, the current section covers the results computed with the refined mesh of 33 million hexahedral cells. Figure 8   The air jets completely cover the upper part of the room, spreading almost freely at this region, the interaction between the jets is only at the distance from the diffusers. Despite of the identical arrangement of the supplied openings in all diffusers, the symmetry in the flow structure is not observed, at least in the instantaneous flow field analyzed. For example, detailed analysis of the velocity plot presented in Figure 8 points to different angles of the jet turn at the distance from the diffuser for the jets issued from two opposite supply sections of diffuser #2 in the positive and negative ydirection.
The mean velocity fields point to non-symmetrical distribution of the air jets even in the time-averaged flow structure. For example, it can be seen from Figure 7b (lines b2 and c2) that the velocity profiles in the range of coordinates 0 m < x < 2 m and 2 m < x < 4 m differ much. That means that the propagation and, in particular, the degree of decay of the opposite jets from diffuser #2 is slightly different.
Detailed analysis of the airflow pattern (with the attraction of detailed animations of flow structures) shows that all jets undergo quasi-periodic low-frequency oscillations. Figure 9 shows two instantaneous fields of the velocity magnitude at the same cross-section at two successive time instants. The section chosen for postprocessing crosses diffusers #2 (in the left part of the plot) and #4 (in the right part of the plot). Visible difference between two fields proves that there are intensive oscillations of the jets propagating from diffuser #2; weaker oscillations of the jets from diffuser #4 are also observed. Animations of flow structure in time prove that these oscillations are significantly threedimensional.  Figure 10 illustrates z-velocity time-evolution and corresponding power spectra; data for three monitoring points P1, P2 and P3 are given (the data for the same points were presented in Figure 5); the plots present the entire sample corresponding to the statistically developed self-oscillating regime.
The evolution of z-velocity at point P3 clearly detects the low frequency oscillations that correspond to the pulsations of the jet propagating from diffuser #3 in the negative y-direction (see the location of the point at Figure 5 a,b). It is also visible in figure 10b where the plots of the power spectral density (PSD) vs. frequency are shown for the same monitoring points. The main period of the oscillations is approximately equal to T = 750 s (that corresponds to the leading frequency of 1.33×10 -3 Hz visible in Figure 10b). The conclusion is that to detect this URANS-predicted frequency in experiments it is necessary to provide measurement samples of 30 minutes and even more. Figure 11 shows an instantaneous isosurface of velocity magnitude that reflects position of each jet core. The jets propagate periodically rotating with respect to the initial airflow direction. For example, jet propagating from the diffuser #4 (diffuser located at the left side of Figure 11) in the negative direction of the y-axis. Local high-velocity regions in the occupied zone impose restrictions on the spatial resolution there on the measurement stage. The distribution of the measurement points must be enough to detect the local velocity peaks, and it is recommended to keep the distance between the neighboring points less than 0.2 m.

Conclusions
Results of the Unsteady Reynolds-Averaged Navies-Stokes (URANS) simulation of ventilation in an isothermal room of a simple rectangular form with numerous jets supplied from ceiling diffusers are presented and discussed. Simplified boundary conditions with equal flow rates at each supply opening and uniform velocity distribution at the supply opening surface were set. Calculations of airflow at the Reynolds number of 2×10 4 were carried out using the ANSYS Fluent 18.2 software with the standard k- turbulence model chosen.
Results of time step and mesh sensitivity analyses are presented. It was shown that time-averaged velocity fields do not depend on the time step in the range used. Hexahedral meshes of more than 30 million cells or polyhedral meshes of about 5 million cells provide the same quality of the solution.
The calculations revealed low-frequency large-scale oscillations of jets propagating from the supply diffusers, with typical periods of about 600-800 seconds. Pronounced spatial non-uniformity of the velocity field in the occupied zone was detected. The information on the computed air velocity distributions and their temporal behavior could be useful for air velocity measurements planning and, in particular should be considered on the stage of the sensor positions distribution.