|Home | About | Journals | Submit | Contact Us | Français|
This study examines the usage of computational fluid dynamics (CFDs) for estimating the time-elapsed decay of contaminants within a chamber experiencing high Reynolds flow. CFD results were compared with measurements taken at a controlled facility. In addition, parameters of the CFD simulation were examined; namely the effects of turbulence and inertial transport at high Reynolds number ventilating flows, as well as inlet duct configuration and its effect on the inlet velocity profile. The agreement between the computational and experimental clearance times was quite good, with percent errors as low as −5.32% at high flow rate and −11.8% at the lower flow rate. This study determined that for high Reynolds flow, diffusive transport effects may be ignored as the majority of mass is transported via the bulk stream, i.e. momentum transport. In addition, resolving the inlet velocity profile was of prime importance for accurate simulation of ventilating flows and prediction of contaminant washout. This was done by including the inlet duct geometry in the computational domain. In addition, it was found that despite different flow rates, the predicted contaminant washout took ~12–13% longer than predicted assuming instantaneous mixing. Furthermore, percent error between computational and experimental data as low as −5.32% shows that CFD is a useful tool for studying ventilation phenomena.
The study of the transport and decay of hazardous gases within ventilated spaces is an important problem that has wide implications in the areas of public and occupational health. Understanding the behavior of contaminants within ventilated spaces and the rates at which they decay to a ‘safe’ exposure level is integral to maintaining occupational health. Typically, contaminated indoor air is diluted with fresh air to improve air quality, which results in the exponential decay of the contaminant species. Examination of the transport and decay of hazardous contaminants is vital to understand exposure risk as well as to design effective ventilation systems for these environments. Contaminant sources can originate with a room or be transported through ventilation systems. For these situations, it is important to understand how effectively a ventilation system replenishes contaminated air with clean air. A typical office space may have an air exchange rate of ~6 air changes per hour (ACH), Peyton and Abdou (1995), but many businesses and industries maintain airflows above 20 ACH such as bars, repair garages, commercial kitchens, and manufacturing plants. The expression ‘air changes per hour’ is essentially a measure of how quickly clean air replaces contaminated air within a space, Peyton and Abdou (1995). In the case of contamination of an occupational setting, a rapid washout of the contaminants is highly desirable. The rapid washout is achieved by washing in fresh air at high flow rates, 19 ACH and 34 ACH, resulting in high Reynolds number flow. Reynolds number linearly increases with increasing flow rate, thus for the two flow rates examined the Reynolds number is ~3.2 and 5.7 times typical room Reynolds numbers. In addition, the critical Reynolds number for a round free jet is about Re = 30, Schinlichting and Gersten (2000), which is easily exceeded in this study. For practical purposes, the flow rates studied can be deemed high Reynolds flow.
Traditionally, predictions of ventilation trends within an enclosure were obtained via experimental procedures that typically involve the use of a tracer gas and measurement devices to monitor its concentration at various locations. This practice of prediction by experimentation can be costly and time consuming, as it requires construction of a suitable experimental chamber or scale models, instruments to measure the tracer gas, and personnel to monitor the process. However, with the advent of enhanced computing technologies and the availability of commercial codes, computational fluid dynamic (CFD) techniques have been presented as cost-effective alternatives for predicting airflow trends in ventilated spaces; commercial codes have an annual cost of ~$620 per license. With the use of CFD methods, one can quickly and inexpensively look at a range of possible parameters and inputs to the problem and determine how they affect the overall ventilation.
Several earlier studies have utilized numerical modeling to predict airflow and resulting contaminant transport within an enclosure such as Chung and Dunn (1997), Longest et al. (2000), Hayashi et al. (2001), Lee et al. (2002), Zhao et al. (2002), Khan et al. (2004), Lin et al. (2004), Lee et al. (2005), Stamou and Katsiris (2005), and Tian et al. (2005). Lee et al. (2002) stressed the importance of using a realistic representation of the inlet velocity profile. Their study compared numerical results for two different velocity inlet profiles, one of which was a uniform profile across the inlet, the other of which was an experimentally measured inlet velocity profile, like a realistic profile. Their results concluded that the measured inlet profile produced better results both qualitatively and quantitatively, thus the assumed inlet velocity profile should be based on inlet duct configuration, Lee et al. (2002). The inlet configuration directly influences the velocity jet formed within the chamber, which affects the momentum flow and transport, Lee et al. (2005). For high Reynolds flow, the inertial effects of the flow can be the dominant factor in determining contaminant transport and decay; therefore, accurate representation of the velocity inlet boundary condition is necessary.
This study also examines turbulence effects on the airflow patterns within the chamber. For high Reynolds and Schmidt number flows, inertial transport mechanisms can overtake species transport via diffusion, leading to momentum dominated flow. Principally, for high enough Reynolds numbers, turbulent effects remain similar for different flow rates, leading to similar transport patterns, i.e. airflow characteristics within the chamber. This asserts that at a certain level of turbulence, airflow profiles of the turbulent jet will only change in their magnitude, not their shape.
This paper presents a comprehensive study of contaminant decay within a ventilated space using both CFD and experimental methodology. The main aspect of this study was an attempt to mirror contaminant concentration decay over time, as predicted by CFD, with the monitored contaminant decay. The effect of turbulence on species transport within the chamber was studied to determine the relative importance of diffusive transport. Furthermore, instead of using an experimentally measured inlet velocity profile as done in Lee et al. (2002), the upper ductwork for the chamber inlet was included in the computational domain and mesh, thus the resulting inlet velocity profile was directly computed via the inlet duct geometry. These goals were pursued because these parameters may influence the overall contaminant distribution and resulting decay. A clear understanding of these parameters is integral to accurately predict exposure and the duration of dangerous contaminant levels within a ventilated space. Comparison against experimental measurements was used for simulation validation.
In this study, experimental and computational methods were used to monitor and predict the transient concentration decay of a tracer gas within a ventilated space. Carbon dioxide was used as the tracer gas, as its presence is relatively easy to monitor over time and is often used to determine the air exchange rate of ventilated spaces. Furthermore, carbon dioxide is considered to be a common indoor contaminant, thus it is a prime candidate for a tracer species in this study.
The experimental data used was gathered in a specially built chamber located at the University of Iowa Hospitals and Clinics. The ventilation chamber, shown in Fig. 1, is ~3.40 m (L) × 3.44 m (W) × 2.29 m (H) and possesses one circular inlet, located at the center of the ceiling with a diameter of 0.673 m. An outlet is located at each corner of the room, giving the chamber four exhaust locations; a schematic is shown in Fig. 2. Air to the chamber was supplied by a dedicated ventilation system that automatically controlled flow rate, temperature, and relative humidity. Sensors in the main air duct leading to the chamber provide feedback to a system that controls the chamber airflow rate. Carbon dioxide was injected into the air duct and mixed with inlet air with a static mixing device prior to the ceiling inlet. Carbon dioxide was metered from a compressed air tank with a valve and calibrated rotameter. Initially, the air contained 330–370 p.p.m. of carbon dioxide, which was elevated to ~1000 p.p.m. before entering the chamber. The gas was fed into the chamber for at least 10 min in order to obtain a steady state, well-mixed initial condition before the beginning of the washout process.
Concentrations of carbon dioxide in the experimental chamber were measured simultaneously at four locations as shown in Figs 3 and and4.4. A Q-trak direct-reading instrument, TSI® Inc., was supported on poles in each location. The Q-trak measures carbon dioxide over a range of 5000 p.p.m. with an accuracy of ±3%, or 500 p.p.m., whichever is greater, and a resolution of 1 p.p.m. Instrument response time is 20 s. Each instrument was factory calibrated before use. After reaching steady state, the tank of carbon dioxide was turned off and the ensuing decay curve was measured by each Q-trak.
For comparison of the experimental data to numerical results, the time constant and ventilation performance of the measured results were quantified. This was achieved via the following methodology. Using the gathered experimental data, the time constant for the contaminant decay was computed by applying an exponential curve fit to the concentration values relative to time. To quantitatively fit the transient decay of the species, a material balance on a molar basis was applied.
The term on the left-hand side is the contaminant accumulation within the chamber; the terms on the right-hand side are the contaminant generation rate and contaminant removal rate. In this framework, V is the space volume, Q is the volumetric flow rate, C is the time-dependent concentration, G is the contaminant's generation rate, and KEXP is defined as a ventilation coefficient of performance and is represented as,
where τactual is the measured time constant, which is computed from the inverse of the slope obtained from a linear regression of the natural log of the transient contaminant concentration values. The term in the denominator of equation (2), τCMSTR, is calculated by dividing the chamber volume V by the flow rate Q, as seen in equation (3). The subtext ‘CMSTR’ denotes a completely mixed stir tank reactor, which assumes instantaneous mixing at all times. The initiation of the decay in carbon dioxide did not occur until after the duct work leading to the chamber was flushed with clean air, so the volume of that duct work was not considered when determining the value of V. Manipulation of equation (1), namely setting the species generation rate G to zero, rearranging the terms, and integrating yields the result of,
Taking the exponential of both sides of equation (4) yields the final result of,
where Co is the initial room concentration and Ct is the room concentration at any point in time. The contaminant decays in an exponential manner, with a characteristic time measured as,
which is the slope of equation (4). Again, KEXP is an experimental ventilation coefficient of performance that measures how close the experimental decay time is to the theoretical decay time. Using a spreadsheet, one can obtain an exponential curve fit with the same form as equation (4) and whose slope is the negative reciprocal of τactual. This was done for both the experimental and computational results.
Two different flow rates from the experimental study were examined, one at 509 m3 h−1 and one at 930 m3 h−1; these flow rates correspond to ~19 and 34 ACH. A typical office space may have an air change rate of about 6 ACH, Peyton and Abdou (1995), thus the flow rates under study are significantly higher. The space within the chamber did not contain any flow obstructions in order to simulate the airflow characteristics of an empty room. GAMBIT® 2.4, was used to create the computational domain that matched chamber dimensions. A tetrahedral mesh of 2167189 elements was created for CFD modeling as shown in Fig. 5.
FLUENT® 6.3 was used for the CFD modeling. The governing equations that were solved consisted of the continuity, momentum, and species transport equations. The k-ϵ turbulence model was used, mainly because it is commonly used as a standard turbulence modeling technique, thus results for comparison are widely available. The SIMPLE pressure–velocity segregated algorithm was applied as well, with a first-order difference scheme to approximate continuity, momentum, energy, and species equations.
The governing equations, known as the incompressible Navier–Stokes equations, were iteratively solved for an unsteady case in order to obtain a series of transient data sets.
where Ui is the velocity components in tensor notation and xi is the coordinate direction in tensor notation.
where P is the pressure and ν is the kinematic viscosity and νT is the turbulent kinematic viscosity, which is defined as follows,
Turbulence is modeled using the Reynolds-Averaged Navier–Stokes (RANS) model. In this model, two equations are solved. One of the equations is for k, the turbulent kinetic energy, which is written as:
where σk is the turbulent Prandlt number for k. The other equations solves for ε, the turbulent dissipation rate, which is represented as:
where σε is the turbulent Prandlt number for ε. The terms Cμ, Cε1, and Cε2 are constants. These constants and the turbulent Prandlt numbers σk, and σε are .
where μ and μT are the dynamic viscosity and turbulent dynamic viscosity, Sc and Sct are the Schmidt and turbulent Schmidt numbers, and YA is the species mass fraction.
The boundary condition used for all model inlets was a uniform velocity normal to the inlet surfaces. The outlet boundary conditions were set as pressure outlets to atmospheric pressure. For all the duct and chamber walls, the no-slip condition was applied, with zero diffusive flux of all species. In addition, the chamber was assumed isothermal.
The simulation produced fifteen data files that corresponded to a time sequence of species concentration decay curve. The convergence criterion for the problem was set to 10−6. After the simulation was completed, FLUENT® was used to export data for x-velocity, y-velocity, z-velocity, velocity magnitude, species mass fraction, and species concentration into a data file format readable by Tecplot® 360. Using Tecplot®, the contaminant concentrations were extracted at each of the four locations and recorded into a spreadsheet for analysis.
Table 1 presents a comparison of computational and experimental results, specifically for the species decay time constant and ventilation coefficient of performance at each location. The difference between KEXP (experimental coefficient of performance) and KCFD (computational coefficient of performance) was computed as a percent difference using KEXP as the standard. This essentially measures the difference between the predicted characteristic decay of the species versus the monitored species decay; the theoretical characteristic decay remains the same.
The results show that the k-ϵ CFD model agreed quite well with the experimental measurements. The largest error of the simulation was approximately −20.0%, which was found at location 4 for 19 ACH. Computational results were much better at other locations within the chamber, with error as low as −11.8% for location 2 at 19 ACH and −5.32% for location 2 at 34 ACH. Except for location 4 at 19 ACH, agreement between experimental and computational results resided in the range of −5.0 to −15.0%, which was acceptable considering the approximations employed by the k-ϵ model. However, it is important to note that in all cases, KCFD was lower than KEXP, meaning that the actual decay rate was less than the predicted decay rate, i.e. the CFD solution predicts a quicker washout of the contaminants than what was seen through experimentation. A possible reason for the CFD results under-predicting the decay time is the use of a RANS turbulence model, which does not resolve small turbulent eddy motion, but instead averages turbulent effects. In addition, the CFD solution assumed that the flow rates in the opposing inlet ducts were precisely the same, a condition that would have been difficult to exactly produce in the experimental chamber. Nonetheless, it is very encouraging that the model agreement was so accurate in the region located underneath the chamber inlet, denoted as Location 2. This is the area of highest turbulence, as well as the region most influenced by the incoming velocity jet. These results suggest that the computational model accurately resolved the shape of the inlet velocity profile and indicated that the numerical model provides a good representation of the actual flow fields and species decay within the chamber. A plot of the computational decay versus the measured experimental decay at location 2 for 19 ACH is shown in Fig. 6, with the y-axis as the natural log of the time-dependent concentration. The plot shows that the predicted decay is very close to the monitored experimental results with an error of only −11.8%.
It is interesting to note that the predicted KCFD for both flow rates was relatively the same. For all four locations at both 19 ACH and 34 ACH, the predicted clearance time ranged from 1.12 to 1.13 as a factor of the infinite mixing time. Thus, the computational time constant of the species decay was ~12–13% longer than computed assuming instantaneous mixing. This suggested that this particular room has a ventilation coefficient of performance of ~1.1 across a range of flow rates, namely for those >19 ACH and <34 ACH. This showed a limitation of the RANS turbulence model as KEXP varied from 1.2 to 1.4, suggesting that the ventilation coefficient of performance varies slightly with flow rate.
The result of a relatively constant KCFD for both flow rates was most likely a result of similar turbulent flow fields due to the high Reynolds number of the flow. For both flow rates, the inlet Reynolds number was quite large, ~18000 for 19 ACH and 32600 for 34 ACH. These Reynolds numbers were calculated using the average velocity at the inlet and the duct diameter. This result indicated that the air entering the chamber formed a highly turbulent velocity jet that had a large influence on contaminant transport below the inlet. Furthermore, the room Reynolds number was ~3500 for 19 ACH and 6400 for 34 ACH, which indicated that some level of turbulent flow persisted throughout the other regions of the room, namely those located away from the inlet. The rooms Reynolds numbers were computed using the average velocity in the room and the hydraulic diameter of the room in the xy-plane. This evidenced that turbulent/inertial transport may be the dominant mechanism for species transport within the chamber. Due to the turbulent eddy averaging technique of the k-ϵ model, the predicted turbulent effects within the chamber were relatively the same. This dominance of inertial effects is explored further in the next section.
Inertial flow was the dominant mechanism for fluid movement within the chamber due to the high Reynolds number, i.e. viscous effects are negligible. However, the importance of diffusive effects is not known, thus it is important to consider whether or not momentum/inertial transport is the main mechanism by which contaminant gases are exorcised from the chamber.
One method to examine the effects of turbulence or high Reynolds flow on the contaminant transport was to perform a general scaling analysis on the species transport equation. The generalized species transport equation is given by:
Assuming that the density of the species remains constant, ρ can be pulled outside of the derivative. The species diffusivity was a constant as well, allowing it to be pulled from inside the derivative. In addition, the solver utilized by FLUENT® for this problem was segregated, meaning that the flow velocities Uj were computed using the momentum equations, not the species equation, thus the velocities may be pulled from the derivatives as well. The species equation can now be simplified to:
Now the scaling parameters can be introduced as:
Substituting these into the species transport equation yields:
Now, the scaled equation is divided by the term ρDYmax/L2, yielding a result of:
The first term in brackets in front of the second term is recognized as the dimensionless Reynolds number, which typically describes the importance of inertial effects to viscous effects and is given by,
where Umax is the flow velocity, L is the characteristic length of the flow area being considered, μ is the kinematic viscosity, and ρ is the density. The second term in brackets is the Schmidt number, which describes the ratio of momentum and mass diffusivities, it is defined as,
Dividing through by the Reynolds and Schmidt numbers yield the equation:
where tc was chosen as:
This result is very similar to the result for characteristic time, τCMSTR, presented earlier in the experimental methods section. All the variables of the transport equation are on the order of one as a result of the scaling. The order of the Schmidt number is ~0.1–1 and the order of Reynolds number is ~104. This indicates that the species transport due to diffusive effects is negligible. The inertial term of the species transport equation has the largest effect on the predicted transient contaminant decay, which is implied by the computational results.
To examine the similarity of turbulent flows at various Reynolds number, Fig. 7 presents normalized velocity magnitude profiles of the turbulent jet for both 19 ACH and 34 ACH at 5, 10, and 15 min into the decay. The plots indicate that at 5 min, there are slight differences between the velocity profiles; however, at 10 and 15 min the curves collapse onto each other as the flow fields reach a quasi-steady state. The same behavior is seen in plots comparing the turbulent kinetic energy profiles. Fig. 8 shows the similarity between the turbulent jet behavior at the two different flow rates. Again, the profiles are markedly similar; indicating that for both flow rates inertial transport is the key transport mechanism that ventilates the chamber contaminants. This implies that the shape velocity profile of the turbulent jet will remain the same despite changes in velocity magnitude and hence Reynolds number. This is because the turbulent jet is a free shear flow; the turbulent flow behavior arises because of mean velocity differences formed away from the wall. One of the important features of the turbulent jet is self-similarity. The shape of the mean velocity of a round jet normalized by its centerline velocity remains unchanged and is independent of the Reynolds number, thus being self-similar. It is because the Reynolds stress in the turbulent jet is dominant over the viscous stress (Pope, 2003).
The effect of a varying geometrical inlet configuration was studied to determine what effect, if any, the inlet geometry has on the flow fields within the chamber. This analysis focused on the presence of the upper ventilation airways leading to the inlet and how their presence and shape influence airflow. For the previously presented computational results, a detailed inlet configuration was used that closely mirrored the actual chamber inlet geometry, as shown in Fig. 5. This geometry was somewhat time consuming to construct, thus it was desirable to determine whether or not that level of detail was necessary to obtain an accurate solution.
Two additional configurations were studied to determine whether or not the upper ventilation should be modeled and if so, how much geometrical accuracy was required. The first additional configuration did not possess any upper ductwork; the velocity inlet was a simple circular opening in the ceiling with a uniform velocity profile, identical to the previously presented Fig. 2. The second additional configuration had a simplified version of the upper ductwork attached to the inlet, which included two square opposing horizontal jets that converged to a vertical inlet to the room, shown in Fig. 9.
Examination of normalized velocity profiles for the three configurations exhibits vast differences between each, as shown in Fig. 10. Velocity profiles were normalized by extracting velocity magnitude line data across the domain of interest, and then dividing by the maximum velocity found, such that the normalized velocity would range from zero to one. Normalization of the x-coordinate position was done in a similar manner. The entering velocity jets for the two simplified configurations in Fig. 10a are much broader, while the detailed configuration produces a much thinner velocity jet. Differences are more recognizable in Fig. 10b, where the two additional configurations again produce broad velocity profiles, whereas the detailed configuration creates a much more diffuse, fork-like, velocity jet. Because the detailed inlet configuration closely resembled the physical dimensions of the actual ventilation chamber, it is reasonable to assume that it produces the most realistic approximation of the velocity profile. This is very important due to the high Reynolds number nature of the flow; transport via momentum flow would be highly sensitive to the inlet velocity profile, which directly influences the flow fields generated within the chamber. This conclusion is further supported through comparisons of the predicted ventilation coefficient of performance for the three different configurations to experimental data, presented in Table 2.
The results show that for the first additional configuration, the predicted KCFD for Locations 1, 3, and 4 was quite good with agreement as good as 1.04%; however, the result for Location 2 was highly inaccurate. In fact, the species decayed to ~5% of its starting value after the first minute of the washout process, thus it was not possible to determine an exponential slope for the decay at this location. This result was in no way similar to the monitored experimental results and evidenced that the first additional configuration decreased the accuracy of the solution. The second additional configuration, with its simplified inlet vent configuration produced an extremely inaccurate solution, which indicated that misrepresentation of the inlet geometry may cause severe errors in the computational solution. The same rapid decay below the inlet was encountered as well, which again was not physically accurate when examining the monitored experimental decay at that location. This indicated that excluding the upper ventilation or simplifying it reduced the accuracy of the solution. This decrease in solution accuracy is due to the unrealistic representation of the inlet velocity profile, which is a direct function of inlet configuration.
The detailed inlet configuration produced an interesting fork-like shaped profile that was not present in the other models. This shape is prominently seen in the normalized velocity magnitude plots in Fig. 10b. By rotating a slice of the plane used in Fig. 10b and adding a plot of velocity isosurfaces, this intriguing velocity profile is more clearly seen in Fig. 11. The fork-like shape was due to the airflow impingement of the opposing vents and inlet vent expansion. This impingement and expansion created a downwardly divergent flow structure. Without accurately modeling the upper ventilation geometry, this velocity inlet profile would not be resolved by the computational model. Thus, when trying to accurately model flow fields within a chamber or space, the configuration of the inlet and any upper ventilation ductwork may be very important, especially in turbulent flows.
Bends or curvatures present within upper ventilation ductwork can affect the flow field due to inertial effects. In highly turbulent flows, inertial effects can cause the flow to impinge upon a wall and in a sense ‘hug’ that portion of the duct, creating secondary flows. This causes asymmetric flow that results in an asymmetric velocity profile. With the use of two converging ducts as in this model, the velocity profiles are alternatively much more symmetric. Without modeling the upper ventilation ducts, one would have to enter x, y, and z components of the velocity inlet profile, as well as its magnitude. This is difficult to obtain analytically and time consuming to measure experimentally; however, by modeling the upper ventilation, the physical response of the flow to the duct geometry is included in the numerical solution. One drawback to this method is that it adds to the number of elements and nodes to the model, increasing computational requirements. However, in this instance it is advantageous to extend the model because the benefits in computational accuracy outweigh the added computational effort.
Some of the limitations present within this study are the small number of locations sampled within the room, as well as using only two different flow conditions, 19 ACH and 34 ACH. Due to the high Reynolds number nature of the flow, inertial effects overshadow other transport phenomenon, which limits the scope of conclusions that may be drawn concerning diffusive transport and natural convection. In addition, a RANS model was used, which is an approximate approach to resolve flow turbulence. Despite these limitations, the results support the use of CFD methods as a cost-effective means for studying occupational ventilation.
The results of this study show that the k-ϵ turbulence model can be successfully applied to the problem of room ventilation. The agreement with experimental results was good at both of the studied flow rates. The results suggest that for high Reynolds ventilating flows, the predicted ventilation coefficient of performance is relatively constant despite differing flow rates. It is expected that this result is due to the dominance of inertial transport mechanisms and the similarities of the turbulent profiles between two differing flow rates.
Furthermore, the predicted contaminant decay is highly sensitive to the incoming velocity profile at the inlet. It was found that accurate modeling of the ductwork above the inlet produced a more accurate solution as opposed to using a uniform velocity profile inlet condition. In addition, simplifications to the inlet duct geometry increased the error of the computation. In the end, it was found that accurate representation of the incoming velocity profile at the inlet was a key factor in modeling the flow fields and contaminant decay within the chamber. Experimental measurements of the inlet velocity profile can be used for modeling, as done in Lee et al. (2002), or the geometry of the upper ductwork can be included in the computational domain, as done in this study.
Key conclusions of this study are that for high Reynolds number ventilating flows, diffusive transport is small compared to inertial transport and in some cases may be neglected. Furthermore, the computational results are highly dependent on the entering velocity profile, which is a strong function of the inlet duct configuration. In addition, the predicted computational coefficient of performance indicated that the contaminant washout required ~12–13% more time than theoretically predicted. Moreover, model error as low as −5.32% indicates that CFD is an applicable tool to ventilation studies.
Environmental Health Sciences Research Center (NIH P30 ES05605).