J. agric. Engng Res. (1997) 66, 221 – 228
Comparison of Predicted and Measured Air Flow Patterns in a Mechanically Ventilated Livestock Building without Animals B. B. Harral; C. R. Boon Silsoe Research Institute, Wrest Park, Silsoe, Bedford MK45 4HS, UK (Receiy ed 5 September 1994; accepted in rey ised form 11 December 1996)
Air flow in a section from a full scale livestock building has been simulated using PHOENICS, a computational fluid dynamics program, to give the mean velocity and turbulence energy distribution under isothermal conditions. The section has sloping ceilings, an air intake at the apex, and side-wall extractor fans. A narrow inlet gap creates the horizontal air jet which is the basis of the ventilation system. A summer ventilation rate of 44 air changes per hour was simulated, with a Jet Momentum Number of 0?0076 and an inlet Reynolds Number of 5x104 based on a gap size of 0?1 m and a mean velocity of 6?1 m s21. Calculated mean air velocites within the section were in close agreement with values measured by a sonic anemometer. Despite its known shortcomings, use of a standard turbulence model to simulate the turbulence energy generation and dissipation gave good predictions of turbulence energy distribution. ÷ 1997 Silsoe Research Institute
1. Introduction Adequate air flow is essential in livestock buildings to ensure that air quality and thermal conditions are satisfactory for both animals and workers. Animals can occupy the building for 24 h each day whilst workers can be in the building for up to 8 h. Mean air speed and temperature are known to affect the thermal comfort of animals (Boon1). Mean air flow and turbulence affect the dispersion of gaseous and particulate pollutants which are harmful to the health and welfare of both animals and workers. It is important therefore to understand and quantify the air distribution patterns and turbulence characteristics of ventilating air flows, so that the limited amount of incoming fresh air can be distributed as efficiently as possible to create a healthy micro-climate around animals and humans. 0021-8634 / 97 / 030221 1 08 $25.00 / 0 / ag960140
Air movement in buildings is caused by ventilation services, metabolic heat, the external weather, and the movement of equipment, animals and people. The disturbance caused by moving objects is usually random and intermittent and its effect is assumed to be small compared with continuous processes. Mean air velocities and turbulence characteristics have been measured in experimental studies at model-scale (Hoff et al .2) and full-scale rooms (Zhang et al .3). These experiments involved rooms with flat ceilings and wall-mounted inlets. Only recently have measurements been made (Heber and Boon4) in a full scale enclosure with a sloping ceiling, a central ceiling inlet and side-wall exhaust fans, a common configuration in intensive livestock buildings in the UK. However, experimental studies are expensive and time consuming. Computational Fluid Dynamics (CFD), a numerical technique based on the solution of mean flow conservation equations coupled with suitable turbulence models, has been used in studies of room air movement (Awbi;5 Haghighat et al .;6 Lamers and van de Velde7) to predict the velocity and turbulence fields. However, these studies were either at model scale or had limited experimental data for comparison, and all involved simple rectangular geometries. The work reported here is the first stage in a study using CFD to predict the movement of gaseous and particulate pollutants in and around livestock buildings. Accurate prediction of the velocity and turbulence fields is a prerequisite. The paper compares the internal air velocities and turbulence energies predicted by a CFD simulation with measured values from full-scale experiments in a building section without livestock under isothermal conditions.
2. Experimental building and CFD grid A typical building used for intensive livestock production in the United Kingdom is shown in Fig. 1 .
221
÷ 1997 Silsoe Research Institute
222
B. B. HARRAL ; C. R. BOON
Livestock building ceiling
ceiling air intake
inlet
inlet
Air intake baffle
y z y x z
Experimental section
CFD grid Exhaust fan
Fig. 1. Experimental section of a typical liy estock building , with details of the air intake (inset) and a coarse CFD grid (superimposed)
Fresh air enters under the ridge, passes through openings on the ceiling centreline and is then directed to both sides of the enclosure by a horizontal baffle immediately below each ceiling opening (see inset in Fig. 1 ). The narrow gap between the baffle and the ceiling creates the jet that is the basis of the ventilation system. In the paper, the ceiling opening will be referred to as the air intake, while the air discharge gap between the baffle and the ceiling will be referred to as the inlet. Exhaust fans in both side-walls create the pressure differential necessary for the desired fresh air exchange rate. A full-scale section of a typical building (shown shaded in Fig. 1 ) has been built for experimental purposes, and the principles underlying its design were described by Carpenter et al .8 Main dimensions of the section are: length (x direction), 3?09 m; width (z direction), 11?73 m; side-wall height (y direction), 1?85 m; maximum distance from floor to ceiling (y direction), 3?34 m; intake length (x direction), 0?985 m; intake width (z direction), 0?4 m; inlet length (x direction), 0?985 m, and exhaust fan diameter, 0?45 m. The CFD analysis considers the flow domain, in this case the interior of the building section, to be divided into an array of cells, Fig. 1 . Since the section is assumed to be empty the central xy plane of the section is a plane of geometric symmetry. Since the
volumetric flow through both fans is assumed to be identical only half the section needs to be considered. The central yz plane is not a plane of symmetry because the outlet fans are not centrally mounted in the section side-walls. Thus, further reduction of the volume to be analysed is not possible. Figure 1 also shows a typical CFD grid superimposed on half the section. This is a non-orthogonal, body-fitted coordinate grid which enables the shape of the building and the geometry of the inlet and outlets to be followed closely.
3. CFD Analysis The Navier-Stokes equations describe the turbulent flow of a compressible Newtonian fluid. These equations can be found in fluid mechanics texts, for example Raudkivi and Callander,9 and will not be repeated here. At the velocities encountered in the building section, compressibility effects can be neglected. The CFD technique solves numerically the Reynolds-averaged form of these equations within each cell in the domain. The Reynolds-averaging process considers the instantaneous fluid velocity to be the sum of a mean and a fluctuating component, the turbulence. In general, the high-frequency and small scale fluctuations of turbulent flow cannot be quan-
223
COMPARISON OF PREDICTED AND MEASURED AIR FLOW PATTERNS
tified directly and modelling is used. Turbulence modelling relates some or all of the turbulent velocity fluctuations to the mean flow quantities and their gradients, in a physically plausible manner. The best compromise at the present time between acceptable computing cost and accuracy in the simulation of turbulent flow phenomena is the k 2 » turbulence model (Launder and Spalding10). This model consists of two additional differential equations for the specific turbulence kinetic energy (k ) and the turbulence energy dissipation rate (» ), and includes a turbulent viscosity analogous to laminar viscosity. The model considers the turbulent viscosity to be isotropic. The k 2 » model is valid only for fully turbulent flows. Close to solid walls, there are boundary layer regions where the local Reynolds Number is so small that viscous effects predominate over turbulent effects. To account for this effect and for the large gradients of variables near the wall, the log-law function method of Launder and Spalding10 has been used. A relatively high summer ventilation rate of 44 air changes per hour was modelled. Bulk jet velocity was 6?1 m s21 through an inlet gap of 0?1 m, total ventilation rate was 1?2 m3 s21 and Jet Momentum Number (Barber et al .11) was 0?0076. The inlet Reynolds Number was 5x104 based on a gap of 0?1 m and bulk jet velocity of 6?1 m s21, indicating fully turbulent flow. The turbulence kinetic energy and dissipation rate at the intake were given values based on recommendations by Hjertager and Magnussen.12 The air intake was defined as a fixed pressure boundary and the pressure set equal to zero. The plane of the exhaust fan was defined as a fixed flux boundary and the mass flow set equal to the product of ventilation rate and density. These boundary conditions were chosen to reproduce the experimental conditions as closely as possible. The CFD code PHOENICS (Rosten and Spalding13) was used to create the cell geometry and solve the governing equations. The solution provided by the discrete form of the flow equations approaches the continuous form solution as the number of cells is increased, particularly where pressure and velocity gradients are high. To increase resolution and check for grid dependence, three grids were created containing 2496, 21216 and 84000 cells. The finest grid had a greater proportion of cells concentrated in the inlet zone, the smallest being 42 3 14 3 33 mm in real size. PHOENICS solves discretised forms of the fluid flow equations in which the differentials are replaced by numerical approximations. To cope with the non-linear coupling of variables in the equations, PHOENICS uses an iterative scheme and underrelaxation to reach a solution. Convergence was
assumed when total residual errors fell to approximately 1% of average fluxes.
4. Experimental measurements Air velocities were measured on three orthogonal planes, shown in Fig. 2 , at 180 points in total. At points close to the inlet, air speed was measured with a hot-wire anemometer which gave a time-averaged mean value. Elsewhere, a sonic anemometer was employed which gave sampled time-histories of the three orthogonal velocity components. These velocities were the space-averaged values over the sonic path length, in this instance 150 mm. Mean and variance values were calculated from the timehistories during post-processing. A specific turbulence kinetic energy (k ) was also calculated by assuming that the variance of the measured velocity components represented turbulence alone rather than part of a time-varying flow. k 5 1–2 (s 2x 1 s 2y 1 s 2z)
(1)
where s x , s y , s z are the standard deviation of velocity in the x, y, z directions, respectively. The quoted accuracy of the sonic anemometer, measuring the mean velocity in a free stream, is 1?5%. At the edge of the jet where velocity gradients and turbulence frequencies are high, the accuracy is dependent on two factors. First, it gives volumeaveraged data rather than point data owing to its size (path length 150 mm), the effect on accuracy will depend on the shape of the velocity profile across the device. Secondly, it can underestimate the magnitude of velocity fluctuations because its frequency response Air intake
y x z
Exhaust fan
Fig. 2. Velocity measurement points ( j) on xy , xz and yz planes (crossed-hatched)
224
B. B. HARRAL ; C. R. BOON
is not only a function of the mean value, but is also limited to 28 Hz because of sampling frequency.
5. Results and discussion Incoming air does not travel directly to the exhaust fan but circulates in primary and secondary recirculation patterns before finally leaving through the fan. Primary and secondary patterns relate to the entire three-dimensional simulated space, not the smaller zones occupied by animals or workers. Predicted velocity vectors in Fig. 3 show the predicted primary circulation in the yz plane on the centreline of the inlet. Air speed and flow direction are given by arrow size and arrow direction respectively. Primary circulation is driven by the viscous drag forces between the high speed incoming jet and the main body of air within the building. Secondary circulation in xy and xz planes is a consequence of the primary circulation and the offset position of the exhaust fan. Figure 4 shows the predicted secondary circulation in a horizontal xz plane through the centre of the outlet. For
clarity, a regular array of vectors is plotted in Figs 3 and 4 , obtained by interpolating the velocities from the simulation using the finest grid. A complex structure of recirculation zones is evident. Figure 5 compares the predicted and measured x, y and z direction velocity components on a vertical line (shown dotted in Fig. 3 ) 0?325 m from the inlet. As the number of cells is increased the results indicate that little improvement is gained above 21216 cells. A similar trend is shown at other positions on the measurement planes. It is only in the jet, at 3?3 m from the floor, that the z direction velocity is significantly changed by increasing the number of cells from 21216 to 84000, and would probably continue to rise towards the measured value as the grid was further refined. However, the results from the simulation using the finest grid can be assumed to be grid-independent throughout most of the domain. Figure 5 appears to show the rather surprising result that the closest agreement between predicted and measured values is obtained using only 2496 grid cells. Since the predicted values in this case are not gridindependent, this must be regarded as coincidental.
Air intake
Baffle
y Reference velocity 3.0 m s–1 z
Fig. 3. Predicted y elocity y ectors on the central y ertical yz plane showing the primary air circulation pattern; air speed and flow direction are giy en by arrow size and arrow direction respectiy ely
225
COMPARISON OF PREDICTED AND MEASURED AIR FLOW PATTERNS
Outlet
x Reference velocity 3.0 m s–1 z
Fig. 4. Predicted y elocity y ectors on a horizontal plane through the outlet duct , height (y direction) , 0?95 m , showing the secondary air circulation patterns; air speed and flow direction are giy en by arrow size and arrow direction respectiy ely
Although grid-independent solutions have been obtained, the predicted mean velocities do not converge on the measured values. This may be due to physical inconsistencies between the simulation and reality, and to numerical approximation and simplifications in the models used in the simulation. For example, a small difference in alignment between the real and simulated inlet, which in practice would not affect the overall operation of the ventilation system but is too small to be included in the grids, could affect the jet direction and modify recirculation patterns locally. Sub-grid scale features such as instrumentation supports and traversing mechanisms will also affect flow locally. Figure 6 shows the variation in the difference between measured and predicted velocities in the three cartesian directions and over the three measurement planes (shown in Fig. 2 ). Percentage difference was calculated from 2 measured velocity Spredicted velocity D 3 100 bulk jet velocity
Differences have been expressed as a percentage of the bulk jet velocity rather than the respective measured value to avoid the large percentages that
can arise when velocities are insignificantly small. Predicted velocities were taken from the simulation using the maximum number of cells. The magnitude of the difference is less than 5% in most areas. This is good considering the complex nature of the flow and the wide range of velocities, illustrated in Figs 3 and 4 . Some larger difference values do occur, notably, in the z direction, in the jet and at opposite corners of the xy plane. The maximum percentage differences between the predicted and measured velocities for the three grids, both at the inlet and in regions away from the inlet, are given in Table 1. The finest grid under-predicts the velocity at the inlet by 14.8% based on a bulk jet velocity of 6?1 m s21. This value is comparable to the 14?6% difference found by Lamers and van de Velde7 near the inlet of a model-scale enclosure with a horizontal ceiling and inlet bulk jet velocity of 15 m s21. The 20% error in the z-direction velocity, given in Table 1, occurs in a recirculation zone at the edge of the incoming jet. It remains the same irrespective of the number of cells and is due to local recirculation resulting in the predicted velocity being in the opposite direction to the measured value. It may be a feature of the modelling technique, in particular the
226
B. B. HARRAL ; C. R. BOON
3
Distance from floor, m
Distance from floor, m
3
2
1
–1
0
1
2
1
–1
0
Velocity (x direction), m s–1
1
Velocity (y direction), m s–1
Distance from floor, m
3
2
1
–1
0
1
2
3
4
5
6
7
8
–1
Velocity (z direction), m s
Fig. 5. Measured and predicted y elocities in the x , y and z directions , on a y ertical line 0?325 m from the inlet (shown dotted in Fig. 3); j measured y alues; numerical grid: ? ? ? ? ? 2496 cells; – – – 21216 cells; —— 84000 cells
k 2 » turbulence model which is known to overpredict k in regions of high shear and hence change the flow pattern locally. Figure 7 shows the measured and predicted specific turbulence kinetic energy distributions on the central yz plane. The comparison is useful for two reasons. First, it indicates the suitability of the turbulence model for the calculation of the velocity fields. Secondly, it gives an indication of the accuracy of the predicted turbulence energy, which is used by random walk models in the prediction of dust and gas dispersion. Near the inlet, contours of measured turbulence, that is turbulence energy derived from air velocity
measurements, have been restricted to the zone covered by measuring points because extrapolation in this highly turbulent region is not justifiable. The distributions of measured and predicted turbulence energy are similar, with high values near the ceiling where high shear strain rates arise from large velocity gradients across the jet, and a minimum in the centre where velocity gradients are low. Measured energy in the jet is likely to be underestimated because of the limitations of the instrumentation outlined earlier, while the predicted energy (k ) may be too high because of factors inherent in the equation describing the production of k .
COMPARISON OF PREDICTED AND MEASURED AIR FLOW PATTERNS
x direction
y direction
227
z direction
xy plane
yz plane
xz plane
Fig. 6. Contours of the difference between measured and predicted y elocities , expressed as a percentage of the bulk jet y elocity , 6?1 m s 21
Close to the floor, tubulence energy is underestimated by a factor of about two. It is unlikely that the sonic anemometer itself is contributing significantly to the turbulence. A probable cause is that the standard k 2 » turbulence model is inadequate because it assumes an isotropic turbulent viscosity. This possibility can be tested by comparing measured and calculated velocity variances. In some areas, notably near the floor, the measured velocity variances for the three orthogonal directions differ from the mean by up to a factor of two. If velocity variances are calculated from predicted velocities and turbulent viscosities, then the variation with direction in the calculated variances near the floor is found to be
Inlet Baffle
Inlet Baffle
Table 1 Maximum difference between predicted and measured velocities as a percentage of the bulk jet velocity (6?1 m s21) At inlet , At measurement points other than % inlet , % ————– ———————————————– Grid cells z direction x direction y direction z direction 2496 21 216 84 000
244?6 232?7 214?8
8?0 9?0 7?2
11?2 14?8 15?4
220?4 222?8 221?8
Fig. 7. Measured (top) and predicted (bottom) specific turbulence kinetic energy contours (J kg 21) on the central y ertical yz plane ( d measurement point)
228
B. B. HARRAL ; C. R. BOON
small, indicating that the assumption of isotropic turbulent viscosity is invalid in this area. A turbulence model such as the Reynolds-stress-transport model that solves six additional equations for turbulent velocity fluctuations, rather than two, may be more appropriate but computational cost will be greater.
References 1
2
3
6. Conclusions It has been shown that the isothermal ventilation flow in a livestock building can be represented by the results of CFD calculations. A relatively coarse grid and the isotropic viscosity k 2 » turbulence model, give results which are quantitatively in good agreement with measured values. Refining the grid improves the agreement in the jet. Once grid effects are overcome the determining factor for the accuracy of simulated fluid flow phenomena is the validity of mathematical models, such as the turbulence model and log-law wall function model, and their constants. Greater accuracy may be obtained by tuning the constants appearing in the models to suit the particular problem or by resorting to a more complex turbulence model. The results from using the isotropic viscosity k 2 » turbulence model in the livestock building analysis indicate that the model is adequate when calculating mean air velocity, despite anisotropy in the measured turbulence field. These factors will need to be examined when analysing more complex situations that include obstacles in the building, metabolic heat input, and dust generation and transport. The full-scale experimental facility will enable CFD predictions to be verified against measured results as more parameters are included in the simulation.
Acknowledgements This research is funded by the Pigs, Eggs and Poultry Division of the Livestock Products Policy Group, Ministry of Agriculture, Fisheries and Food.
4
5
6
7
8
9
10
11
12
13
Boon C R The effect of air speed changes on the group postural behaviour of pigs. Journal of Agricultural Engineering Research 1982, 27(1), 71 – 79 Hoff S J; Janni K A; Jacobson L D Three-dimensional buoyant flows in a scaled model, slot ventilated, livestock confinement facility. Transactions of the American Society of Agricultural Engineers 1992, 35(2), 671 – 686 Zhang J S; Christianson L L; Wu G J; Riskowski G L Detailed measurements of room air distribution for evaluating numerical simulation models. Transactions of the American Society of Heating, Refrigerating and Air-Conditioning Engineers 1992, 98, Part 1 Heber A J; Boon C R Air velocity characteristics in an experimental livestock building with nonisothermal jet ventilation. Transactions of the American Society of Heating, Refrigerating and Air-Conditioning Engineers 1993, 99, 1139 – 1151 Awbi H B Application of computational fluid dynamics in room ventilation. Building and Environment 1989, 24(9), 73 – 84 Haghighat F; Wang J C Y; Jiang Z Three-dimensional analysis of airflow pattern and contaminant dispersion in a ventilated two-zone enclosure. Transactions of the American Society of Heating, Refrigerating and AirConditioning Engineers 1990, 96(1), 831 – 839 Lamers A P G G; van de Velde R Air flow patterns in ventilated rooms. The PHOENICS Journal of Computational Fluid Dynamics and its applications 1989, 2(2), 219 – 238 Carpenter G A; Moulsley L J; Randall J M Ventilation investigations using a section of a livestock building and air flow visualisation by bubbles. Journal of Agricultural Engineering Research 1972, 17(4), 323 – 331 Raudkivi A J; Callander R A An Introduction to Advanced Fluid Mechanics. London: Edward Arnold, 1975 Launder B E; Spalding D B The numerical computation of turbulent flows. Computer Methods in Applied Mechanics and Engineering 1974, 3, 269 – 289 Barber E M; Sokhansanj S; Lampman W P; Ogilvie J R Stability of air flow patterns in ventilated spaces. ASAE Paper No. 82 – 4551, American Society of Agricultural Engineers, St. Joseph, USA, 1982 Hjertager B H; Magnussen B F Calculation of turbulent three-dimensional jet-induced flow in rectangular enclosures. Computers and Fluids 1981, 9(4), 395 – 407 Rosten H I; Spalding D B The PHOENICS beginners guide, CHAM TR / 100, Concentration, Heat and Momentum Limited, London, UK, 1987