Evaluation of mean velocity and mean speed for test ventilated room from RANS and LES CFD modeling

The paper presents and discusses data for the ventilation airflow in an isothermal room corresponding to the Nielsen et al. (1978) test computed with Large Eddy Simulation (LES) and Reynolds-Averaged Navier-Stokes (RANS) approaches. As LES computations provide directly both the speed and velocity components data, the difference between the mean speed and mean velocity values is computed and discussed. For the RANS computations that give the mean velocity data only, application of the velocity-to-speed conversion procedure based on the turbulence kinetic energy field provided by a turbulence model resulted in accurate mean speed evaluation.


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 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]). 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 three-dimensional (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 relatively moderate computational resources required. 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 [4]. 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 [5]. Consequently, RANS solutions demand comprehensive validation, and to validate CFD results, well-documented benchmark experimental data obtained in simplified room configurations are required.
Another principal problem with RANS modeling in HVAC design is that even properly validated RANS results may still be incorrect due to improper interpretation. The reason for that is due to RANS computations provide information on mean components of velocity, <Vx>, <Vy> and <Vz>, that define the mean velocity magnitude, or simply "velocity", Vm  (<Vx> 2 +<Vy> 2 +<Vz> 2 ) 0.5 , here and after <…> denotes time averaging. To evaluate thermal comfort indices, such as Predicted Mean Vote (PMV) and Draught Rating (DR), the mean speed distributions, Va  <(Vx 2 +Vy 2 +Vz 2 ) 0.5 >, must be used, because the air speed is related to the cooling effect of the airflow on the skin, but not the mean velocity magnitude [6,7]. It is known that in case of large velocity fluctuations, especially when the flow direction changes from time to time, that is typical for ventilation problems, the Va values could differ significantly from the Vm values [6].
In HVAC design, to use RANS results that provide information on Vm only, a proper conversion of the Vm data into Va fields is necessary. Empirical correlations derived from ultrasonic anemometer measurements were suggested for mean speed evaluation in [6]. Correlations for RANS velocity processing based on the Laser-Doppler Anemometry (LDA) measurements were proposed later in [7].
Opposite to RANS, the Large Eddy Simulation (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. The overview and outlook of LES could be found in [8]. LES, being an eddy-resolving technique, requires for large computational resources. However, it gives more accurate prediction of turbulent flows and with the parallel computations development becomes more and more popular not only in fundamental studies, but in industrial applications as well.
LES application gives information on instantaneous components of velocity, Vx, Vy and Vz. (sure, the term "instantaneous" covers the resolved frequencies only). An appropriate averaging procedure gives the mean components of velocity (this quantity is available in RANS solutions). Based on the same instantaneous velocity components, another averaging procedure results in the distribution of mean speed, Va, required for thermal comfort indices evaluation. Therefore, LES computations provide both Vm and Va fields that could be used for Vm data processing procedures development and testing. In [9,10] a theoretically developed procedure of Vm processing has been proposed and examined using both the LES and the cabin ventilation qualification test data for the International Space Station pressurized module Columbus.
Recent contributions [11,12] presented wallmodeled LES (WMLES) data on mixing ventilation in a test isothermal room with a sidewall jet for which welldocumented measurement data by Hurnik et al. are available [13,14]. The computational data demonstrated the level of the difference between Vm and Va values both in the jet zone with relatively high velocity values [11] and in the low-velocity occupied zone [12]. As well, three Vm data correction correlations available in the literature [6,7,9] were tested in [12]. For all three correlations, the Va fields processed from the Vm distributions were close to each other and to the mean speed distributions extracted directly from the LES data almost everywhere in the airflow domain, except small particular regions (e.g., the mixing layers).
The present paper considers the most popular isothermal benchmark test by Nielsen et al. [15]. In this test, mean velocity components and fluctuations were measured with the LDA technique, with the uncertainty less than 0.5%; the measurement data are available in [15] and in full on the website http://www.cfdbenchmarks.com/. The test 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 [16] and Voight [17] provide the most complete sets of RANS data under conditions of the test [15]. These theses describe the influence of various numerical parameters (e.g., mesh sensitivity) and the effect of a turbulence model choice on the results obtained. Several contributions present LES data for the Nielsen test, i.e. papers by Davidson et al. [18], Emmerich et al. [19] and Zhang et al. [20]. Though satisfactory agreement with experimental data was reported in the papers cited, relatively coarse computational meshes of 0.5-2 million cells were used there, and a revision of early LES data is anticipated.
Wall-modeled LES (WMLES) of airflow for the Nielsen test conditions has been performed recently with finer meshes [21,22]. In [21] a simplified formulation with the periodicity boundary conditions in the transverse direction was used. It was proved that widthto-height ratio W/H = 1 is enough to provide domainindependent periodic solution. The mesh sensitivity study demonstrated that for the Nielsen test conditions it is enough to use 250 uniform cells per room height to get accurate solution Results of computations for the real test geometry with the transversal side walls were presented in [22]. Comparison of WMLES results with the experimental data on mean and rms velocities resulted in conclusion that the mesh used was sufficient to predict accurately the jet spread, though there was slight overprediction of the jet mixing. At the same time, there was some difference between the computed and measured data in the low-velocity occupied zone, especially for the essential 3D case with narrow inlet slot.
The current contribution presents results of WMLES for the Nielsen test conditions discussed partially in [22] in comparison with the RANS solutions. The focus of the paper is on the velocity-to-speed conversion procedure: the literature correlations are applied to the Vm fields extracted from the WMLES and RANS data. Figure 1 shows the geometry model adopted for the computations of the indoor airflow that corresponds to the test facility [15]. The measurement data [15] are available at eight dashed lines shown in Figure 1: vertical lines A-A located at x/H = 1.0 and B-B at x/H = 2.0 and horizontal lines C-C located at y/H = 0.972 (starting at the midheight of the inlet slot) and D-D located at y/H = 0.028 (that is equal to hin/2). The subscript «1» in the line notation corresponds to the section z/H = 0.5, while subscript «2» -to z/H = 0.1.

Turbulence modeling
Turbulent airflow in the room was computed using different approaches: the first one based on the solving Reynolds-Averaged Navier-Stokes equations (steadystate RANS with standard k-ε turbulence model and unsteady RANS -i.e. URANS -with RNG k-ε model), and the vortex-resolving LES approach that solves the filtered Navier-Stokes equations resolving large scales of motion, while small scales are modeled with SGS model. For the incompressible fluid, the equations for LES calculations are written as follows: where V  is the velocity vector with components (Vx, Vy, Vz), S is the strain rate tensor for the resolved motion, and  SGS is the SGS stress term arising from the spatial filtering procedure. The filtering operation [6] for a variable f determines the filtered (resolved) and smallscale (non-resolved) components f~ and f as follows: where G (xx', ) is a filter function that determines the size and structure of the small scales, x is a coordinate of the point under consideration, and ∆ is the filter width.
To determine the SGS stress term, the generalized Boussinesq hypothesis is used: where SGS is the SGS viscosity, which must be determined by a SGS model. The algebraic Wall-Modeled LES S-Omega SGS model available in ANSYS Fluent was applied in the current study. The model realization in the code is based on [23], and the SGS viscosity is calculated with the use of a hybrid length scale and the wall-damping function: where S and Ω are the strain rate and vorticity magnitude, CS = 0.2 is the Smagorinsky constant, к = 0.41 is the von Karman constant, dw is the distance to the nearest wall, y + is the normal to the wall inner scaling. The grid scale is defined as follows: where Δmax and Δwn are the maximum local grid spacing and the grid step in the wall-normal direction, and Cw = 0.15 is the empirical constant. To specify the time-evolving inflow conditions with superimposed perturbations on inlet mean velocity profiles, the vortex method available in ANSYS Fluent was used. For RANS computations, the ratio of the turbulent to molecular viscosity turb/ = 2.3 was set at the inlet, while the turbulent intensity was equal to 4% (that corresponds to the experimental inlet turbulence characteristics reported in [24]).

Computational grids
For the LES computations, the uniform mesh was used. The mesh was created with the ANSYS ICEM CFD 16.2 mesh generator and consisted of 751252250 (about 48 million) cubic cells with the cell length  = 12 mm. (Note that here and after, in the results description, abbreviations WMLES and LES are used as synonyms.) Detailed mesh-sensitivity analysis was performed earlier in [21]. The maximum y + -values (dimensionless wall distance) did not exceed 0.5 in the occupied zone and were about 20 in the jet zone near the inlet slot. The ratio of the SGS to molecular viscosity is less than 5 everywhere in the computational domain (Figure 2a).
3D mesh for RANS calculations consisted of 300×140×100 (4.2 million in total) hexagonal mesh elements clustered to the walls and to the inlet and outlet sections. It was controlled that the y + values at the walls are below unity almost everywhere, so that the enhanced wall treatment approach was activated. The maximum turb/ value is more than 700 if the standard k-ε turbulence model is used (Figure 2b). Though in the recirculation zone turb/ values are also rather high, up to 400, if the k-ε RNG turbulence model is used (see an instantaneous field given in Figure 2c), it was not possible to get a converged steady-state solution for this turbulence model with the relatively fine mesh used.

Solver settings
Numerical solution was obtained with the CFD package ANSYS Fluent 16.2 based on the finite volume method with the cell-centered variable arrangement. For the LES computations, the non-iterative time advancement scheme (NITA) based on the fractional step method was used. The spatial discretization was performed with the central-differencing scheme for convective terms; the second order pressure interpolation was used. The second order implicit time integration was used. The value of a time step, Δt, for the LES calculations is Δt = 0.01 s; it was chosen to provide the Courant number in the computational domain less than 1. Note that an auxiliary time-step sensitivity study was performed with Δt = 0.006 s, and a decrease in the time step value did not lead to any changes in the time-averaged airflow characteristics. To accumulate representative statistics, the sample of about 150,000 time steps (1500 s) was computed. This period for averaging was assumed to be sufficient to get statistically steady data.
To illustrate local fluctuations, Figure 3 shows velocity magnitude evolution in time at three monitoring points placed at the middle section. The first point, P1, with the coordinates of (3 m, 2.916 m, 1.5 m) is located directly in the jet region. As was expected, for this point the most intensive fluctuations were detected. The second point, P2 (6 m, 1.5 m, 1.5 m), is located in the main recirculation zone, and the corresponding plot demonstrates lower amplitude and frequency of pulsations. Point P3 (8 m, 2.916   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 19 nodes (for WMLES calculations) and 7 nodes (for RANS/URANS calculations) 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, 512 cores were used for WMLES computations and 84 cores were used for RANS/URANS computations.

Mean velocity and mean speed discussion
Instantaneous velocity fields shown in Figure 1 and mean velocity distribution given in Figure 5 illustrate the global airflow pattern in the room. The room airflow pattern could be separated into two zones: the relatively high-velocity zone of the near-wall jet under the ceiling with the mixing layer region and the low-velocity occupied zone with the clockwise recirculation flow. Note that besides the room airflow, Figure 5 demonstrates the flow peculiarities in the outlet channel, where a pronounced separation zone is visible. The length of the outlet channel chosen ensured that the outlet section in the WMLES computations is downstream the attachment point.  [15] profiles of mean longitudinal velocity and its pulsations at the mid-section and at the side-section of the room. It is visible that the jet spreading is reproduced quite accurate (the only slight deviation is in the <Vx> values at large x). At the same time, there is a pronounced difference between the computed velocity profiles and the experimental data in the recirculation zone ( Figures  6c, 7c), especially for the RANS results with the standard k-ε model. Another visible difference between the experimental and computational data is in the rms plots. The pulsations in the low-velocity recirculation zone are slightly under-predicted (see the rms velocity profiles at lines D1-D1 and D2-D2, as well as in the lower half of lines B1-B1 and B2-B2). Note that the rms velocity values in the recirculation zone (especially near the floor) are under-predicted in all previously published LES studies, see, e.g., results and discussion in the recent paper [25].
Time-averaging of LES-data provides directly both the mean speed and mean velocity distributions. For the middle cross-section of the room, these quantities are illustrated in Figures 5 and 8a respectively. The mean speed is larger than the mean velocity everywhere in the domain. The difference between Vm and Va (plotted in Figure 8b for the same mid-section) is less than 10% in the jet core and in the region of the recirculation flow formed just after interaction with the opposite wall where the velocity values are high, and the flow has a predominant direction. The difference between Vm and Va is about 50% in the mixing layer and almost everywhere in the recirculation zone where velocities are low; it is up to 80% and even higher in the regions of the lowest velocity values (near the left side wall -below the inlet, and in the center of the recirculation zone).     Figure 9 shows the difference between the LESobtained Va field and the Va distribution computed using the Vm field (averaged from the same LES data) with the correction procedure proposed by Smirnov et al. [9]: where <k> = (<Vx′ 2 >+<Vy′ 2 >+<Vz′ 2 >)/2 is the resolved turbulence kinetic energy.
The distribution given in Figure 9 demonstrates the reasonability of Va estimation if Vm and <k> data are available (note that both the quantities are usually available in the numerical solution when RANS modeling is applied). Relatively large difference between the Va data and the result of the correction procedure application is visible in the mixing layer only: locally it achieves 5% and even more, up to 10%. For the low-velocity recirculation (occupied) zone, the correction procedure resulted in accurate Va estimation: the difference does not exceed 5%.

|V a -V a
Smirnov | / V a  100% Fig. 9. Mean speed evaluation with correction procedure; data are given at mid-section, z/H = 0.5 Figure 10 shows the Vm fields obtained in RANS/URANS calculations using two turbulence models. A comparison with the LES-obtained field given in Figure 5 allows to conclude that in general both RANS and URANS solutions reproduce the global flow structure (it is in accordance with the velocity profiles discussed above, Figures 6 and 7). However, steady-state computations with the standard k-ε model do not reproduce some peculiarities of the room recirculation zone (compare the low-left sub-region of the velocity distributions in Figure 5 and Figure 10a). Therefore, the velocity-to-speed conversion procedure application is presented for the URANS data only. Figure 11 compares the profiles of Vm and Va for the WMLES and URANS approaches. In both the cases Vmfields are directly extracted from the solutions, while Vafields are evaluated using the correction procedure according to equation (6) [9]. Note that for WMLES and URANS data processing different k distributions were used (resolved k for WMLES and modeled k for URANS). It is visible that at line B1-B1 the URANS Vm profile almost coincide with the WMLES data; k fields in this region are almost similar, and the correction procedure results in the same Va profiles. On the contrary, at line A1-A1 WMLES-and URANS-obtained Vm profiles differ much. However, due to reasonable difference in k fields, even at this section the conversion procedure applied to the URANS data results in the Va profile that is close to the result of the WMLES data processing.

Conclusions
Numerical simulation of indoor airflow in the isothermal ventilated room with the under-ceiling jet at conditions of the benchmark test experiment at Re of 5233 have been performed using the ANSYS Fluent 16.2 CFD package. Comparison of the LES and RANS computational results with the LDA test measurement data on the mean values of the longitudinal velocity components demonstrates good agreement almost everywhere in the jet zone. At the same time, there is visible difference between the computed and measured data in the low-velocity recirculation zone (the "occupied" zone).
Both the mean speed and mean velocity data were extracted from the LES solution. According to the computational data, in the jet core the difference between the mean speed and mean velocity does not exceed 20%, while in the occupied zone on the average it is about 50%, and locally it exceeds 80%. Based on the LES data, evaluation of the mean velocity processing procedure that takes into account the turbulence kinetic energy field was performed. The processed Va fields are very close to the mean speed distributions extracted directly from the LES data. The difference does not exceed 5% everywhere except the region of the mixing layer where it locally achieves the level of 10%. It was shown also that application of the correction procedure to the RANS data on Vm and k provides adequate Va fields for the whole room volume.