A Numerical and Experimental Investigation of a Confluent Jets Ventilation Supply Device in a Conference Room

Abstract

In this study, confluent jets ventilation (CJV) supply devices with three different nozzle arrays (1 × 19, 2 × 19, 3 × 19) were investigated both numerically and experimentally at two different airflow and supply air temperature set-ups. The performance of the CJV supply devices was investigated concerning thermal comfort, indoor air quality (IAQ), and heat removal effectiveness in a conference room environment. A comparison between the experimental and numerical results showed that the $\overline{{\vartheta }^2}-f$ model had the best agreement out of the investigated turbulence models. The numerical results showed that the size of the array had a great impact both on near-field development and on the conditions in the occupied zone. A larger array with multiple rows and a lower momentum conserved the inlet temperature and the mean age of the air better than a single-row array with a higher momentum. A larger array with multiple rows had a higher IAQ and a greater heat removal effectiveness in the occupied zone because the larger array conserved the mean age of air better and the buoyancy driven flow was slightly better at removing the heat. Because of the lower inlet velocities, they also had lower velocities at ankle level, which decreased the risk of draft and thermal discomfort.

1. Introduction

Heating, ventilation, and air conditioning (HVAC) systems account for about 50% of the total energy use in buildings and more than 10% of total national energy use [1]. The aim of ventilation is to maintain a comfortable indoor climate and indoor air quality (IAQ). It is therefore important to design ventilation supply devices that utilize energy-efficient ventilation strategies that are appropriate for the constraints of the ventilated space.

The most documented and well-defined ventilation strategy is mixing ventilation [2]. The principle behind mixing ventilation is to dilute the contaminated air in the room by supplying fresh air at higher velocities (often at ceiling level) to mix the air and lower the concentration of contaminants. The mixing conditions generally generate good thermal comfort and a uniform indoor environment. The Air Change Effectiveness (ACE) and Heat Removal Effectiveness (ε_T) for well-mixed conditions in a mixing ventilation system is 100% [3], which means that the occupied zone has a uniform temperature and an even dilution of contaminants.

Displacement ventilation is an air distribution system that uses natural buoyancy from the heat sources in the occupied zone to create stratified conditions. This is ensured by supplying cool air (3–5 °C below ambient air temperature) [4] at the floor level with low inlet velocities (<0.5 m/s) [2]. The cooler air displaces the warmer air and the local contaminants, which follow the natural plumes to ceiling level where they are removed via the exhaust [5]. Displacement ventilation generally has an ACE and ε_T above 100% [2]. This increase in efficiency has been calculated to reduce energy usage by as much as 20–34% in test chambers [6,7], 21–25% in offices [8,9], 41% in classrooms [8], and 33% in retail shops [8] compared to mixing ventilation. However, some studies have shown that displacement ventilation has some problems with thermal comfort [2,4,10,11], as the stratification can be disrupted by movement or obstacles in the occupied zone [4,12] and the penetration of the air supply is limited in cases with long distances and high heat loads (>35 W/m²) in the occupied zone [2,13].

There are many so-called hybrid ventilation strategies that combine different aspects of both mixing and displacement ventilation, such as impinging jet ventilation [2,4,14,15,16], under floor air distribution [4,17,18], wall attached ventilation [4,19], and confluent jet ventilation (CJV) [2,4,20]. The aim of the hybrid systems is generally to overcome some of the limitations of displacement ventilation by increasing the momentum of the supplied air while still maintaining semi-stratified conditions [4].

Confluent jets (CJ) can be defined as multiple interacting jets issuing from different nozzles [3]. Numerous studies have been completed on the near-field development of confluent jets both experimentally [21,22,23,24,25] and numerically [26,27,28,29]. Generally, the near-field development of CJ is divided into initial, converging merging, and combined regions. After a certain distance downstream from the nozzles, the combined jets behave as a single jet, where the individual jets can no longer be identified; this region is called the confluent core zone. In an array of free-moving jets where the number of columns is the same as the number of rows, the length of the regions downstream can be estimated based on the jet diameter (d) [23,24,26,28]. In the initial region (0d to 2d) the jets show no sign of interaction with each other or any decrease in velocity. In the converging region (2d to 4d) the jets bend towards the central line of the array and the jets at the edges bend faster than those close to the center. In the merging region (4d to ~9d–14d) the jets start to merge, which happens faster for the edge jets and later for the central jets. In the combined region (~9d–14d to ~20d) the minimum velocity of the combined flow field is no less than 90% of the maximum velocity. In the confluent core zone (~20d to ~11 × d × n, where n is the number of rows) the jets have combined to a single jet with uniform streamwise velocity and low turbulence in the central line of the array. The two parameters that have the highest influence on the near-field development are the jet diameter (d) and the spacing between the jets (S). The combined factor (S/d) has the highest influence on the jet development. An increase in (S/d) increases the length of the converging and merging regions. A higher (S/d) also prolongs the confluent core zone but lowers its velocity and turbulence. The factor (S/d) does not however affect the entrainment rate of the combined jet. CJ on a curved surface behaves similarly, but will produce a slightly more oblong jet instead of a uniform one [27,30].

When CJ are implemented in the ventilation supply device, the array of jets is placed parallel to a wall directed towards the floor, called a wall confluent jet (WCJ). The jets attach to the wall due to the Coanda effect and become a wall jet. This occurs in three regions, which can be described based on the jet diameter (d), the distance between the jets and the wall (x), and the distance from the nozzles (r) [21]. The process begins with the free jet region (0 < √(r/d) < 5.9) where the jets bend towards the wall and decelerate due to the adverse pressure gradient. How fast the airflow decelerates in this region depends on the distance to the wall. Airflow close to the wall (√(x/d) < 1.6) decelerates less than airflow in the centerline of the array due to less entrainment. In the Coanda effect region (5.9 < √(r/d) < 11) the combined jets behave as a wall jet with similar velocities both close to the wall and in the centerline. At the start of the wall jet region (5.9 < √(r/d) < 11), the airflow behaves as a wall jet with similar velocities and decays near the wall and in the centerline. However, at the end of the wall jet region (15.5 < √(r/d)) the air in the centerline is more influenced by outside forces (impinging pressure from floors, etc.) [21]. Studies have shown that WCJ conserves the momentum better than other types of jets (free confluent jets, free plane jet, and plane wall jet) [21,31].

The application of CJ for hybrid ventilation strategies has been studied since the early 2000s [20,32]. Numerous experimental and numerical studies have been performed on CJV in classrooms [10,32,33,34,35,36,37], offices [17,25,31,38,39,40,41,42,43,44] and assembly halls [45]. The results show a span of ε_T from 95% to 131% and an ACE between 99 and 124% depending on the type of room and the configuration of the supply device.

Since the design parameters for CJV (nozzle spacing, nozzle diameter, jet exit velocity, and jet direction) greatly affect the airflow development and conditions in the occupied zone, it is vital to investigate how these parameters interact under different operational conditions (airflow rate and supply temperature) and how these interactions affect the CJV performance. This study is part of a series of studies [10,33,46] investigating a novel CJ supply device that uses a lower momentum than most other CJ devices and is mounted in the ceiling rather than a wall. The aim of this series of studies is to optimize the configuration of this kind of supply device for energy efficiency, IAQ, and thermal comfort.

The key objectives of this work are:

−

To validate a numerical turbulence model for simulating low momentum confluent jets in a conference room environment;

−

To investigate how the number of nozzle rows affects the near-field of the confluent jets’ (0–50d) temperature and flow development as well as thermal comfort, indoor air quality (IAQ), and heat removal effectiveness.

2. Experimental Set-Up

2.1. The Studied Supply Devices

In the present study, the supply device was a circular duct (0.25 m diameter and 2 m long) mounted 0.05 m below the ceiling level which had 57 nozzles on each side (d = 0.028 m) placed in three rows with 19 nozzles each. The configuration of each side can be seen in Figure 1, the spacing between columns was 3 jet diameters and the spacing between the top and bottom rows was 1 jet diameter, see [33] for more details, shown in Figure 1 and Figure 2. Supply device nozzles can be fitted with caps that block the nozzles.

2.2. Test Facilities

The study was performed in a climate chamber (7.2 × 4.2 × 2.7 m³) at the University of Gävle in Sweden, representing a typical conference room environment. The climate chamber is versatile and has previously been used to investigate different indoor environments, e.g., open-plan offices [17,41], classrooms [33,34,35,47], and small offices [15,48].

2.3. Equipment

Thermocouples and constant temperature anemometers (CTA) were used to measure temperatures. In total, 38 thermocouples (type K with an accuracy of ±0.2° C) were used. They were used to measure the wall, ceiling, and floor temperature, and to measure the vertical temperature gradient (VTG) at two different locations in the occupied zone (front and back), see Figure 3. Additionally, two thermocouples measured the supply and outlet temperature.

CTA probes were used to measure the air speed, root mean square (RMS) of the velocity and temperature. Eight probes were used to measure two velocity profiles in the ceiling in front of the middle of the supply device, one velocity profile was located 0.6 m from the center of the supply device (called velocity profile 0.6 m) and one velocity profile was located 1.2 m from the center of the supply device (called velocity profile 1.2 m). Twenty-two probes were placed in eleven pairs with a distance of 0.45 m at eleven different heights (0.1, 0.3, 0.5, 0.7. 0.9, 1.1, 1.3, 1.5, 1.7, 1.9 and 2.1 m) and traversed over a distance of 2.5 m (stopping at 11 locations) to measure the wall, manikin, and symmetry zone on both sides of the room, see Figure 3. At each location, the traverser stopped, and 2 × 15 min measurements were carried out. After the measurements, the traversing system was moved to measure a new series of 2 × 11 × 11 measurement points for a total of three series. This gave a total number of 734 (11 × 11 × 2 × 3 + 4 × 2) measurement points for the CTA probes. The CTA probes were calibrated for the velocity range 0.05 and 1.00 m/s and had an air speed accuracy of ±0.05 m/s anda temperature sensor accuracy of ±0.2 °C. The CTA probes had a response time of 0.2 s to 90% of a step change and the sampling interval was 60 s. The air speed was averaged over the 2 × 15 min measurements for each point.

To calculate the nominal time constant (τn) and the local mean age of air (τp), the decay method was used by measuring the local mean age of air at the locations A, B, C, D, E, and F (outlet) with tracer gas (sulfur hexafluoride (SF6)). The human heat load was simulated using six manikins consisting of a 0.32 m diameter steel tube covered with fabric, with a sensible heat load of 100 W each, see [49]. The measurement equipment added 50 W of extra internal heat load and no lighting was used during the measurement. Figure 3 shows the positions of the measurement points.

2.4. Case Set-Up

Six experimental cases were investigated, see

Table 1. Case set-up.

	Nozzles	U0 (m/s)	Heat Load (W/m2)	ACH (-)
Case 1—17.5 °C-50 L/s-1R	38	2.2	22	2.3
Case 2—17.5 °C-50 L/s-2R	76	1.1	22	2.3
Case 3—17.5 °C-50 L/s-3R	114	0.7	22	2.3
Case 4—16.5 °C-70 L/s-1R	38	3.0	22	3.2
Case 5—16.5 °C-70 L/s-2R	76	1.5	22	3.2
Case 6—16.5 °C-70 L/s-3R	114	1.0	22	3.2

. The cases were named after supply temperature, airflow rate, and number of rows. Three cases had T_S = 17.5 °C and an airflow rate of 50 L/s. Each of these cases had one, two, or three rows of nozzles open on each side of the supply device, while three cases had T_S = 16.5 °C and 70 L/s. The airflow rate and number of nozzle rows affected the jet inlet velocity (U₀) for each supply device and case, see Table 1. The two airflow rates were chosen together with the number of nozzles to give a wide span of inlet velocities. The airflow rate 50 L/s conformed to Category (III) in the EN 16798–1:2019 standard (low polluting building, LPB-2) and the airflow rate 70 L/s conformed to Category (II) in the EN 16798–1:2019 standard (low polluting building, LPB-2) [50]. The Nozzle configuration for the 3 rows can be seen in Figure 1, for the 2-row cases the top row was blocked, for the 1-row cases the two top rows were blocked.

2.5. Measurement Procedure and Analysis

Each case was set up with airflow rate, supply temperature, and number of open rows. The measurement was running for at least 16 h to ensure quasi-steady state conditions. The temperature data were collected and analyzed to make sure that steady-state conditions had been reached. The CTA probes were then traversed in three series to measure the air speed and temperature. Afterwards, the CTA measurements were taken.

Trace gas measurements were taken by distributing the gas throughout the room and mixing room air with movable fans (for three minutes) to ensure a homogeneous concentration of about 400 ppm. Tracer gas measurements were left running until the concentration levels became negligible. Previous laboratory studies using the same climate chamber and measurement procedure estimated the tracer gas measurements to have an accuracy of ±7% ([15,41,48]).

The local mean age of air (ACE_p) was calculated according to:

$${\mathrm{ACE}}_{\mathrm{p}}=\raisebox{1ex}{${\tau }_{\mathrm{n}}$}\!\left/ \!\raisebox{-1ex}{${\tau }_{\mathrm{p}}$}\right.$$

(1)

where τ_n is the nominal time constant and τ_p is the local mean age of air. The local heat removal effectiveness (ε_T) was calculated according to:

$${\epsilon }_{\mathrm{T}}=\frac{T_{\mathrm{E}}-T_{\mathrm{S}}}{T_{\mathrm{P}}-T_{\mathrm{S}}}$$

(2)

where T_E, T_S, and T_P are the temperatures for the exhaust, supply, and the specific measurement point, respectively. The Archimedes number (AR) was calculated according to [5]:

$$\mathrm{AR} =\frac{\beta g\mathsf{\Delta }T\sqrt{A_0}}{U_0^2}$$

(3)

where ΔT is the difference between the inlet temperature and the wall temperature, β is the volume thermal expansion coefficient, g is the gravitational constant, A₀ is the area of the inlet, and U₀ is inlet velocity [5].

The draft rating was calculated using the equation [51]:

$$\mathrm{DR} =\left(34-T_a\right)\times {\left(U_a-0.05\right)}^{0.62}\times \left(0.37U_a\times TI+3.14\right)$$

(4)

where T_a is air temperature, U_a is the air speed, and TI is turbulence intensity.

3. Computational Set-Up and Numerical Schemes

3.1. Geometrical Set-Up and Boundary Conditions

The geometrical set-up had the same dimensions as the physical set-up, but was divided by a symmetry plane in the middle of the 4.2 m side. The model therefore only modeled one half of the room (the right side), see Figure 3, in order to reduce the computation time. The manikins had the same dimensions as the physical ones and were simulated as internal heat loads with 100 W each. The velocity components (u, v and w) for each individual nozzle were calculated from a numerical model (Andersson et al., 2018) of the supply device, where k and ω were used to determine the k, ε and $\overline{{\vartheta }^2}$ at the inlet. The experimental measurements for the inlet/wall/ceiling/floor temperatures were used to set the corresponding temperature in the numerical model (see Section 4.1). The pressure outlet was chosen for the outlet boundary condition and the outlet was extended by 8d in the numerical model by axial extrusion in order to ensure convergence. All the surfaces were assumed to be grey with an emissivity of 0.95, except for the highly polished surface of the supply device which was set to 0.7.

3.2. Governing Equations

The model was assumed to be steady-state and three-dimensional. The buoyancy effect was included in the momentum equation, where the density was modeled by the incompressible ideal gas law. The radiation heat was accounted for by the discrete ordinates (DO) model. Based on these assumptions, the Reynolds-Averaged Naiver–Stokes (RANS) equations were given by:

$$\frac{\partial (U_i)}{\partial x_i}=0$$

(5)

$$\frac{\partial \left(U_j{U}_i\right)}{\partial x_j}=-\frac{1}{\rho }\frac{\partial P}{\partial x_i}+\upsilon \frac{\partial }{\partial x_j}\left(\frac{\partial U_i}{\partial x_j}+\frac{\partial U_j}{\partial x_i}\right)+\frac{\partial }{\partial x_j}\left(-\overline{u_i^{\prime }{u}_j^{\prime }}\right)+\left(1-\frac{{\rho }_o}{\rho }\right)g_i$$

(6)

$$\frac{\partial \left(U_{j}T\right)}{\partial x_j}=\frac{\partial }{\partial x_j}\left(\alpha \frac{\partial T}{\partial x_j}\right)+\frac{\partial }{\partial x_j}\left(-\overline{u_i^{\prime }{T}^{\prime }}\right)$$

(7)

where $\overline{u_i^{\prime }{u}_j^{\prime }}$ and $\overline{u_i^{\prime }{\theta }^{\prime }}$ are extra terms and unknown, called Reynolds stresses and turbulent heat fluxes, respectively. By using the Boussinesq hypothesis, ($\overline{u_i^{\prime }{u}_j^{\prime }}$) and ($\overline{u_i^{\prime }\theta }$) are defined by:

$$\overline{u_i^{\prime }{u}_j^{\prime }}=-{\upsilon }_t\left(\frac{\partial U_i}{\partial x_j}+\frac{\partial U_j}{\partial x_i}\right)+\frac{2}{3}{\delta }_{ij}k$$

(8)

where ${\upsilon }_t$ is the kinematic turbulence viscosity with dimensions m²/s. ${\upsilon }_t$ can be expressed as the product of a turbulence velocity scale and length scale. k is turbulent kinetic energy defined as $k=\left(\overline{u_i^{\prime }{u}_i^{\prime }}\right)/2$ and δ_ij is the Kronecker delta:

$$S_{ij}=0.5\left(\frac{\partial U_i}{\partial x_j}+\frac{\partial U_j}{\partial x_i}\right)$$

(9)

The turbulence heat fluxes are modelled in a similar way by:

$$\overline{u_i^{\prime }{\theta }^{\prime }}=-\frac{{\upsilon }_t}{{\sigma }_t}\frac{\partial \Theta }{\partial x_j}$$

(10)

where σ_t is turbulence Prandtl number and assumed to be constant in this study.

3.3. Turbulence Modeling

The choice of turbulence model needed to be investigated cautiously. The physics of the flow, knowledge about the pros and cons of the turbulence models, the outcome from previous similar studies [43], the CPU and time resources as well as the computational accuracy are examples of the issues one needs to be aware of. In a series of simulations, five different turbulence models, i.e., the k—ε model, the k—ε model RNG, the k—ω model, the SST k—ω model, and the $\overline{{\vartheta }^2}-f$ model were investigated. The $\overline{{\vartheta }^2}-f$ model was found to have the best accuracy and had the smallest number of numerical instabilities.

In the $\overline{{\vartheta }^2}-f$ model, the additional transport equation of $\overline{{\vartheta }^2}$ is solved together with an equation for elliptical relaxation function (f). This model solved the governing equation of k and ε all the way down to the wall surface without using wall functions. The velocity fluctuation normal to the streamlines was used as the velocity scale to calculate the near-wall turbulence eddy viscosity. The turbulence kinetic energy k and its dissipation ε, the wall normal stress $\overline{{\vartheta }^2}$, and the elliptic relaxation function, f are given by:

$$\frac{\partial \left(\rho U_{j}k\right)}{\partial x_j}=\frac{\partial }{\partial x_j}\left(\left(\mu +\frac{{\mu }_t}{{\sigma }_k}\right)\frac{\partial k}{\partial x_j}\right)+P_k-\rho \epsilon$$

(11)

$$\frac{\partial \left(\rho U_j\epsilon \right)}{\partial x_j}=\frac{\partial }{\partial x_j}\left(\left(\mu +\frac{{\mu }_t}{{\sigma }_{\epsilon }}\right)\frac{\partial \epsilon }{\partial x_j}\right)+\frac{C_{\epsilon 1}^{,}{P}_k-C_{\epsilon 2}\rho \epsilon }{T_{ts}}$$

(12)

$$\frac{\partial \left(\rho U_j\overline{v^2}\right)}{\partial x_j}=\frac{\partial }{\partial x_j}\left(\left(\mu +\frac{{\mu }_t}{{\sigma }_k}\right)\frac{\partial \overline{v^2}}{\partial x_j}\right)+\rho kf+6\rho \overline{v^2}\frac{\epsilon }{k}$$

(13)

$$f-L^2\frac{{\partial }^{2}f}{\partial x_i\partial x_i}=\frac{\left(C_1-1\right)}{T_{ts}}\left(\frac{2}{3}-\frac{\overline{v^2}}{k}\right)+C_2\frac{1}{\rho k}{\mu }_t\left(\frac{\partial U_i}{\partial x_j}+\frac{\partial U_j}{\partial x_i}\right)\frac{\partial U_i}{\partial x_j}+\frac{5}{T_{ts}}\frac{\overline{v^2}}{k}$$

(14)

With:

Eddy viscosity:

$${\mu }_t=\rho {\mu }_t\overline{v^2}{T}_{ts}$$

(15)

Turbulence production:

$$P_k={\mu }_t{S}^2 \to \left(S=\sqrt{2S_{ij}{S}_{ij}}\right)$$

(16)

Turbulence time scale:

$$T_{ts}=\min \left[\max \left(\frac{K}{\epsilon },C_t{\left(\frac{v}{\epsilon }\right)}^{\frac{1}{2}}\right),\frac{\alpha }{\sqrt{3}}\frac{k}{\overline{v^2}{C}_{\mu }\sqrt{2S_{ij}{S}_{ij}}}\right]$$

(17)

Turbulence length:

$$L=C_L\max \left[\min \left[\frac{k^{\frac{3}{2}}}{\epsilon },\frac{1}{\sqrt{3}}\frac{k^{\frac{3}{2}}}{\overline{v^2}{C}_{\mu }\sqrt{2S_{ij}{S}_{ij}}}\right], C_{\eta }{\left(\frac{v^3}{\epsilon }\right)}^{1/4}\right]$$

(18)

The constants are:

$$\begin{gathered} C_{\epsilon 1}^{\prime }=C_{\epsilon 1}\left(1+0.045\sqrt{k/\overline{v^2}}\right) , C_{\epsilon 1}=1.4 , C_{\epsilon 2}=1.9 ,C_1=1.4 ,C_2=0.3 , C_{\mu }=0.22 , \\ {C}_L=0.23 , C_{\eta }=70 , {\sigma }_{\epsilon }=1.3m , {\sigma }_k=1.0, \alpha =0.6 , \end{gathered}$$

3.4. Numerical Details

The commercial finite volume solver Fluent 19.5 (Ansys 2019) was used to numerically solve the governing equations. The installed memory was 128 GB and 24 × 2.50 GHz processors (Intel Xeon E5–2680 v3) were used. The SIMPLE algorithm was used to control the pressure–velocity coupling. The gradient was solved with the least squares cell-based method. The pressure term used the PRESTO! scheme and an under-relaxation factor of 0.5. The momentum term used a third order-scheme (3rd MUSCLE) and an under-relaxation factor of 0.5. Both terms for the turbulent kinetic energy and the specific dissipation rate used a third-order scheme (3rd MUSCLE) and both used an under-relaxation factor of 0.8. The parameters for the DO Radiation model were set to 10 energy iterations per radiation iteration, 2 Theta Divisions, 3 Phi Divisions, 1 Theta Pixels, and 3 Phi Pixels. The solution was considered converged if the residual was less than 10⁻³ for the continuity and less than 10⁻⁴ for u, v, w, k, ε, f and $\overline{{\vartheta }^2}$. At least 20,000 iterations were used for each case.3.5. Mesh Strategy

Different mesh densities and grid configurations were tested in order to establish grid independence with regards to spatial convergence. The difference in the velocity and temperature results for different mesh densities was calculated according to:

$$Difference U \% =\frac{{\sum }_{i=1}^{ 365}\sqrt{{\left(\frac{U_{i,Grid A}-U_{i,Grid B}}{U_0}\right)}^2}}{365}$$

(19)

$$Difference T \% =\frac{{\sum }_{i=1}^{ 365}\sqrt{{\left(\frac{T_{i,Grid A}-T_{i,Grid B}}{T_E-T_S}\right)}^2}}{365}$$

(20)

The mesh size was increased by roughly doubling the number of cells between each mesh refinement from 6.5 to 24.5 million in two steps. The difference in results was 6.6% for the velocities and 0.6% for the temperatures in the first step, and 0.5% for the velocities and 0.5% for the temperatures in the second step. Since the second step between 12.2 million cells and 24.5 million cells changed the result by only 0.5%, the mesh with 12.2 million cells was used.

A range of different mesh strategies and densities were used in the different parts of the geometry. The highest mesh density was located in and around the nozzles. An unstructured tetrahedral mesh with cell sizes between 1.0 mm and 15 mm was used in the critical areas around the nozzles to better capture the sharp velocity gradients (see Figure 4 right side). The total number of cells in this part of the mesh was 2.7 million. The ceiling and wall part of the mesh downstream from the supply device where the near-field development of the confluent jet was assumed to occur was meshed using uniform quad strategy with 3.9 million cells with a maximum cell size of 10 × 16 × 40 mm (see Figure 4 left side). An inflation of 20 layers was used at the walls and ceiling to keep y+ less than 1. The part around the three manikins had an unstructured mesh with an inflation layer at the floor and around the manikins and the maximum cell size was 40 mm. Likewise, the outlet had an unstructured mesh with inflation layers and the outlet was extended to 8 diameters of the outlet to prevent backflow. The other parts of the mesh were all structured mesh with a maximum cell size of 40 × 40 × 16 mm with inflation layers at the walls, ceilings, and floors. The mesh was the computationally efficient compromise between grid independence, aspect ratio, skewness, and y+, see

Table 2. Mesh skewness, orthogonal quality, aspect ratio, and y+.

Skewness
Quality	Excellent	Very good	Good	Acceptable	Poor	Unacceptable
Value	0–0.25	0.25–0.50	0.50–0.80	0.80–0.94	0.95–0.97	0.98–1.00
Number of cells	9,763,910	2,091,122	356,092	26,405	1	0
% of cells	79.8%	17.1%	2.9%	0.2%	0.0%	0.0%
Orthogonal Quality
Quality	Excellent	Very good	Good	Acceptable	Poor	Unacceptable
Value	0.95–1.00	0.70–0.95	0.20–69	0.15–0.20	0.001–0.14	0–0.001
Number of cells	6,594,682	3,825,128	1,676,456	15,498	717	0
% of cells	54.4%	31.6%	13.8%	0.1%	0.0%	0.0%
Aspect Ratio
Value	1–10	10–25	25–35	35–45	45–55	55–68.6
Number of cells	10,928,255	1,061,709	119,207	117,007	850	319
% of cells	89.4%	8.7%	1.0%	1.0%	0.0%	0.0%
y+
Case	1	2	3	4	5	6
Avg	1	0.9	0.8	1.1	1	0.9
Max	2.9	2.5	2.1	3.1	2.7	2.5

.

4. Results and Discussion

4.1. Experimental Results

4.1.1. Air Speed and Temperature Measurement

The results from the air-speed measurements show the same pattern for all six cases. The airflow attached itself to the ceiling due to the Coanda effect as shown by the two velocity profiles, see Figure 5. The airflow had the highest velocities (0.25–0.71 m/s) at 2.6 m close to the ceiling and dropped to 0.1–0.13 m/s at the heights of 2.3 and 2.4 m. Due to the entrainment of the surrounding air, the confluent jet stream became wider and the maximum velocity decreased with distance from the inlets. The airflow followed the ceiling until it hit the wall opposite of the supply device and entered the wall zone, where it had higher velocities in the center of the confluent jet. The velocity profile of the confluent jet smoothed out as the jet approached the floor, see Figure 5. A distinct discontinuity in the pattern occurred between height 1.1 m and 1.3 m for all six cases where the velocity profiles dropped by roughly 0.05–0.10 m/s. This difference in velocity also occurred in the manikin zone and the symmetry zone. An explanation for this discontinuity can be found in Section 4.3.2 of this paper. Some of the air followed the floor into the occupied zone, which was indicated by the higher velocities close to the floor at a height of 0.1 m, see Figure 5 (middle). The velocities in the occupied zone (manikin and symmetry zones) were very low (<0.2 m/s) even for the case with the highest inlet velocity (Case 4, see Figure 4). If one compares the velocities in the wall and symmetry zones, one can see that the velocities were slightly higher at the back of the measured distance in the wall zone, whereas in the symmetry zone they were slightly higher in the front, close to the outlet. The velocity profiles on the left and right side of the room were very similar and the patterns were repeated on both sides with only minor variations in amplitude in some points.

As

Table 3. Measured average velocity in different zones.

	Case 1	Case 2	Case 3	Avg.	Case 4	Case 5	Case 6	Avg.
Wall Zone avg.	0.16	0.14	0.14	0.15	0.21	0.19	0.16	0.19
Man. Zone avg.	0.09	0.09	0.10	0.09	0.10	0.09	0.09	0.09
Sym. Zone avg.	0.11	0.11	0.10	0.11	0.12	0.11	0.10	0.11

shows, the average velocity in the wall zone was dependent upon the number of rows as the cases with only one row of nozzles (Cases 1 and 4) had the highest velocities in this zone. However, the average velocities in the manikin and symmetry zones were independent of the inlet velocity as well as the airflow rate and supply temperature.

Table 4. Average temperature measurement for all cases.

	Case 1	Case 2	Case 3	Avg.	Case 4	Case 5	Case 6	Avg.
Inlet(°C)	17.49	17.60	17.42	17.50	16.49	16.43	16.37	16.43
Outlet (°C)	24.70	24.80	24.83	24.78	23.38	22.76	22.78	22.97
Diff-in-out (°C)	7.21	7.20	7.41	7.27	6.89	6.33	6.41	6.54
Wall Zone (°C) avg.	24.46	24.52	24.39	24.46	22.95	22.33	22.41	22.56
Man. Zone (°C) avg.	24.67	24.73	24.63	24.68	23.28	22.67	22.73	22.89
Sym. Zone (°C) avg.	24.60	24.66	24.64	24.63	23.36	22.75	22.54	22.88
VTG front (°C) avg.	24.78	24.87	24.86	24.84	23.49	22.86	22.86	23.07
VTG back (°C) avg.	24.21	24.26	24.28	24.25	22.99	22.35	22.33	22.56
Left wall (°C) avg.	24.50	24.67	24.71	24.63	23.42	22.74	22.75	22.97
Right Wall (°C) avg.	24.62	24.68	24.79	24.70	23.37	22.73	22.74	22.94
Ceiling (°C) avg.	23.82	23.94	23.94	23.90	22.47	22.07	22.05	22.19
Floor (°C) avg.	24.50	24.67	24.64	24.60	23.38	22.73	22.71	22.94
Back wall (°C) avg.	24.37	24.53	24.95	24.61	23.34	22.71	22.69	22.91
Front wall (°C) avg.	25.04	25.24	25.26	25.18	24.00	23.35	23.32	23.55
LHRE	1.02	1.02	1.04	1.03	1.03	1.03	1.04	1.03

shows the average temperatures in the different zones in the room, which shows that the temperature in the wall zone was on average 0.2–0.3 °C colder than in the two other zones because of the higher velocities and the lower infiltration of the surrounding air. The cases with a lower supply temperature and a higher airflow rate (Cases 4–6) had on average a 1.8 °C lower room temperature than the other three cases. The outlet temperature was slightly higher than the temperature in the occupied zone for all six cases, which resulted in an average heat removal effectiveness of 1.02–1.04, which was close to a mixing system. All measured vertical temperature gradients (VTG) were very small (<1.0 °C) and the overall temperature distributions in the occupied zone were uniform with only small variations. There was no considerable difference in the wall temperature between the left and right side for any of the six cases. There was a consistent difference in temperature between the front side and the back side of the room. The temperature of both the front VTG and the front wall was 0.5 °C warmer than the back VTG and back wall for all six cases. This was probably because the heat load (manikins) was slightly closer to the front of the room. The floor was slightly warmer than the ceiling because of the cooling effect from the jets and the heat from the manikins.

Figure 6 shows the temperature measurement for Case 4, which had the largest difference in temperature difference between left and right side. All cases had slightly higher temperatures on the right side of the room. The pattern of the temperature curves was the opposite of the velocity profiles with the lowest temperature in the middle of the horizontal stretch in the wall zone. Higher temperatures were in the front of symmetry zone and lower temperatures were in the back of the symmetry zone, compare Figure 4 and Figure 5. The amplitude of the pattern was somewhat higher on the left side of the room in Case 4, see Figure 5. It is worth noting that the difference in velocities between the heights of 1.1 m and 1.3 m did not significantly affect the VTG in any of the zones. The velocities increased above 1.1 m, but the VTG remained relatively flat above 0.7 m in all three zones on both sides of the room, see Figure 5. See the explanation in Section 4.3.2.

The average absolute difference in air speed between each of the measurement points on the left and the right side of the room can be calculated according to the equation:

$$U-Difference \% =\frac{{\sum }_{j=1}^{363}\sqrt{{\left(\frac{U_{Rigth.j}-U_{Left.j}}{U_0}\right)}^2}}{363}$$

(21)

Likewise, the average absolute difference in air speed between each of the measurement points on the left and the right side of the room can be calculated according to the equation:

$$T-Difference \% =\frac{{\sum }_{j=1}^{363}\sqrt{{\left(\frac{T_{Rigth.j}-T_{Left.j}}{\mathrm{Texhaust}-\mathrm{Tsupply} }\right)}^2}}{363}$$

(22)

Table 5. Difference between the left and right side of the room.

	Case 1	Case 2	Case 3	Case 4	Case 5	Case 6	Avg.
Air Speed
(m/s)	0.02	0.01	0.02	0.02	0.02	0.01	0.02
(%)	0.9%	0.9%	2.9%	0.7%	1.3%	1.0%	1.3%
Temperature
(°C)	0.31	0.11	0.15	0.41	0.10	0.20	0.21
(%)	4%	2%	2%	6%	2%	3%	3%

shows that with regards to velocity, the differences between the left and right side were on average much smaller than the accuracy of the CTA equipment. With regards to the temperature, the difference was smaller than the accuracy for four of the six cases. The two cases with one row of nozzles, Cases 1 and 4, had slightly more differences than the accuracy and a slightly uneven temperature distribution. Although the two cases had a slightly uneven temperature distribution between the two sides of the room, the shape of the temperature and velocity profiles was very similar and the overall statistical analysis shows that the air distribution can be regarded as symmetrical. This means that a numerical model with a symmetry plane along the middle of the supply device should be able to predict the velocity and temperature distribution by only simulating one side of the room using a symmetry plane and therefore reduce the simulation cost by 50%.

4.1.2. Tracer Gas Measurement

Table 6. ACE_p and mean age of air for all six cases.

	ACEP						MEAN AGE OF AIR [MIN]
CASE:	1	2	3	4	5	6	1	2	3	4	5	6
A	1.02	0.99	0.96	0.95	0.95	1.01	25.1	25.4	27.1	20.2	19.2	19.2
B	1.04	1.02	1.02	0.98	1.03	1.03	24.6	24.6	25.6	19.5	17.8	18.8
C	1.01	1.03	1.05	0.96	0.97	1.00	25.3	24.4	24.8	19.8	18.8	19.4
D	1.17	1.12	1.11	1.01	1.09	1.13	21.9	22.5	23.5	19	16.8	17.2
E	1.15	1.10	1.09	1.09	1.08	1.20	22.2	22.9	24.0	17.6	16.9	16.1
OUTLET	1.00	1.00	1.00	1.00	1.00	1.00	25.6	25.2	26.0	19.1	18.3	19.4

shows some variation in the total airflow between the cases. The nominal time constant τ_n (F outlet) was 25.6 min for cases 1–3 and 18.9 min for cases 4–6. This means that mean age of air in the outlet for all the cases was very close to the nominal case (<4%). The ACE_p values in the occupied zone (A–C) were close to mixing conditions (0.95–1.05) and the points in the ceiling close to the supply device were higher at 1.01–1.20 with an average of 1.12. The flow was not fully mixed when the air passed by the measurement points D and E.

4.2. Validation of the Numerical Model

The numerical model was validated via a statistical analysis by comparing each case against the measured results. Since there was a slight difference in the temperature profiles of the two sides of the room, an average (over the two sides) of both the velocity and temperature measurement was used for the validation of the numerical model. In the cases where there was an oscillation in the results (Cases 2 and 3), average values from the CFD results over one cycle (4000 iterations) were used for the statistical analysis (the case with the best statistical match is presented in the figures). The error between the measured and predicted result was calculated according to:

$$T-Difference \% =\frac{{\sum }_{j=1}^{363}\sqrt{{\left(\frac{T_{CTA.j}-T_{CFD.j}}{T_E-T_I}\right)}^2}}{363}$$

(23)

$$U-Difference \% =\frac{{\sum }_{j=1}^{363}\sqrt{{\left(\frac{U_{CTA.j}-U_{CFD.j}}{U_0}\right)}^2}}{363}$$

(24)

$${ \tau }_p-Difference \% =\frac{{\sum }_{j=1}^5\sqrt{{\left(\frac{{\tau }_{p-M.j}-{\tau }_{p-CFD.j}}{{\tau }_{n-M}}\right)}^2}}{5}$$

(25)

Table 7. Difference between the measurements and simulation for all cases.

	Case 1	Case 2	Case 3	Case 4	Case 5	Case 6	Avg.
Air Speed
(m/s)	0.036	0.033	0.044	0.054	0.037	0.037	0.040
(%)	1.7%	3.2%	6.1%	1.8%	2.4%	3.7%	3.1%
Temperature
(°C)	0.29	0.35	0.42	0.27	0.30	0.28	0.32
(%)	3.6%	4.2%	6.0%	4.0%	4.9%	4.4%	4.5%
Mean Age of Air (τp)
(min)	2.3	0.9	0.7	0.8	1.0	1.0	1.1
(%)	9.0%	3.5%	2.6%	4.2%	5.3%	5.4%	5.0%

shows the difference between the CFD and the measurements for each case. The difference in the mean age of air was on average 1.8 min (5.4%), which was within the estimated accuracy of the tracer gas measurements. The numerical predictions showed the same pattern as the measurement of high ACE_p in the ceiling and mixed conditions in the occupied zone and without stagnation zones. The difference in the nominal time constant was between ±6% with an average of −1.5%, which was less than the accuracy of the tracer gas measurements. The error for velocity was on average 0.04 m/s (3.1%) and 0.32 °C (4.5%) for the temperature, which was close to the accuracy of the CTA anemometer. The difference became larger as the inlet velocity decreased and the cases that had oscillations also had the lowest inlet velocities (Cases 2 and 3). This was due to the difficulty involved in model jets passing a buoyancy driven flow with a non-transient model [51]. The greatest source of error in both Case 2 and 3 was the wall zone where the CFD model slightly overpredicted the velocities. The model did however show a good correlation in the manikin and symmetry zones for velocity as well as a good correlation of mean age of air for both Case 2 and 3. Figure 7 and Figure 8 show the velocity and temperature profiles close to the ceiling for both the CFD and the measurements. The model predicted the values and the slope of the curve very accurately in most cases. The profiles for the one-row cases (1 and 4) were identical to the behavior of WCJ. The velocity profile 0.6 m (√(r/d) = 4.7) looked similar to the WCJ in the free jet region (0 < √(r/d) < 5.9) [21] and the velocity profile 1.2 m (√(r/d) = 6.7) was very similar to a WCJ in the Coanda region (5.9 < √(r/d)<11) [21], see Figure 7. The profiles for the three row cases (3 and 6) were somewhat similar to profiles that have been measured for an impinging jet [15,43]. This could be because the third nozzle row was aimed at 30° angle (from the horizontal plane) towards the ceiling, which meant it could create a small impinging effect, whereas Cases 2 and 5 seem to have had velocity profiles somewhere between an impinging jet and a WCJ.

Figure 9 and Figure 10 show the results from case 5 where the CFD model predicted the two velocity profiles at ceiling level more accurately and therefore had better results in the wall zone. The numerical results captured the profile and average velocity for each horizontal measurement line in all three zones fairly well. It did however slightly overpredict the amplitude for each horizontal line, with more variation between each measurement point. This could be because the CFD slightly underpredicted the dissipation of momentum.

For all six cases, the CFD model underpredicted the temperature in the same two areas, at the floor level (0.1–0.9 m) in the manikin zone and in the symmetry zone. The underpredictions in these two areas were above the accuracy of the CTA (0.2 °C) and ranged from 0.3 to 0.7 °C. This underprediction was the reason why the average error of the temperature ranged between 0.3 and 0.4 °C, see Table 5. In all other areas the model captured the measurement results well, both in amplitude and pattern. There was, as with the velocities, a slight overprediction of the amplitude for each horizontal profile, see Figure 10 for an illustration of these trends for Case 5.

The numerical model slightly overpredicted the velocities in the wall zone in cases with low inlet velocities and the measurements showed that the room did not have a completely symmetrical temperature distribution. The overprediction of velocities and the underprediction of temperatures in some zones was most likely due to an underprediction of the dissipation of momentum, which led to slightly less thermal mixing in the model causing slightly lower temperatures. However, the model predicted the ceiling velocity and temperature profile well enough to investigate the effects of the number of nozzle rows on the near-field development. This and the good correlation both statistically and qualitatively show that the model was proficient enough to predict the indoor environment factors in the occupied zone necessary to perform a parametric study with regards to energy efficiency and IAQ.

4.3. Simulation Results

4.3.1. Jet Development

The number of rows had a profound effect on the jet development, especially on the entrainment of the surrounding air into the jet. Figure 11 shows the velocity profiles at different distances downstream. Case 4 had an inlet velocity that was three times higher than case 6 and therefore had higher velocities downstream. However, since the color in Figure 11 was normalized based on inlet velocities, it shows that the velocities declined at a faster rate in case 4 than in case 6. At a distance of 50d, case 6 had a maximum velocity of 0.44 m/s (44% of U₀), while case 4 had a maximum velocity of 0.85 m/s (28% of U₀).

Case 6 also conserved the inlet temperature (16.5 °C) better than case 4. At a distance of 50d, the average temperature of the confluent jet was roughly 20 °C in case 6, while in case 4 it was almost 2 °C higher (22 °C), see Figure 12. This was due to the lower entrainment of the surrounding air, because of a lower TI (due to lower velocities) because in Case 6 a coherent confluent jet was created faster (because of the smaller distance between the nozzles) and therefore the jet had a smaller “surface area” against the surrounding air. Both the temperature and velocity plots show that the confluent jet attached itself faster to the ceiling in case 6 than in case 4. This was because the two rows closest to the ceiling were aimed slightly upwards because of the curvature of the channel. The airflow from these two rows changed the angle for the confluent jet and the jet reached the ceiling faster. The confluent jet was also wider in case 6 than in case 4 because the starting array of the confluent jet was wider.

Figure 13 shows the mean age of air in cases 4 and 6; the entrainment of air into the confluent jet was much slower in case 6 than case 4. At a distance of 5d in case 6, part of the jet still had a mean age of air less than 0.25 τ_n, while the jet in case 4 was mostly above 0.50 τ_n. At 10d, most of the jet in Case 4 was close to 0.75 τ_n, so most of the air entrainment in Case 4 occurred before a coherent confluent jet was formed. At 10d in Case 6 however, a coherent jet was formed and because of the attachment to the ceiling a “protective layer” below the boundary layer of the ceiling was formed which slowed down the entrainment of the air attached to the ceiling. This effect continuedto 50d where the air close to the ceiling was still only 50% mixed in Case 6, while in Case 4 most of the jet’s air was above 0.75 τ_n.

4.3.2. Indoor Air Quality, Energy Efficiency and Thermal Comfort

The effects on the jet development also affected the temperature, velocity, and IAQ on the occupied zone of the room. Figure 14 shows the jet with the highest velocity in Case 4 when it reached the occupied zone. However, because of high entrainment the jet had been mixed more with the air from the occupied zone. The jet in cases 5 and 6 had a lower penetration (see velocity isoform at velocity = 0.3 m/s) and much lower entrainment with the air from the occupied zone, which resulted in the isoform for the mean age of air reaching much further into the room.

This increase in the penetration of air into the room resulted in differences in ACE for the different cases.

Table 8. ACE for each manikin, ε_T and AR in Case 1–6.

	ACE Manikins *				ACE BZ **	εT OZ ***	ACE-0.95	AR
Case	Left	Center	Right	Avg.
1	100%	102%	100%	100%	102%	104%	100%	0.009
2	105%	102%	99%	102%	103%	105%	99%	0.049
3	107%	107%	100%	105%	103%	105%	88%	0.139
4	99%	94%	88%	94%	94%	103%	40%	0.004
5	101%	100%	91%	98%	97%	104%	53%	0.021
6	101%	103%	98%	101%	102%	105%	100%	0.056
Avg.	102%	101%	96%	100%	100%	104%	80%	0.046

* = volume of hemisphere (r = 0.2m) at 1.1 m in front of manikin. ** = average ACE at horizontal plane at 1.1 m. *** OZ = 0.075 m from floor −1.8 m; 0.6 m to internal wall, ASHRAE Standard 62. ACE-0.95 percentage of area at 1.1 with ACE > 0.95. AR = Archimedes number.

shows that a larger array with lower momentum was more efficient in terms of IAQ. The three cases with a lower airflow, cases 1–3, had higher ACE than the other cases. It should however be noted that the mean age of air in cases 4–6 was still lower than in cases 1–3 because of the higher air change rate.

Figure 15 illustrates that the case with the highest inlet velocity (Case 4) had the lowest energy efficiency in terms of both heat removal and ACE in the part of the room away from the outlet. Case 6 with the same airflow, but one-third of the inlet velocity, had a much more even ACE in the BZ. Since the ACE in some of the cases and ε_T in all of the cases were above 100% it means that the system theoretically could use a lower air volume than a mixing system and still achieve the same comfort in the occupied zone. This would reduce the energy demand of the ventilation system since a lesser air volume would have to be heated/cooled and the power of the fans would be reduced.

As previously mentioned in Section 4.1.1 there were lower measured velocities below 1.3 m in the manikin zone and to a lesser degree in the symmetry zone. The vector field in Figure 16 shows that there was a stagnation point in the velocity behind the back of the manikins between the downwards airflow at the wall and the upwards airflow from the manikins. In the multiple row array cases (2, 3, 5 and 6), buoyancy forces (AR > 0.049) from the manikins were the main driving force in the occupied zone, moving the air up from floor level similar to displacement ventilation. In the one-row cases (1 and 4), however, the momentum from the jet was the main driving force (AR < 0.01), bending the vector field towards the jet, distorting the plumes from the manikins and drawing in the air from the occupied zone. This resulted in cases 5–6 having much lower temperatures (0.7 °C) than Case 4 in the occupied zone where the heat and old air were forced upwards in the plumes, but the heat and air were not directly entrained into the center of the confluent jet as show in Figure 12 and Figure 13. Cases 1 and 4 had no clear plumes because of the distortion from the momentum of the confluent jets, but the low airflow Cases 2 and 3 had more well developed plumes than high airflow cases 4 and 5 because of higher Archimedes numbers.

The high inlet velocity case (Case 4) also underperformed in terms of local thermal comfort compared to the five other cases. The five other cases had similar values for draft rate (DR) [52] and velocities, see

Table 9. ADPI, DR, and velocities (Vel. [m/s]).

Case	U0 [m/s]	Occupied Zone			Volume of Occupied Zone with
		ADPI	DR-Max	Vel.-Max	Vel. < 0.1	DR < 10 *	Vel.< 0.2	DR < 20 **
1	2.2	93%	22	0.26	69%	89%	98%	99%
2	1.1	95%	15	0.23	76%	95%	100%	100%
3	0.7	92%	17	0.23	80%	95%	100%	100%
4	3	97%	40	0.40	51%	56%	88%	89%
5	1.5	94%	22	0.24	71%	78%	98%	99%
6	1	90%	23	0.24	75%	81%	99%	99%

* = Category A in ISO-7730:2005. ** = Category B in ISO-7730:2005 [52].

. All cases had ADPI values above 90%, which indicates uniform thermal conditions in the occupied zone. Cases 1–4, 6 had very low velocities (<0.2 m/s) and a low DR (<20) in at least 98% of the occupied zone; their highest velocity was 0.26 m/s (Case 1) and the highest DR was 23% (Case 6) in the occupied zone. Case 4 had high velocities and a high DR close to the floor because the momentum of the high inlet velocity was conserved into the occupied zone. The maximum velocity was 0.40 m/s and the highest DR was 40% at ankle level in the occupied zone, which would cause thermal discomfort.

The variation of temperature at ankle level in case 4 was very small (<1.0 °C); therefore, the high DR values were not caused by uneven temperature distribution. There was however a high correlation between the DR and variation in velocities at ankle level in case 4, see Figure 17.

5. Conclusions

The main conclusions from the present study are:

• The $\overline{{\vartheta }^2}-f$ -model in combination with the DO radiation model is an adequate turbulence model to study the confluent jet ventilation (CJV) supply devices and showed few tendencies for numerical instabilities. The numerical predictions with $\overline{{\vartheta }^2}-f$ model had the best agreement with the experimental results, both statistically and qualitatively, and is computationally cost-efficient enough to run numerous cases needed for a parametric study.
• A larger array with multiple rows and a lower inlet momentum conserves the inlet temperature and mean age of air better than a single-row array with higher inlet momentum due to the confluent effect. Since the single row has a higher momentum because of a smaller inlet area, the velocities downstream are higher which leads to higher velocities at floor/ankle level in the occupied zone, even though the velocities decline at a faster rate for a single row array at 0–50d.
• The results show that the size of the array has a great impact both on near-field development and on the conditions in the occupied zone. The multiple row arrays had a higher IAQ in the occupied zone because the larger array conserves the mean age of air better. Because of lower inlet velocities, they also had lower velocities at ankle level, which decreases the risk of draft and thermal discomfort.

6. Future Work

Based on this work a parametric study is proposed to optimize the nozzle array with regards to heat removal effectiveness, indoor air quality, and thermal comfort. The parametric study should consider such factors as number of nozzles, length and width of nozzle array, inlet temperature, and airflow rate.