Airflow characteristics under planar opposed ventilation jets in a controlled indoor environment

. Healthy, comfortable and intelligent indoor environment is a key objective in comprehensive well-being. This is also a main target of advanced HVAC-technology. In indoor environments, air distribution plays a major role while providing clean air to occupants. Therefore, investigating ventilation jets is an essential matter. In this study, the main objective was to improve knowledge on numerical modeling and airflow characteristics. In addition, the reliability of modeling methods were investigated. The experiments were carried out in a test room by using omnidirectional anemometers. The planar air jets were supplied below the ceiling from the opposite long-side walls. The exhaust openings were correspondingly near the floor. Isothermal and incompressible viscous airflow was simulated by using RANS, URANS, DES (SST-k- ω LES) and SBES (SST-k- ω - LES) methods. The results show that modeling method has considerable effects on the predicted airflow field. However, the study indicates that correctly implemented numerical modeling may predict well the averaged airflow characteristics. Furthermore, the unsteady simulation allows airflows to fluctuate reasonably. In addition, SBES and DES methods were more sensitive in generating the airflow fluctuations than URANS. A recommendation is to carefully test and choose a modeling method for indoor airflows.


Introduction
Air distribution is one of the major factors for health, comfort and performance [1]. The main purpose of air distribution is to provide clean air for individuals, maintaining thermal conditions, and reducing pollutants from breathing zone. Therefore, airflow interaction is an important matter to investigate. Numerical modeling is a common method while investigating indoor airflows. The simplified mathematical model is discretized, in which the partial differential equations are approximated by a system of linear algebraic equations at discrete locations in space and time [2]. The boundary conditions are then set to define the given flow phenomena, whereas the initial conditions are described for a time-dependent computation. A numerical solution of linear algebraic equations is mainly obtained by using iterative methods. However, experimental methods are necessary for validating numerical modeling. For instance, the standards EN ISO 7726 and EN 13182 provide a basis for measuring instruments of indoor climate.
Turbulence has been a great challenge in numerical modeling which has attracted many studies also in the recent years. Traditionally, Reynolds-averaged Navier-Stokes (RANS) modeling has been applied to simulate a time-averaged result. In RANS, the Reynolds-averaged momentum equations provide nonlinear stress terms that have to be modeled because of insufficient number of equations for the unknown variables. Consequently, e.g. the eddy viscosity models have been developed, in which the turbulent viscosity is added to reduce kinetic energy from the predicted flow field. This is mainly because the turbulence motion can be highly related to inertial and viscous forces. It follows that e.g. the two-equation eddy viscosity models have been widely applied because of their robustness and stability.
The Standard k-ε model, the SST k-ω model and the RNG-k-ε model are common two equation turbulence models. However, more advanced methods have been arising due to increased computing resources. In the Standard k-ε model, the turbulent kinetic energy k and the turbulence dissipation ε are calculated to estimate turbulent viscosity . The standard k-ε model was developed especially for the free shear flows and it is widely used for turbulence modeling. The SST-k-ω model, in turn, is a combination of the k-ε and the k-ω models controlled by a blending function [3]. The model was developed particularly for the boundary layer simulations and strong adverse pressure gradients and separation. In theory, the SST-k-ω model behaves similarly to the k-ε model in the free shear flows and similarly to the k-ω model near the wall boundaries. Furthermore, a limiter for turbulent viscosity and turbulence production has been implemented to prevent abnormal changes in the flow.
The turbulent flow is time-dependent and Reynoldsaveraging dampens the fluctuations. Therefore, an obtained numerical solution is often not very reasonable for the detailed airflow structures. However, RANS simulation can predict easily the airflow trends under the given boundary conditions if the uncertainty in steady simulation has been taken into account. Consequently, when modeling the airflow structures, the unsteady simulations are usually preferred, e.g. the Unsteady Reynolds-averaged Navier-Stokes (URANS), Detached Eddy Simulation (DES), Stress-Blended Eddy Simulation (SBES) or Large Eddy Simulation (LES). In an unsteady simulation, also the time-derivative term is discretized, hence a numerical solution is progressed step by step for every time-step. Therefore, a numerical solution is highly affected by the time-step size, so that the residual levels can remain small.
In the LES simulation, the large-scale flow is resolved and the small-scale flow is modeled. The filtered variables are sort of local averages, e.g. around a cell centre, that can be typically a control volume. Hence, the LES simulation resolves only those portions of turbulence that are greater than the filter width [4]. Consequently, the requirements for dense computational grid may exist, particularly in the boundary layers and in the mixing layers, in which the large gradients occur. A time step size can then be limited by calculating the Courant number in preventing an information loss from computation. However, the requirements for LES simulation can be relieved by implementing hybrid methods, e.g. DES or SBES methods. In those methods, the large eddies are resolved in the middle of the flow field and the boundary layers are modeled by URANS.
In this study, the objective was to investigate air velocity fluctuation and assess reliability of modeling methods on airflow characteristics. The novelty of the study lies in the comparison of SBES, DES, URANS and experiments in a controlled environment. In addition, the fluctuation energy of air velocity has been investigated as a function of frequency.

Methods
Numerical modeling was conducted to simulate airflow structures. The spatial and time discretization was implemented to describe an airflow field. The boundary conditions were set to define flow phenomena, and the initial conditions were determined for time-dependent computation. A structured hexahedral computational grid was used and the gradually increased time-step size was limiting unphysical instabilities at the beginning of the computation. Numerical solution of linear algebraic equations was obtained by using the iterative methods. RANS, URANS, DES and SBES methods were applied to discover airflow characteristics. The validation was carried out by using omnidirectional anemometers. The experimental set-up and the indoor airflow measurements were made in the study by Kandzia and Müller (2016;2018) [6][7].

Test room
The model room is 5 m long, 4 m wide and 3 m high (Fig.  1). The walls are made of wooden plates (10 mm) and polyurethane plates (80 mm). Each interior aluminium object is 5.0 m long, 0.4 m wide and 0.6 m high. The opposed planar jets are supplied at the ceiling zone. The slot diffusers are 0.02 m high and covers the total length of wall. The exhaust openings are correspondingly 0.15 m high near the floor. The x-axis refers to width, y-axis refers to length and z-axis refers to height. The symmetry yz-plane is at x=0 and the corresponding xz-plane is at y=2.5m.

Experimental set-up
The experimental setup contained the air distribution system and the interior objects. The supply air velocity was 2 m/s at the inlet and the supply air temperature was 20 °C. The anemometers were attached on the horizontal measuring pole at the vertical mast of the traverse system. The measuring instrument is shown in Table 1. The boundary conditions provided a flow field, in which two  large-scale eddy zones occurred above the internal objects ( Fig. 2).

Numerical modeling
Fig. 3 shows the computational grid in the middle of the room. The total grid included 10 249 000 nodes and 10 037 790 elements. Thus, the grid density was on average 1.7 × 10 5 nodes per m 3 . The CFD simulations were made by using the ANSYS CFX software. The simulations were performed by applying the CFX finiteelement-based control volume method and implicit pressure-based multigrid coupled solver [9]. Furthermore, the high resolution spatial discretization scheme was used. The temporal discretization scheme was the second order backward Euler.

RANS and LES simulation
Incompressible, isothermal and viscous airflow with constant flow properties was simulated. In RANS modeling, the Reynolds averaged fluid flow is expressed as where u is the velocity, ′ is the fluctuating component, p is the pressure,  is the density, ν is the kinematic viscosity, x is the Cartesian coordinate axis, u i ' u j ' is the Reynolds stress tensor, overbar indicates the timeaveraged variable and subscripts i,j,k denotes the Einstein summation convention. The Boussinesq assumption for the Reynolds stress stensor is defined by where ν t is the turbulent kinematic viscosity and δ ij is the ̅̅̅̅̅ is the turbulent kinetic energy. The given terms of assumption can be rewritten as The momentum equation is then modeled numerically by using the effective viscosity ν eff = ν + ν t that is the sum of molecular and turbulent viscosity. Consequently, the time-averaged momentum equation becomes It follows that RANS simulation resolves only the mean flow and models the turbulence motions. In the Standard k-ε turbulence model, the turbulent viscosity is defined as where k is the turbulent kinetic energy, ε is the turbulence dissipation and is the coefficient, e.g. = 0.09. Furthermore, the transport equations for the turbulent kinetic energy k and the turbulence dissipation ε are determined including several constants. Hence, Eq. (6) shows generally the relationship between the turbulent kinetic energy and the turbulence dissipation.
In the SST-k-ω turbulence model, the corresponding turbulent viscosity is calculated as ( 1 , 2 ) , in which 1 is the constant, is the specific turbulence dissipation rate from transport equation, S is the invariant measure of strain rate and 2 is the blending function [3]. Furthermore, the constants for the k-ε model and the k-ω model are calculated by the blending factor where represent a constant in the SST-k-ω model, 1 is a constant in the k-ω model, 2 is a constant in the k-ε model and 1 is the blending function that is zero for the the k-ε model and one for the k-ω model. In large eddy simulation (LES), the volume averaged filtering is usually made for each equation. Thus, the filtered quantity is a function of space and time. General volume filtering may be expressed as where quantity ̅ is the resolvable filtered variable, V is the filtering volume limiting the resolvable turbulence length scale, and s is the turbulence length scale. While the scale is smaller than the filtering volume, the given equation is zero and this portion of turbulence is modeled by adding turbulent viscosity to the flow, which in turn reduces kinetic energy. The modeling is justified, because turbulent motion approaches an isotropic form at small scales. The incompressible filtered momentum equation becomes where , = i j ̅̅̅̅̅ − ̅ ̅ is the Subgrid-scale stress tensor that is modeled. The Boussinesq approximation is made by where ν sgs is the kinematic viscosity of a subgrid scale model. Consequently, Eq. (4) and Eq. (11) have some similarities, although RANS simulates the time-averaged equations and LES simulates the spatial averaged equations at each time-step. It follows that LES simulation with a control volume filter is a simple and convenient simulation method. However, small time-step size may be necessary to limit the Courant number to an acceptable level at high grid density zones. The Courant number can be defined by where ∆ is the time-step size and ∆ is the size of control volume.

URANS simulation
In URANS, the time-derivative term is discretized and included to the computation. Time-averaging dampen a turbulent fluctuation and, as a consequence, the turbulence spectrum cannot be solved. However, URANS may resolve time variations in the unstable flow. In the kε turbulence model, the scalable wall functions were used [9]. In SST-k-ω turbulence model, the automatic wall treatment was used. In those both models, the high resolution advection scheme and turbulence numerics was computed. In addition, four iteration loops were calculated in each time-step to decrease the residuals at an acceptable level. Moreover, the transient scheme was 2 nd order backward Euler. The indoor air was modeled as an ideal gas. The inlet boundary conditions were constant velocity with medium turbulence intensity (Tu=5%) and medium eddy viscosity ratio (r=10). The outlets were zero relative pressure openings. Other boundaries had isothermal no-slip conditions.

DES simulation
The detached eddy simulation (DES) is a hybrid RANS-LES simulation method. The purpose is to model boundary layer with RANS and switch to the LES in the detached areas based on the grid resolution [4]. In principle, the classification can be made by choosing a limiting length scale e.g. ~3 /2 ⁄ and comparing the given limit with the maximum length scale of each control volume, that is ∆ = (∆ , ∆ , ∆ ) such as if ≥ ∆ , the method will switch from RANS to LES mode. However, this means that the method is highly grid dependent. In this study, SST-k-ω turbulence model was used as a RANS model. The transient scheme was 2 nd order backward Euler. The discretization scheme for the advection term was high resolution scheme, in which the CDS blending was used between the high resolution and the central difference schemes. Moreover, the high resolution advection and transient schemes were implemented on the turbulence equations. The smallest time step size was chosen such as the root mean squared Courant number was one in the flow field. Another computational parameters were similar to that of the URANS simulation.

SBES simulation
The stress-blended eddy simulation (SBES) is also a hybrid RANS-LES simulation method. However, the switch from RANS to LES is done by using a blending function for shielding a RANS region. The blending is expressed as where SBES between 0 and 1 is the shielding function, is the RANS portion of the predicted stress tensor and is the LES portion of the stress tensor. It follows that in a free shear area, the shielding function approach zero and the method switch from RANS to LES. In this case also, the SST-k-ω turbulence model was used as a RANS method. Furthermore, the other parameters were similar than in the DES simulation.

Test cases
The test cases are summarized in Table 2.

LES-RANS
In the airflow characteristics, the fluctuation energy ratio is expressed as that will show the local portion of fluctuation energy compared to kinetic energy in a dataset record. Each instantaneous air speed can be divided to the fluctuation part and to the mean part. The level of fluctuation can be considered e.g. by calculating the standard deviation (STD) and other descriptive statistics. Fig. 4 shows the numerical solution of URANS (SST-kω), SBES, DES and experiments in the middle of the room. Generally, SBES (Fig. 4b) and DES (Fig. 4c) provided more airflow structures than URANS (Fig. 4a) that can be obtained as a diversity in the airflow field. In the experiments, the downstream airflow inclined slightly to the other side of the room. This asymmetry was also found with CFD. However, the flow direction was case dependent. This is mainly because a large-scale variation of downstream airflow is three-dimensional and the given cross-section represents only a two-dimensional flow field.  5). This is because the average downward flow inclined differently in the experiments than in the predictions (Fig.  4). In the experiments, the average of local mean air velocities was 0.33±0.12 m/s (±standard deviation) at the investigated locations. The corresponding velocity was 0.34±0.11 m/s with SBES, 0.33±0.11 m/s with DES and 0.31±0.11 m/s with SST method. Hence, the average air velocity was at the same level with those methods. However, the range of local mean air velocities was 0.14-0.67 m/s in the experiments and 0.18-0.53 m/s, 0.17-0.50 m/s and 0.15-0.52 m/s with the SBES, DES and SST methods, respectively. Thus, the SST method generated a slightly larger range than SBES and DES that, in turn, was 29% smaller than in the experiments. Furthermore, the lowest predicted maximum was 25% smaller in the DES method than in the experiments. In the middle zone at x=±0.17 m, the mean fluctuation energy was higher in the predictions than in the experiments (Fig. 6). The highest vertical deviation was obtained with SBES method and more similar deviation with DES and SST methods. In addition, lower sampling rate affects the results of the experiments. In this sense, the SBES provided more changes in airflow fluctuation than those other studied methods.  Fig. 7 shows that the mean energy of air velocity (̅̅̅̅) was at a similar level in the middle zone. However, the experiments had higher energy at x=-0.17 m near the ceiling than those other methods (Fig. 7a). This is most probably because the average flow field inclined differently in the experiments than that of the predictions (Fig. 4d). In most locations, the fluctuation was higher with SBES and DES methods than with SST ( Fig. 8 and Fig. 9a). This is reasonable, because URANS do not resolve turbulence scales but merely model the turbulence by increasing the viscosity. In addition, SBES method generated more changes in the fluctuation energy and in a wider range of frequencies than those other studied methods (Fig. 8).

Results
However, this depends highly on the location. This may imply larger sensitivity of RANS-LES methods on producing turbulent motions in the flow. Furthermore, the filtered functions show that the predicted fluctuation was at a reasonable level compared with experiments (Fig. 9b).  The average of local mean air velocity is summarized in Table 3. The results indicate that the CFD simulation can predict reasonably well the averaged air velocity conditions. Furthermore, the average velocity was larger in the middle of the room than in the other regions.

Discussion
More reliable indoor airflow simulations can be made by using an unsteady CFD simulation method than by a steady CFD simulation methods, because the turbulent flow is time-dependent and the unsteady simulation allows airflows to fluctuate. Generally, the results show that numerical modeling method has a considerable effect on a predicted airflow field. However, time-averaged airflow patterns can be rather well predicted, but local differences naturally exist. This means that probable high velocity areas can be coarsely assessed but their accurate locations are difficult to indicate. The modeling can be evolved by using the measured data as the boundary conditions. Currently, the air velocity conditions are not used while controlling HVAC systems in indoor environment. However, local thermal discomfort and thermal comfort models require information on air velocity levels in indoor air. Therefore, the measured and modeled data are suggested to be combined in the future while controlling indoor environments.
In this study, the averaged airflow characteristics were mainly well predicted. The SBES and DES methods provided higher temporal fluctuation than SST method. This is most probably due to RANS-LES combination that generates more reliable airflow structures than the eddy viscosity models. Furthermore, the sensitivity in producing turbulence structures was better with SBES method than with DES or SST methods, particularly near the colliding region where high velocity and gradients exist.
More locations and time-step sizes have to be considered while concluding the comprehensive results of the study. However, this first part shed light on main differences between the studied CFD simulation methods and experiments. Most probably, the CFD simulations can be further improved by adding artificial turbulence to the inlet boundaries. In addition, computing resources are growing all the time in the future. Therefore, further studies are necessary for investigating reliability of advanced modeling methods. The study emphasizes the importance of combining the measurements and numerical modeling of indoor airflows in controlling modern HVAC technology.

Conclusions
The results show that correctly implemented numerical modeling may predict reasonably well the indoor airflow characteristics. The sensitivity on generating airflow fluctuation was greatest with SBES method. The downward flow near the colliding area varied depending on the method. The modeled airflow characteristics differed statistically from the experiments. A recommendation is to carefully test and choose a modeling method for indoor airflows.