As seen in Figure 1, the geometry is generated to represent the office with cubicles. Note that the furniture, desks, and workers are modeled with very simplified geometry. This is common in indoor CFD simulation, because the primary interest is often the bulk air movement and temperature, and the microenvironment near objects and human bodies is of less importance. Consequently, the small-scale features and curvatures of the room objects can be overlooked. This particular office space geometry is generated using a CAD tool included in a commercial CFD package. However, it might be simpler for a practitioner to generate 3D geometry in an external CAD tool, especially if complex shapes and curvatures are required.

Figure 1. (A) Geometry and (B) Mesh

Part B of Figure 1 displays the mesh generated for the office once the geometry is finalized. The base size of cells in the room might be relatively large depending on the complexity of the flow regime, but a grid independence test should be conducted to determine that a balance between calculation time and accuracy has been achieved. In this case, a base cell size of 11.8 in. (roughly the size of the diffuser opening) is used for the automatic mesh generator. Polyhedral mesh is used here to reduce amount of mesh manipulation required (compared to hexahedral mesh). There are additional refinements near the region of interests, such as the inlet, the heat flux sources (lights, computers), and near the occupants. In this case, a 10% base cell size is used to refine the region near the occupants, in the interest of comfort level. The 11.8 in. base size was also checked for grid independence, where smaller base size did not yield a change in results in monitor points near the occupants (due to the refinement, the cell size is roughly 1 in. near the occupants, hence it is sufficient for the simulation). This cell size results in roughly 50,000 cells in the simulation domain, a number that can be accommodated by typical personal computers.

RANS with a k-ε turbulence model (Chen et al. 2010) is selected because bulk flow near the occupants is the primary interest. Due to the heat flux in the room (see the section on Boundary Conditions), the flow energy coupled solver is used to determine the buoyancy aspect of the driving force in addition to the forced convection from the diffuser. For the same reason, the simulation assumes ideal gas (or real gas if a specific property is known) due to the small temperature and pressure variance in the domain.

The specific boundary conditions for this example are the inlets, outlets, and heat sources in the room. The other boundary is a simple wall assumed to be no-slip. The outlet is a simple pressure outlet condition with a relative pressure of 0 in. of water compared to the room pressure (just upstream of the outlet) This is due to the inlet already providing the flow driving force, and results in a negative pressure at the outlet. The heat sources (people, computers, and lights) are given specific heat fluxes on their corresponding surfaces (see Chapter 18 of the 2017 ASHRAE Handbook—Fundamentals):

The supply diffuser model requires special consideration, because the inlet has a dominant influence on the room airflow pattern and velocity. Unlike the return/outlet, the supply/inlet must have both direction and magnitude in order to establish a correct flow field for the simulation. Figure 2 shows an example of using the momentum method (Srebric and Chen 2002) to model the square louver diffuser in this problem.

The diffuser momentum method was chosen because, if the square diffuser supply for the office was left as a square hole without any additional information beyond the designed flow rate, this inlet boundary condition would be assigned to a specific pressure differential or velocity inlet normal to the opening. This normal assumption is correct in terms of providing airflow rate into the room, but the initial velocity direction and subsequent flow pattern are completely incorrect: a diffuser would spread out the supply air for good mixing, whereas a square-hole inlet model generates a column of air with downward velocity. In Figure 2, the square inlet is divided into 16 smaller inlets, each assigned with a different direction and velocity magnitude, in an attempt to replicate the dispersion airflow field from the diffuser. The resulting airflow pattern should follow the diffuser manufacturer data (throw/spread) or any experimental data available. Supply air temperature was set to 57°F, a typical value in U.S. office buildings.

Figure 2. Illustration of Momentum Method for Inlet Model

In this method, because the air velocity and direction of multiple inlets are used, it is important to verify whether the total airflow rate from the inlets still represents the intended flow rate from the diffuser. If the momentum method is not applicable, Srebric and Chen (2011) also provide alternative diffuser modeling methods.

Convergence of the CFD solution was achieved by monitoring the solution residuals as well as several temperature monitor points in various part of the room (Figure 3). In this case, the solution’s mass and momentum equation residuals stabilized and dipped below 1 × 10⁻⁴ at roughly 450 iterations, and the temperatures stabilized at about 700 iterations. Selection of 1 × 10⁻⁴ in the continuity, momentum, and turbulence equation residuals is a typical rule of thumb for indoor CFD simulations, but 1 × 10⁻⁶ and 1 × 10⁻⁵ are recommended for the energy and species equations, respectively. These convergence tolerance metrics should be accompanied by stabilization of the solution variables. This stabilization indicates that a steady-state flow solution has emerged and the solution has been determined to be converged.

Figure 3. Temperature over Iteration Indicates Steady-State Convergence

To display the simulation results, two planes were used to show the solution contours near the occupants, displaying the air temperature and air velocity (Figure 4). A horizontal cut plane at head level of the sitting occupants demonstrates breathing zone quantities, and the vertical cut plane through the diffusers demonstrates inlet flow regime.

Figure 4. Temperature and Velocity Results from Simulation Shown in Two Planes Bisecting Region of Interest

The goal of the simulation was to determine whether the temperature and air velocity meet the ASHRAE comfort requirements. Per ASHRAE Standard 55, it is evident the velocity did not reach the level of a noticeable cold drift, and temperature near occupants is well within the ASHRAE comfort zone. Hence, this particular design fulfills the comfort criteria.

The modeled section of the open office is a 1250 ft² space with a 9 ft high ceiling, and exterior windows. The office has cubicle desks with partitions and cabinets (as well as computers). This section of the office is conditioned with six two-way chilled beams in the ceiling, and return air is assumed to leave through an opening on the back wall (Figure 5).

Figure 5. Office CFD Model: Simplified Geometric Model

Figure 5 is a simplified representation of the gross geometry of the office. The geometric features are

Because the purpose of the model is to study bulk airflow and temperature distributions, simplicity is key when representing these features. Although box representations are used for all furniture and computers in the space, extra detail has been added to increase the fidelity of the simulation. For example, the spaces under the desks are included in the model, and the computers are split into parts (monitor and tower). This is not generally required, but was done to more accurately capture the thermal plumes from the occupants and the equipment in the space. Because there are only 12 occupants, more refined models of the occupants are used, although simplified models such as cylinders would have also been acceptable. These specific cuboid models include the arms, legs, and a head, while retaining the total surface area of an average-sized, seated adult male (see Chapter 9 of the 2017 ASHRAE Handbook—Fundamentals). When modeling heat/contaminant sources (e.g., occupants), the surface area of the simplified representation is very important, because it determines the temperature of the surface when a heat flux boundary condition is specified and vice versa. The windows in the model include some details of the mullions, because one of the purposes of the CFD model is to study condensation risk. Geometric representations of the chilled beams and returns depend on the modeling approach as well as the desired level of accuracy. The chilled-beam model is described in detail in the section on Boundary Conditions.

A snapshot of the computational mesh in a horizontal and vertical plane is shown in Figure 6. A tetrahedral mesh is used in this example; it can efficiently handle complex geometries and is less sensitive to random changes in flow direction like the hexahedral mesh. It is also less labor intensive to generate for the geometry used in this example. Local mesh refinement is used around the various geometric features mentioned in the Overview of CFD Simulation Process section, with particular attention paid to local refinement around the heat sources (to capture thermal plumes), as well as the discharge slots of the chilled beams (to capture the jets) and, to a lesser extent, the returns. Mesh refinement perpendicular to all surfaces in the model is used to capture the thermal and hydrodynamic boundary layers. For the current turbulence model, 12 layers of prismatic elements are used at the walls, with a dimensionless wall distance of y⁺ < 2 for the first layer. In this simulation, a global mesh size of 5 in. is used, whereas mesh sizing of 2 in. is used around the computers, and 1 in. for the return/induction opening and the occupants. A mesh size of 1 in. is used on and around the discharge slots of the chilled beams to better capture the spread of the supply jet along the ceiling. A mesh growth ratio of 1:1 is used around the chilled-beam discharge slots, occupants, and computers, to ensure that this mesh refinement smoothly merges into the larger cell size of the room.

Figure 6. Classroom CFD Model: Computational Mesh Resolution at Vertical and Horizontal Planes

The boundary conditions for this example are

All radiation heat fluxes are assumed to be isotropic (independent of direction). The chilled-beam model itself is worth discussing in more detail. Specification of this boundary condition is divided into two parts: the discharge opening and the induction opening. To specify both, the following performance metrics of the chilled beam are required:

The induction opening is specified as an outlet boundary with a mass flow rate equal to the total volume of air induced by the chilled beam q_pK_in. Velocity measurements across the induction face of chilled beams show that the induction velocity varies very little, making a uniform velocity assumption reasonable. The discharge opening is specified as an inlet boundary with the total volume of air (induced + primary, q_p[K_in + 1]) set at the specified discharge angle. In the current case, the discharge angle is set to 15° and is based on smoke videos of the chilled beam from the manufacturer. The correct discharge is important; if it is set too low, it will cause the simulation to show throw (projection of discharge jet) as longer than in reality. Note that the velocity distribution along the length of the discharge slots is assumed to be uniform, which may not be realistic but is sufficient for the purposes of this simulation. Setting the appropriate discharge air temperature of the chilled beam requires some additional arithmetic to be implemented in the boundary condition. By applying a simple energy balance to the chilled beam, an equation for the discharge temperature T_m can be derived as a function of the preceding performance metrics:

(1)

where T_in is the temperature of the induced air, which is calculated by the simulation from the induction opening boundary. The sign convention for coil capacity is a negative value in cooling and a positive value in heating. Because the specific heat capacity of air varies very little between 55 and 80°F, it can be taken at standard conditions. Air density, however, should be taken as an average between the induced air and primary air temperature. It is also assumed that the cooling capacity of the coil is fixed when, in reality, it is a function of the water flow rate and the temperature difference between the water and the induced air (as well as the primary airflow rate). Depending on how much information is available from the manufacturer, this dependence can be incorporated into Equation (1).

A RANS approach is used with a finite-volume, fully coupled solver. The shear stress transport (SST) turbulence model is used to capture thermal effects near the walls with good accuracy for the velocity field as well (Zhang et al. 2011). A Monte Carlo model is used to capture the thermal radiation heat transfer in the space from all the heat sources. The buoyant force is computed directly from the density field and is important for modeling the cold downdraft from the windows, as well as accurately capturing the nonisothermal throw of the chilled beams.

Convergence is judged by a combination of the normalized RMS residuals of mass and momentum, as well as the global energy balance in the domain. A solution is considered converged when the mass and momentum residuals fall below 1 × 10⁻⁴ and the global energy imbalance is less than 5%. The global energy imbalance (total energy coming in versus going out of domain) is a very important convergence parameter, and adding the mixed flow temperature equation generally tends to slow down convergence. The criterion of 1 × 10⁻⁴ for the normalized residuals is a frequently used rule of thumb that gives good convergence and sufficient accuracy for the purposes of this simulation. The 5% criterion for the global energy balance is based on experience and is a compromise between solution accuracy and the ability to converge the solution to a smaller tolerance.

Figure 7. Vertical Temperature Contour Showing Cold Downdraft near Windows

Figures 7 to 10 show results for one of the configurations in the example. Contour plots of temperature, velocity, predicted mean vote (PMV), and contaminant removal effectiveness are shown. PMV is a thermal comfort metric that correlates the thermal comfort response of a large group of occupants to various flow field variables such as velocity, temperature, humidity, and thermal radiation, as well as occupant properties (e.g., clothing insulation, activity level). This metric is based on the seven-point comfort scale described in detail in ASHRAE Standard 55. The band between 0.5 and –0.5 on the PMV scale is considered to be acceptable. Contaminant removal effectiveness is defined as a ratio of contaminant concentrations (see Chapter 15 of Zhang [2004]) from a chosen surrogate tracer gas. It is one way of measuring how efficiently the HVAC system distributes the ventilation air; for other approaches, see Etheridge and Sandberg (1996).

The results show some nonuniformities in the space temperatures, with a noticeable cold zone near the window, most likely caused by the downdraft of cold air. The thermal comfort varies throughout the space as well, with a cold zone near the windows due to the downdraft. The streamline plot shows that the chilled beams do a good job of washing the windows, as well as the collision zone above the occupants. The contaminant removal effectiveness (CRE) is reasonable, although some stagnant zones (CRE < 1) are evident in the back cubicles.

Figure 8. Velocity Streamlines Showing Supply Air Velocity and Direction

Figure 9. PMV Contour Plot 42 in. Above Floor

The classroom is an 861 ft² space with a 10 ft high ceiling, an exterior window, and a radiant slab with two controlled zones (interior and exterior). The classroom has a projector and a teacher’s desk, as well as other furniture, including multiple cabinets and tables. The classroom is conditioned with three wall displacement diffusers on the side walls and two returns located on the ceiling near the exterior window. The various modeling decisions made for this specific case follow. Figure 11 shows a simplified representation of the gross geometry of the classroom. The geometric features in this CFD model are

Figure 11. Classroom CFD Model: Simplified Geometric Model

The spaces under the tables are included in the model because airflow moves along the floor in a displacement system. It is therefore important to avoid flow obstructions that do not exist in the real space. In addition, the surface of the table affects the dynamics of the thermal plumes of the seated occupants. There are 35 occupants, so simplified models are used to represent them. They are meant to capture the rough geometric shape of a seated occupant while retaining the total surface area of an average-sized seated adult male (see Chapter 9 of the 2017 ASHRAE Handbook—Fundamentals). It is acknowledged that there is some loss in accuracy, because the occupants are actually children. The standing teacher is modeled as a simple rectangular prism, with the appropriate surface area.

A snapshot of the computational mesh in a horizontal and vertical plane is shown in Figure 12. A tetrahedral mesh is used in this example; it can efficiently handle complex geometries and is less sensitive to random changes in flow direction like the hexahedral mesh. This is particularly important in displacement flows where thermal plumes (direction not known a priori) dominate the flow. It is also less labor intensive to generate for the geometry used in this example. Local mesh refinement is used around the various geometric features, with particular attention paid to local refinement around the heat sources (to capture thermal plumes), the displacement diffuser (to capture the three-dimensional floor jet), and, to a lesser extent, the returns. Mesh refinement perpendicular to all the surfaces in the model is used to capture the thermal and hydrodynamic boundary layers. For the current turbulence model, 12 layers of prismatic elements are used at the walls with a dimensionless wall distance of y⁺ < 2 for the first layer. In this simulation, a global mesh size of 9 in. is used, and mesh sizing of 2.5 in. and 2 in. is used around the return and the occupants, respectively. A mesh size of 1 in. is used on the displacement diffuser to capture its waterfall flow pattern as well as the complex three-dimensional thermal wall jet. Refinement around a displacement diffuser is typically coarser than a mixing diffuser at the same flow rate, as a result of the shorter projection distance and the velocity of the jet. A mesh growth ratio of 1:15 is used around the displacement diffusers, the occupants, and the projector, to ensure that this mesh refinement smoothly merges into the larger cell size of the room.

Figure 12. Classroom CFD Model: Computational Mesh Resolution at Vertical and Horizontal Plane

The boundary conditions for this example are

All radiation heat fluxes apart from the window are assumed to be isotropic (independent of direction). The displacement diffuser model is worth discussing in more detail. Although the face velocity of the diffuser is low and flow in the space is driven by the heat sources, it is important to accurately capture the velocity temperature distribution in front of the diffuser for an accurate assessment of thermal comfort. The model used in this example is a simplified version of the momentum method (Chen et al. 1999) that does not take the perforated face into account. In situ measurements have shown this to be a reasonable approximation for this type of diffuser. Note that better accuracy can be achieved with a more refined model.

A RANS approach is used with a finite volume, fully coupled solver. The shear stress transport (SST) turbulence model is used to capture thermal effects near the walls with good accuracy for the velocity field as well (Zhang et al. 2011). A Monte Carlo model is used to capture the thermal radiation heat transfer in the space both from the short-wave (solar) and long-wave (e.g., people, lights) heat sources. Thermal radiation heat transfer from the ceiling to the cool floor also has a strong effect on the vertical temperature gradient in the space. The buoyant force is computed directly from the density field and is a crucial component of the simulation, because the dynamics of displacement ventilation systems are driven primarily by buoyancy.

Figure 13. Vertical Temperature Contour Showing Stratified Temperature Distribution Typical of DV Systems

Figure 14. Velocity Streamlines Showing Supply Air Velocity and Direction

Figure 15. PMV Contour Plot 42 in. Above Floor

Convergence is judged by a combination of the normalized RMS residuals of mass and momentum as well as the global energy balance in the domain. A solution is considered converged when the mass and momentum residuals fall below 1 × 10⁻⁴ and the global energy imbalance is less than 5%. The criterion of 1 × 10⁻⁴ for the normalized residuals is a frequently used rule of thumb and gives good convergence and sufficient accuracy for the purposes of this simulation. The global energy imbalance (total energy coming in versus going out of domain) is a very important convergence parameter for displacement ventilation problems, because the thermal field takes much longer to stabilize than the hydrodynamic field. For example, the mass and momentum equation residuals can stabilize and converge before the global energy balance reaches its converged state, which yields the wrong solution (no thermal stratification).

Figure 16. Contaminant Removal Effectiveness (CRE) 42 in. Above Floor (Seated Breathing Height)

Figures 13 to 16 show results for one of the configurations investigated in the study. Contour plots of temperature, velocity, predicted mean vote (PMV), and contaminant removal effectiveness are shown. The contaminant removal effectiveness is defined as a ratio of contaminant concentrations (Zhang 2004) from a chosen surrogate tracer gas. It is one way of measuring how efficiently the HVAC system distributes the ventilation air. The results show good thermal stratification in the space with a layer of cold air spread across the floor. The streamline plots imply that the specific arrangement of diffusers does a good job distributing the cool supply air across the floor and around all of the flow obstructions in the space. The PMV plot shows a very distinct hot spot under the window, which is a direct consequence of the magnitude and direction (pointed down) of the solar load. Overall, the thermal comfort in the space is slightly outside the desired comfort range (on the warm side). Contaminant removal effectiveness is also not very uniform, suggesting that some improvements can be made to the air distribution system to better distribute the supply air.

In contrast with many general-purpose CFD tools, which require the user to assemble models using simple geometric primitives, data-center-specific CFD modeling tools and electronics-specific modeling tools typically provide high-level building-block objects that only need to be configured to represent the geometry and operating conditions of a given data center. For example, racks, coolers, or tiles can be selected from a library of objects and simply placed into the room.

Figure 17. Data Center Layout

To best capture the physics that govern a given data center’s airflow patterns, it is recommended that the floor plenum be modeled together with the whitespace. However, it may be simpler to model the two spaces separately, and in cases where the floor tiles are fairly restrictive (in the range of 25% open area or less), doing so delivers reasonable accuracy.

The example data center was modeled using a rectangular Cartesian grid, which maps well to typical data center geometry. Generally, cells of side length 6 in. or smaller are recommended (Zhang 2008) for modeling data centers. In the example case, grid cells in the whitespace were set to have a maximum side length of 6 in. and a minimum side length of 1 in. below the dropped ceiling. The variation in grid cell dimensions allows for finer grid cells near objects and complex geometry, and larger grid cells in open space. However, the model uses grid cells of up to 24 in. on a side within the ceiling plenum, which is largely open space and far away from features of interest.

The example case uses a RANS CFD solution methodology, assuming steady-state conditions, and the standard k-ε turbulence model was employed.

Several specific modeling choices resulted in greater accuracy of the example model. Because the room is located inside a larger indoor enclosure and there is little heat transfer across the room walls, walls were modeled as adiabatic solid boundaries. Additionally, rather than modeling each rack as a monolithic box with uniform airflow into the front face and out of the rear (as per Zhai et al. [2012]), each rack was subdivided into a number of slices (Pardey et al. 2015), as shown in Figure 18. Each slice was assigned its own airflow rate and accordant temperature increase based on its power, with the assumption that servers draw 0.04 cfm per Btu/h. Blanking panels were modeled as solid blocks, and open spaces were modeled as open. Rack doors were modeled as two-dimensional flow resistances coincident upon the front plane of the rack. A vertical momentum source was enforced inside a volume region over each tile extending to a height of 4 in., in the method of Abdelmaksoud et al. (2010), to accurately model the jet-like tile airflow. Inside the floor plenum, care was taken to explicitly model the vertical stanchions that support the raised floor. However, stanchions can equally be modeled as a distributed resistance, as discussed in VanGilder et al. (2016).

Figure 18. Rack Model

Individual perforated floor tile airflow rates were measured using the anemometer scaling method described in VanGilder et al. (2016). The floor of the example facility is bolted and gasketed, preventing any significant air leakage, so the sum of floor tile airflows was considered to represent the total room airflow. To simulate airflow patterns inside the floor plenum, the total room airflow was divided among the 21 airflow inlet bays that feed the plenum, and each bay was modeled as a fixed-flow boundary condition. The portion of the total room airflow allocated to each bay was calculated based on measured air velocity through each given bay. It is important to note that total room airflow could also be estimated by summing the nameplate flow rates of all active computer room air-handler (CRAH) units, or by measuring airflow from CRAHs using an anemometer or a flow hood.

The model is considered converged when residual errors in the mass, momentum, energy, and turbulence-model equations reach a predefined level, which is usually set to a default value. Here, in accordance with best practices, a number of simulations were performed using computational grids with varying cell sizes. By comparing the predicted temperatures and airflow patterns produced by simulations with different cell sizes, grid independence was established.

The CFD model was created to help optimize performance of an existing data center, so greater detail was required than would be necessary for preliminary design purposes. As such, the model was calibrated by iteratively increasing the detail and fidelity of its construction based on measurements and observations. Key parameters for accuracy were inlet rack temperature and tile airflow rate. Rack inlet temperatures were measured at four equally spaced points along a vertical line at the center of the face of each rack.

The accuracy of the model was enhanced by dividing each rack into 42 individual one-unit slices, each of which was either occupied by an active piece of IT equipment or blanking panel or was left open. Additionally, because IT equipment often generates jets of warm exhaust air, it was important to account for any IT equipment that produced such jets directed into a given cold aisle (either inadvertently or because of suboptimal equipment design). Including the doors of racks in the model as two-dimensional flow resistances also improved accuracy. In the floor plenum, accuracy was increased by modeling the vertical support stanchions, because it was discovered that if the stanchions were not included in the CFD model of the plenum, the model produced substantially different tile airflow predictions.

Capturing accurate input data is equally as important as making appropriate modeling choices, and to this end, it was crucial to account for the resistance of the flow hood on the measured perforated tile airflow rates. Methods for doing so have been evaluated at length in Zhang (2004). Note that CRAH return temperatures to the cooling system could also be used as a criterion by which to evaluate the accuracy of the CFD simulation or check the consistency of other assumptions.

Figure 19. Comparison of Measured and Predicted Tile Airflow Rates

Figure 20. Comparison of Measured and Predicted Rack Inlet Temperatures

The model was created specifically for managing and optimizing an existing data center. Therefore, it was necessary to model in more detail than is required for preliminary design purposes. Figures 19 and 20 show the results of predicted versus measured perforated tile airflows and rack inlet temperatures, respectively. Figure 19 shows two variants of the plenum simulation: one that includes raised-floor support stanchions and one that does not. The model that includes stanchions predicts 186 of 192 tiles to within ±10% of the individual measured airflows and the remaining six tiles to within ±20%. Figure 20 depicts predicted rack inlet temperatures according to their accuracy and shows that the calibrated model predicts 105 of the 115 active racks with an error rate of less than 10% of the maximum temperature range observed in the data center. To achieve this level of accuracy, it was important to model the contents of each rack at the unit (U) level, rack doors as flow resistances, momentum sources above each perforated floor tile, and any backwards-oriented IT equipment. It will be necessary to recalibrate the model periodically to account for any changes to IT equipment population, room arrangement, or general operating conditions.

Melikov and Kaczmarczyk (2007) discuss the importance of detailed indoor objects such as a human body on indoor airflow characteristics, and indicate the local impacts of most details of indoor objects. Focusing on the general indoor airflow patterns and interactions between patient and medical staff, they simplified the simulation of indoor subjects such as human bodies and equipment (except for surgical lighting) as rectangular geometries with exact heat sources as tested. This practice facilitates the generation of high-quality meshes and therefore improves both speed and accuracy of the simulations.

A rectangular Cartesian grid, which maps well to typical OR geometry, was used. Local grid refinement was implemented near critical spaces and objects such as walls, inlets, and persons. The results of a CFD simulation are highly dependent on the quality of the computational grid. The grid refinement study was conducted on the following grids: 70 × 58 × 45 (180 k cells), 87 × 73 × 57 (362 k cells), 106 × 91 × 70 (675 k cells), 124 × 111 × 86 (1.2 million cells), and 155 × 142 × 108 (2.4 million cells). Figure 22 demonstrates the finest grid distribution.

Both RANS and LES CFD methods were tested for this example case. Although advanced CFD modeling techniques such as LES provide substantial benefits, the currently available RANS technologies have proven to be adequate for modeling the steady-state characteristics of the hospital operating room air distribution. In the RANS CFD solution methodology, the RNG k-ε turbulence model (Yakhot and Orszag 1986) was used, as suggested by Zhang et al. (2007).

Figure 22. Grid Refinement Case: 7.8 ft Cells

Most indoor objects such as persons and equipment were specified straightforwardly using the standard wall/block boundary condition methods. Inlet boundary condition is critical to accurate CFD modeling of indoor environments, because this is the primary source of momentum responsible for overall room air distribution pattern. Srebric and Chen (2002) performed a comprehensive analysis of diffuser boundary conditions to determine appropriate simplified boundary conditions, and found the box and momentum methods to be the most appropriate for the diffusers that were tested. The momentum method was used in this example, based on the recommendation of Chen and Srebric (2012) for the grille diffuser that is similar to the nonaspirating diffuser type.

The simulation was considered converged when the sums of residual errors in the mass, momentum, energy, and turbulence-model equations reached a predefined level (0.1%). The grids of different sizes were evaluated using the normalized root mean squared error (NRMSE) of the CFD model results with different grids (Wang and Zhai 2012). Figure 24 shows the NRMSE of the predicted x and y direction velocity at the four measure poles (1 to 4) across the center axis of the room, 9.45 ft high (Figure 23), between the 180k and 362k meshes and the 675k mesh. It reveals that there is generally a great improvement in error with the 362k mesh; the computational error is typically below 10%, and absolutely below 30%. Based on this, and to minimize the simulation time, the 362k mesh was chosen for various parametric simulations.

Figure 23. CFD Grid Refinement Measurement Locations in Central Cross-Sectional Plane

The simulation replicates the airflow pattern as observed in the lab (see Chapter 13 of the 2017 ASHRAE Handbook—Fundamentals): an inward curvature of the airflow to the center of the jet stream, as seen in Figure 25. This behavior reduces the overall coverage area and could pose a contamination risk to the patient.

Figure 24. NRMSE Comparison Between 180k and 362k Meshes and 675k Mesh

Figure 25. 

The quantitative comparisons of simulation and experimental results were plotted in Figures 26 and 27, for x and y velocity components, respectively. Figures 26 and 27 show that the CFD simulations closely follow the experimental results, with a few exceptions. It also appears that there is, in general, a large difference between the experimental results and the 180k mesh, but a smaller difference between each of the latter meshes.

This example demonstrates the applicability of CFD for modeling and analysis of airflow in the surgical environment. Although CFD can be accurately used for modeling indoor air distribution and contaminant transport in operating rooms, the CFD user must be extremely careful in implementing these models, to ensure accurate simulation of airflow. The sensitivity of airflow to thermal characteristics of the indoor environment makes the model sensitive to heat gain input parameters. The heat gain and inlet boundary conditions must be carefully selected to ensure that the resulting air distribution patterns are correct.

Figure 26. Comparison of U-Velocity in X Direction

Figure 27. Comparison of W-Velocity in Z Direction

The following modeling methods were found to adequately model the represented physics identified in the operating room air distribution.

The general indoor environment conditions place the operating room indoor air distribution in the mixed convection category, but high cooling loads can lead to a strongly buoyancy-driven flow, verified by the parametric study of the Archimedes number of the supply air jet in the OR. The study reveals that the dependence of the room air distribution on the Archimedes number of supply air jets, rather than face velocity of supply diffuser, is of significant importance.

Both the geometry and the mesh were generated using the preprocessing tools available in the software package. When modeling natural ventilation in CFD, it is essential to accurately model the effects of flow through the openings. One approach is to model some proportion of the external space, thus removing the need to specify boundary conditions at the opening boundaries. However, this requires more computation and detailed information about the geometry and components that comprise the inlets and outlets. As an alternative, the boundary here was defined at the point at which air enters and leaves the auditorium, and conditions specified to represent the effects of components such as louvers, dampers, attenuators, and heating elements at openings.

The mesh used in this simulation was a structured Cartesian grid. Such meshes can reduce the computational overhead required by unstructured meshes and avoid some of the numerical instabilities during the iterative solution procedure.

Predicted pressure drops across components at inlets and outlets were used in the orifice flow equation (see Chapter 16 of the 2017 ASHRAE Handbook—Fundamentals) to determine a discharge coefficient C_d to represent resistance to flow at openings. For any openings where pressure loss data were unavailable, a discharge coefficient of 0.6 was used, representing the effect of a sharp-edged orifice. In these cases, it was necessary for the design team to ensure that the combination of free opening area and discharge coefficient in the final design provided at least as much airflow as the combination used in the CFD simulations.

Heat gains are generated by occupants (205,000 Btu/h) and lighting (136,000 Btu/h). These were modeled as convective gains into the domain. However, to account for the convective/radiative split, 50% and 10% of the occupant and lighting gains (respectively) were assumed to be radiated to the surrounding surfaces from where they were convected into the domain. Buoyancy was modeled using the Boussinesq approximation (Gray and Giorgini 1976), which assumes density to be constant everywhere except in the source term of the momentum equation.

Turbulence is modeled using the k-ε model with constants derived from renormalisation group theory (ASHRAE Standard 55). This model was used, because it was found to more accurately predict entrainment into the buoyant plumes, which are known to determine the key characteristics of these flows (e.g., interface height, temperature of stratified layer) (Cook and Lomas 1998). Validation of the CFD techniques described here was conducted using analytical and experimental models as described in Cook and Lomas (1998).

Convergence was deemed to have been achieved when the residual in the enthalpy equation, in watts, was less than 1% of the total heat convected into the domain at the heat sources, and when the absolute values of variables at the user-defined monitoring point did not change by more than about 0.1% over approximately 20 iterations. The monitoring point used in this case was located just below the roof outlets in the auditorium. Convergence was successfully achieved according to these criteria, using under-relaxation in the form of false time steps. This type of convergence control uses the local cell size and a characteristic time scale for the flow evolution to define under-relaxation factors, which vary throughout the domain. False time step values of 0.1s were used for each of the three momentum equations.

The simulation predicted a stable stratification in which the occupied zone is bathed in the air at, or just above, the ambient temperature (see Figure 29). This illustrates one of the key advantages of buoyancy-driven displacement ventilation, whereby the height of the stratified layer is determined by the ventilation opening sizes, not by the amount of heat gain in the space. The system can be thought of as self regulating because any increase in heat gain increases the temperature of air in the stratified layer, which increases the driving force, and hence the flow rate through the building, as required. In this case, the layer of stratified air drives a flow of 8.3 ach through the main auditorium, equivalent to about 45 cfm per person. This is typical for naturally ventilated buildings of this type and is necessary to provide thermal comfort, as well as an adequate supply of fresh air.

Figure 29. Temperature Prediction over Vertical Plane in Auditorium

This example includes two actual warehouses: the QT and KS buildings. The QT building is a storage warehouse with a size of 39 × 30.5 × 22 ft (L × W × H), with about half of the space divided horizontally into two sections (an upper section and a ground section) to maximize storage space, as shown in Figure 30A. The existing QT warehouse was equipped with a 51,000 Btu/h electronic forced-air heater beneath the ceiling, and a 670 cfm ceiling bucket fan. The current analysis added a 160 rpm ceiling rotating mixing fan for the evaluation of potential destratification improvement.

The KS warehouse (Figure 30C) is a larger building, 135 × 57 × 23 ft, heated by a mechanical duct diffuser system and six 3400 Btu/h baseboard heaters at the ground level under the windows. The conditioned air from every duct diffuser is projected horizontally into the space at 100 cfm per outlet, at a constant 104°F. Six 670 cfm bucket fans were installed above the baseboard heaters for destratification.

Figure 30. View Models and CFD Models of Warehouses

General geometrical information was collected, including the dimensions of the warehouses, the size of the doors and windows, connections to neighboring rooms or buildings, and storage rack layouts and dimensions. Here, a commercial CFD software package was used. Most of the storage racks and shelves were modeled by solid blockages, but the mixing devices and HVAC ducts were modeled as their original cylinder shapes. Internal partitions were modeled as zero-thickness blockages because no heat transfer calculation was needed for internal structures. The ceiling rotating fan was modeled in detail to reflect the real-world unit, considering the fan blade diameter of 4.3 ft and a curve surface with a 1.15 ft radius arc; the chord of each blade has a 15° angle with the horizontal plane (Momoi et al. 2004). Figure 30E shows the CFD model of the ceiling rotating fan, and Figure 30F, the bucket fan.

Both warehouse models are meshed mostly by Cartesian cut-cell grids first in the software. For regions with heaters and fans, more detailed localized meshes were created to provide enough resolutions considering the existences of strong gradients. Figures 30E and 30F show the local mesh examples for the ceiling rotating fan and the bucket fan models. To model the rotating fan, the whole model was divided into two parts: the cylindrical section encompassing the fan and rotating at a speed of 160 rpm, and a stationary section, as shown in part E of the figure. The bucket axial mixing fan was modeled by a constant-pressure jump fan model at a set volumetric flow rate (Figure 30F). The forced air heater (not shown) was also modeled by the pressure jump fan model with a constant heat source, following the method proposed by Forrest and Owen (2010). For more detailed CFD setups, see Wang and Li (2017). For the KS building, similar mesh refinement treatments were applied to the baseboard heaters and near the HVAC duct supplies, as shown in Figure 30D. As a result, a total mesh of 7.5 ft was created for the QT and KS buildings.

To model buoyancy-driven flows in such large spaces, a real gas model with temperature varying density was used. For the turbulence effects, the standard k-ε two-equation turbulence model has been found effective by previous studies (Li et al. 2009). Momoi et al. (2004) also found the standard k-ε model adequate for modeling ceiling rotating fans.

The thermal resistance properties of building envelopes were collected from the owners or engineers. For CFD simulations of large spaces, various thermal boundary conditions are available for modeling solid surfaces: constant temperatures, heat transfer coefficients, and heat fluxes. For modeling actual warehouse envelopes, the constant heat transfer coefficient boundary conditions were more reasonable for modeling both interior and exterior surfaces. The estimation of the convective heat transfer coefficients h_c was based on the models suggested by Qi et al. (2013) as a function of temperature differences, surface orientations, and Reynolds and Prandtl numbers. For other surfaces (e.g., the ground), constant temperatures were applied. For the duct diffusers, constant velocity and temperature were modeled using the data available from the on-site HVAC monitoring systems. Table 4 provides the key boundary condition settings.

Figure 31. Temperature Validation for Two Locations in QT, and Velocity Validation for Modeled Ceiling Rotating Fan Using Literature Data (Momoi et al. 2004)

A grid sensitivity study was conducted for the QT warehouse using 7.5 ft cells and 10.1 ft cells. For the KS warehouse, the total grid number was also chosen to be 2.3 million.

Different thermal destratification strategies were simulated, as shown in Table 5. For the QT warehouse, two more bucket fans were added: a second fan, to the left of the existing, first fan, and a third fan, close to the wall, with the supply air pointing upwards towards the ceiling. A ceiling rotating fan was also added, as shown in Figure 30B. Validation results of temperature and velocity profiles are shown in Figure 31. Including the existing system q₁, a total of six different combinations of fan operating strategies were simulated from q₁ to q₆. For the KS warehouse, different bucket fan operating speeds (all off, 30%, and 100%) and air supply directions from the duct diffusers (horizontal throw, vertically downward supply, and 45° downward supply) were simulated by a total of five cases from k₁ to k₅ (Table 5).

To quantify the level of thermal stratification, the non-uniformity coefficient θ is defined with Equation (2): a greater value of θ indicates a stronger stratification.

(2)

where

Figure 32A compares the nonuniformity coefficients for all six cases of the QT warehouse. Comparing the cases with bucket fans running (q₁, q₄, q₅, and q₆) and off (q₂) demonstrated a clear improvement in mixing when the fan was running. For the bucket fan flow direction, pushing air downwards (q₁) was better than sending air upward through the bucket fan (q₄): θ ≈ 0.23 in q₁, and θ = 0.38 in q₄. Therefore, it is more beneficial to use a bucket fan close to the roof heating source with downward air than a wall-mounted bucket fan blowing air upwards, which also causes a potential increase of ceiling heat transfer. The q₁ case was the worst thermal stratification case: the resultant roof heat loss was found to be about 9900 Btu/h, as shown in Figure 32B. Applying a destratification strategy (e.g., the best case of q₃), the energy saving can be 3800 Btu/h in q₃, a reduction of around 40% compared to q₁. Energy savings from other strategies also found that a lower θ result in a less heat loss through the roof.

Figure 32. Comparison of Nonuniformity Coefficients, Heat Losses, and Savings in Warehouses

The calculated nonuniformity coefficients of the KS warehouse are all less than 0.1, ranging from 0.006 to 0.096 (not shown here). This shows that, although the KS warehouse is about seven times larger than the QT, using ceiling duct diffusers and many wall-mounted bucket fans significantly overcomes the thermal stratification. The corresponding energy saving from ramping up fan speed is only 2% for k₁ and 5% for k₃ when using k₂ as the baseline, as shown in Figure 32B. Therefore, using duct diffusers already ensures a good mixing, and increasing mixing fan speed is unnecessary in this case. A more effective way to further reduce thermal destratification is to change the diffuser air supply angles (at zero extra energy costs). When the diffuser supply air is changed from horizontal throw in k₁ (the existing system) to 45° downwards in k₄ and, further, to vertically downward (90°), the corresponding energy saving is 5% in k₄ and over 6% in k₅ (a similar scale as, or slightly better than, adjusting bucket fan speed). Therefore, for both warehouses, a preferred thermal destratification strategy is always to actively deliver warm conditioned air effectively to the lower portion of the building, rather than passively relying on local mixing devices (e.g., wall mounted bucket fans), which are often complicated by issues of location and number of units used. An exception is the ceiling-mounted rotating fan in the QT warehouse, which was found to be a superior destratification method. This example demonstrated that using CFD techniques can effectively aid the design and analysis of thermal stratification and destratification in warehouses.

The building (Figure 33) is based on a representation of a U.S. Department of Energy (DOE) Commercial Reference Building model developed by Ng et al. (2012). The building measures 165 by 108 ft. Floor and plenum heights are 9 ft and 4.125 ft, respectively, for a total of 39.37 ft. It consists of three floors, each served by a return air plenum located above the occupied level.

The multizone representation of the building is shown in Figure 33. Each floor was subdivided to include five occupied zones: a core zone and four perimeter zones. This zoning strategy reflects that of the original DOE energy model upon which it is based, which treats the perimeter zones separately due to orientation-based solar loads and the core zone having a different energy load from the perimeter zones. This zone topology is also practical for the purposes of this example, because it captures the wind-related effects on the perimeter zones and the core/perimeter ventilation system layout. The core zone was slightly modified from the original model so that each floor includes a stairwell, elevator, and restroom. The stairwell and elevator span the entire height of the building to allow consideration of stack-driven flows through these vertical shafts.

Figure 33. Medium Office Building Model: (A) Schematic Floor Plan and (B) 3D Representation

Airflow characteristics of the building envelope and interior partitions are based on the model described in Ng et al. (2012). Exterior envelope leakage is defined using an effective leakage area of 0.08 in²/ft² at a reference pressure of 0.08 lb/ft², a discharge coefficient of 1.0, and a pressure exponent of 0.65. The leakage between the perimeter and core zones are defined as open office plans to reduce interzone resistance to airflows, and the ceiling plenum was connected to the occupied zones with orifice areas of 10 ft² per 1000 ft² of ceiling area. Wind pressure coefficient profiles were specified based on low-rise buildings, as defined by Swami and Chandra (1988), and the building was situated in a suburban location to yield a surface-averaged wind speed modifier of 0.36 (Walton and Dols 2006). Simulations were run using a TMY3 weather file for Baltimore, MD.

Each floor of the building is served by a single air handler, and each restroom is served by a dedicated exhaust system. Two ventilation systems were modeled in this example: a constant-air-volume (CAV) mechanical ventilation system and a CAV system that also incorporates CO₂-based demand-controlled ventilation (DCV). The restrooms were exhausted at a constant rate of 212 cfm in both systems.

Contaminants were simulated to represent both indoor and outdoor sources. The selected contaminants enable the investigation of the effects that various ventilation system and building leakage characteristics have on different types of sources. Occupant-generated carbon dioxide (CO₂) was included to model demand-controlled ventilation, because it can be used as an indicator of people being present in the building and thus as a controlled variable for establishing ventilation airflow rates. Volatile organic compounds (VOCs) were simulated to represent a generic indoor source associated with building materials and occupant activities. Particulates less than 2.5 μm in diameter (PM_2.5) were simulated to represent an outdoor source. The properties of these contaminants, sources, and sinks are based on those provided in Ng et al. (2012).

Multizone modeling provides the ability to evaluate ventilation system designs with respect to building envelope tightness using pressure-based airflow network calculations as opposed to assumed and empirically determined air change rates (Ng et al. 2014). Figure 34 provides plots of the CO₂ concentration in the three return plenums. The plot on the left is for the base configuration with a building leakage rate of 0.08 in²/ft², and that on the right is ten times tighter. As shown in Figure 34, there is less variation in the CO₂ concentrations among the three zones in the tighter building, enabling better control of the DCV system.

Calculating airflow and contaminants together allows comparisons between various ventilation systems that use contaminant-based control strategies, evaluating the controlled contaminant (here, CO₂) and the effect of different control schemes on uncontrolled contaminants (here, VOCs and PM_2.5). Figure 35 shows the building average VOC concentration for four different combinations of building envelope leakage and ventilation system type: leaky envelope with CAV system, tight envelope with CAV system, leaky envelope with CAV/DCV system, and tight envelope with CAV/DCV system. This plot reveals that the internal concentration is elevated by the lower outdoor air intake rates attributed to the DCV control system for both the leaky and tight envelopes. Thus, even though DCV and tighter building envelope may lead to energy savings by reducing outdoor air intake, this approach could lead to increased occupant exposure of internally generated, uncontrolled contaminants.

Figure 34. CO₂ Concentration for DCV System in (A) Leaky and (B) Tight Buildings

Further, as mentioned above, building energy measures such as ventilation strategies and envelope leakage rates can affect building energy usage. Multizone energy simulation tools could be used to evaluate the trade-offs between energy conservation strategies and IAQ. Multizone, whole-building energy analysis software can provide a wide range of building energy usage characteristics and has been coupled with various means, including cosimulation, to allow these analyses to be carried out in an integrated manner (Dols et al. 2016; Ng et al. 2018).

Figure 35. Comparison of VOC Concentrations with Respect to Envelope Leakage and Ventilation System