Airﬂow Characteristics Downwind a Naturally Ventilated Pig Building with a Roofed Outdoor Exercise Yard and Implications on Pollutant Distribution

: The application of naturally ventilated pig buildings (NVPBs) with outdoor exercise yards is on the rise mainly due to animal welfare considerations, while the issue of emissions from the buildings to the surrounding environment is important. Since air pollutants are mainly transported by airﬂow, the knowledge on the airﬂow characteristics downwind the building is required. The objective of this research was to investigate airﬂow properties downwind of a NVPB with a roofed outdoor exercise yard for roof slopes of 5 ◦ , 15 ◦ , and 25 ◦ . Air velocities downwind a 1:50 scaled NVPB model were measured using a Laser Doppler Anemometer in a large boundary layer wind tunnel. A region with reduced mean air velocities was found along the downwind side of the building with a distance up to 0.5 m (i.e., 3.8 times building height), in which the emission concentration might be high. Additional air pollutant treatment technologies applied in this region might contribute to emission mitigation e ﬀ ectively. Furthermore, a wake zone with air recirculation was observed in this area. A smaller roof slope (i.e., 5 ◦ slope) resulted in a higher and shorter wake zone and thus a shorter air pollutant dispersion distance.


Introduction
Natural ventilation systems have the advantages of low capital investment, energy saving [1], and absence of noise [2] compared with mechanical ventilation systems. Thus, they are widely implemented in livestock buildings. Naturally ventilated pig production systems equipped with outdoor exercise yards are receiving increasing interest in Europe, the United States, building was limited. Moreover, the configurations of NVPBs with roofed outdoor exercise yards are very different from street canyons and dairy houses. There is lack of airflow information about this type of buildings with roofed outdoor exercise areas. Therefore, it is required to understand in detail airflow characteristics including both mean (time-averaged) velocity and turbulent fluctuations downwind NVPBs with outdoor access. This knowledge can contribute to a better understanding of the transport and dispersion processes of airborne pollutants and of an optimised roof design.
Wind tunnel and scaled model experiments are widely applied in aerodynamics studies because of their advantages of allowing fully controlled boundary conditions, working with real airflow [29], providing a large amount of data in a short time [30], and the flexibility in the experimental setup [11]. Because of the scaled-down model, similarity criteria e.g., geometric similarity, boundary similarity, Reynolds number similarity have to be met in order to make the wind tunnel experimental results be comparable to the results from full-scaled buildings [31]. By carefully checking the similarity criteria, the wind tunnel tests method was therefore used in this paper in order to conduct the airflow characteristics research. Moreover, the data obtained from the wind tunnel measurements can also be used to validate accuracy and reliability of computational fluid dynamics (CFD) models and therefore to investigate more complex flow (e.g., airflows above emitting surfaces) and dispersion phenomena, to perform comprehensive parametric studies, and finally to mitigate the production of emissions from livestock buildings.
The objectives of this study were to investigate mean and turbulent characteristics of airflows downwind of a NVPB with an outdoor exercise yard covered with roofs with different roof slopes, to predict air pollutant distribution and dispersion properties, and to provide valuable experimental data for CFD validation. The novelties of this study are as follows: • It is the first to provide detailed airflow information around a naturally ventilated livestock building combined with a roofed outdoor area; • It predicts the potential distribution of gaseous pollutants from the building using airflow measurement results; • It provides a large amount of experimental data of both mean velocity and turbulent fluctuations downwind the building with a high resolution that can be used for CFD validation.

Scaled Pig Building Model
A scaled model of a naturally ventilated pig building with an outdoor exercise yard was used in this study. The prototype pig barn, situated in the state of Lower Saxony in northwest Germany, was designed for rearing around 80 fattening pigs. The scaled model was a 1:50 geometric reduction of the full-scale pig barn and was constructed at the Leibniz Institute for Agricultural Engineering and Bioeconomy (ATB), Germany.
The scaled model was made of transparent acrylic glass and consisted of an indoor housing area and an outdoor exercise area. The external dimension of the model was 0.427 m (length) × 0.256 m (width) × 0.130 m (height). The housing area had two sidewall openings with opening heights of 0.020 m and 0.064 m, respectively. Eight pigpens with open pen fronts were placed in the housing area. The height and width of all pigpens were 0.020 m and 0.054 m, respectively. The length of the two pigpens located next to gable walls was 0.056 m, and of other pigpens was 0.050 m. Pigpen walls had a thickness of 0.002 m, and all other walls of the building model had a thickness of 0.003 m. The roof of the building housing area had a fixed slope of 15 • . The free access between indoor and outdoor areas through plastic strips or rotating doors in the prototype pig barn was constructed by a 0.025 m high acrylic sheet in the scaled model. The outdoor exercise area had a flexible roof with a length of 0.108 m and with a fixed top part, which totally covered the outdoor yard. There was a 0.010 m high sidewall but were no gable walls in the outdoor area. The prototype pig building was constructed to direct the excretion behaviour of pigs to the outdoor exercise yard [32], in which solid floors in the housing area and slatted floors with a deep pit in the outdoor area were adopted. Therefore, in this study we only considered the roof slope variations for the outdoor area, where the majority of emissions are expected to come from. Three cases with roof slopes of 5 • , 15 • , and 25 • were studied in this paper. It has been found that the presence of animals has an insignificant impact on the airflow and pollutant dispersion within the barn [33], therefore, pigs and other internal structures that have smaller dimensions than the group of pigs, such as feeders, drinkers, metal bars, and slatted floors are expected to have minor influences on the airflow field downwind of the building. To simplify the model construction, pigs and these internal structures were not constructed. Detailed dimensions of the scaled pig building model and the three roof slope variations are illustrated in Figure 1.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 4 of 19 0.025 m high acrylic sheet in the scaled model. The outdoor exercise area had a flexible roof with a length of 0.108 m and with a fixed top part, which totally covered the outdoor yard. There was a 0.010 m high sidewall but were no gable walls in the outdoor area. The prototype pig building was constructed to direct the excretion behaviour of pigs to the outdoor exercise yard [32], in which solid floors in the housing area and slatted floors with a deep pit in the outdoor area were adopted. Therefore, in this study we only considered the roof slope variations for the outdoor area, where the majority of emissions are expected to come from. Three cases with roof slopes of 5°, 15°, and 25° were studied in this paper. It has been found that the presence of animals has an insignificant impact on the airflow and pollutant dispersion within the barn [33], therefore, pigs and other internal structures that have smaller dimensions than the group of pigs, such as feeders, drinkers, metal bars, and slatted floors are expected to have minor influences on the airflow field downwind of the building. To simplify the model construction, pigs and these internal structures were not constructed. Detailed dimensions of the scaled pig building model and the three roof slope variations are illustrated in Figure 1.

Boundary Layer Wind Tunnel and Measurement Devices
The experiments were carried out in a large boundary layer wind tunnel (BLWT) at ATB, Germany. The BLWT was specially designed to investigate ventilation and dispersion processes in agricultural buildings [34][35][36][37], and has also been used to generate datasets for CFD validation [24,38,39]. The wind tunnel is 28.5 m long, consisting of an air inlet fitted of honey combs, an air outlet equipped with an axial fan, and a 19.5 m-long test section. The cross-sectional area of the test section is 3 m (width) × 2.3 m (height). A combination of six spires and roughness elements was used to create a boundary layer flow. The spires were installed at the entrance of the test section. The roughness elements consisting of two sizes of right-angled steel brackets, with dimensions length × width × height of 0.010 m × 0.004 m × 0.004 m and 0.004 m × 0.002 m × 0.002 m for small and big brackets respectively, were arrayed in staggered rows downstream the spires. The total length of the roughness elements was 9.6 m with 0.2 m space between each row. The 1:50 scaled pig building model was placed at the symmetry line of the wind tunnel at 1.2 m downstream from the roughness elements. The model was oriented with sidewall openings perpendicular to the approaching flow and the outdoor yard at the downwind side. The blockage ratio of the scaled model to the cross-section of the wind tunnel was 0.8%, which is far less than the recommended maximum value of 5% for wind tunnel tests in VDI-guideline 3783/12 [40], and thus the tunnel effect can be neglected. Figure 2 shows photographs of the wind tunnel with the scaled model placed inside.

Boundary Layer Wind Tunnel and Measurement Devices
The experiments were carried out in a large boundary layer wind tunnel (BLWT) at ATB, Germany. The BLWT was specially designed to investigate ventilation and dispersion processes in agricultural buildings [34][35][36][37], and has also been used to generate datasets for CFD validation [24,38,39]. The wind tunnel is 28.5 m long, consisting of an air inlet fitted of honey combs, an air outlet equipped with an axial fan, and a 19.5 m-long test section. The cross-sectional area of the test section is 3 m (width) × 2.3 m (height). A combination of six spires and roughness elements was used to create a boundary layer flow. The spires were installed at the entrance of the test section. The roughness elements consisting of two sizes of right-angled steel brackets, with dimensions length × width × height of 0.010 m × 0.004 m × 0.004 m and 0.004 m × 0.002 m × 0.002 m for small and big brackets respectively, were arrayed in staggered rows downstream the spires. The total length of the roughness elements was 9.6 m with 0.2 m space between each row. The 1:50 scaled pig building model was placed at the symmetry line of the wind tunnel at 1.2 m downstream from the roughness elements. The model was oriented with sidewall openings perpendicular to the approaching flow and the outdoor yard at the downwind side. The blockage ratio of the scaled model to the cross-section of the wind tunnel was 0.8%, which is far less than the recommended maximum value of 5% for wind tunnel tests in VDI-guideline 3783/12 [40], and thus the tunnel effect can be neglected. Figure 2 shows photographs of the wind tunnel with the scaled model placed inside.
The free stream wind speed at the wind tunnel inlet (U inlet ) was measured using a Prandtl tube, connected to a pressure transducer MKS Baratron ® Type 120A (MKS Instruments, Andover, USA). The Prandtl tube was located at the centre of the entrance of the test section at a height of 1.3 m from the wind tunnel floor. Air velocity and turbulence around the scaled model were measured using a 2D fibre-optic Laser Doppler Anemometer (LDA) (Dantec Dynamics, Skovlunde, Denmark) combined with the BSA Flow Software package (Dantec Dynamics, Skovlunde, Denmark). The LDA probe head was 0.06 m in diameter and 0.45 m in length and provided a focal length of 0.25 m. The LDA probe was mounted on a three-dimensional computer-controlled traverse system that allowed automated and precise probe positioning with an uncertainty of <0.1 mm. A fog generator Tour Hazer II (Smoke Factory, Burgwedel, Germany) was placed at the wind tunnel inlet to produce seeding particles for LDA measurements. The temperature, relative humidity, and static pressure of the ambient air were measured using a FHAD 46x sensor with ALMEMO ® D6 plug (AHLBORN Mess-und Regelungstechnik GmbH, Ilmenau, Germany). The data were sampled approximately every hour.
The mean values of temperature, relative humidity, and static pressure for each day were used to calculate the air density, and therefore to calculate the U inlet by using together with the data from the Prandtl tube. The free stream wind speed at the wind tunnel inlet (Uinlet) was measured using a Prandtl tube, connected to a pressure transducer MKS Baratron ® Type 120A (MKS Instruments, Andover, USA). The Prandtl tube was located at the centre of the entrance of the test section at a height of 1.3 m from the wind tunnel floor. Air velocity and turbulence around the scaled model were measured using a 2D fibre-optic Laser Doppler Anemometer (LDA) (Dantec Dynamics, Skovlunde, Denmark) combined with the BSA Flow Software package (Dantec Dynamics, Skovlunde, Denmark). The LDA probe head was 0.06 m in diameter and 0.45 m in length and provided a focal length of 0.25 m. The LDA probe was mounted on a three-dimensional computer-controlled traverse system that allowed automated and precise probe positioning with an uncertainty of <0.1 mm. A fog generator Tour Hazer II (Smoke Factory, Burgwedel, Germany) was placed at the wind tunnel inlet to produce seeding particles for LDA measurements. The temperature, relative humidity, and static pressure of the ambient air were measured using a FHAD 46x sensor with ALMEMO ® D6 plug (AHLBORN Mess-und Regelungstechnik GmbH, Ilmenau, Germany). The data were sampled approximately every hour. The mean values of temperature, relative humidity, and static pressure for each day were used to calculate the air density, and therefore to calculate the Uinlet by using together with the data from the Prandtl tube.

Reynolds Number Independence Study
In order to obtain a fully developed turbulent flow in the wind tunnel, a Reynolds number independence study for the approaching flow was carried out. The wind profile was measured at 20 positions along a vertical line ranging from 0.003 m to 0.6 m above the wind tunnel floor at the upstream edge of the scaled model (i.e., at Line 1 in Figure 3). The wind profile measurements were performed without the presence of the scaled model at free stream wind speeds (Uinlet) of 4, 6, 8, 10, and 12 m·s −1 , respectively. The streamwise (U) and vertical (W) velocity components were measured using the 2D LDA. The measurements at each measurement position were taken continuously until the sampling number reached 40,000 readings or the maximum sampling time reached 820 s before

Reynolds Number Independence Study
In order to obtain a fully developed turbulent flow in the wind tunnel, a Reynolds number independence study for the approaching flow was carried out. The wind profile was measured at 20 positions along a vertical line ranging from 0.003 m to 0.6 m above the wind tunnel floor at the upstream edge of the scaled model (i.e., at Line 1 in Figure 3). The wind profile measurements were performed without the presence of the scaled model at free stream wind speeds (U inlet ) of 4, 6, 8, 10, and 12 m·s −1 , respectively. The streamwise (U) and vertical (W) velocity components were measured using the 2D LDA. The measurements at each measurement position were taken continuously until the sampling number reached 40,000 readings or the maximum sampling time reached 820 s before moving the LDA probe to another position. The average sampling rate during experiments was around 300 s −1 . A 10s pause between each measurement position was set in order to minimise the disturbance of the movement of the LDA probe to the flow field. The above LDA setup was chosen according to a preliminary experiment for the reproducibility of statistic results, in which air velocities at heights of 0.01 m, 0.03 m, 0.06 m, 0.1 m, and 0.2 m along Line 1 were measured and repeated three times at U inlet of 8 m·s −1 . The results showed that with this setup the measurement uncertainty for the time-averaged U and W velocities was 0.3%, and 5.5%, respectively ( Figure 4). The above mentioned LDA setup was used for all air velocity measurements in this study.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 6 of 19 uncertainty for the time-averaged U and W velocities was 0.3%, and 5.5%, respectively ( Figure 4). The above mentioned LDA setup was used for all air velocity measurements in this study.

Stability Study
To establish a stable boundary layer airflow that represents a farmland terrain, the characteristics of wind profiles within the region of interest were assessed. The arrangement of the roughness elements was adjusted accordingly until reaching a desirable boundary layer flow following the VDI-guideline 3783/12 [40]. Air velocities at two vertical lines, one at the upstream edge of the scaled model (i.e., Line 1 in Figure 3) and the other at the downstream end of the region of interest (i.e., Line 2 in Figure 3), from 0.003 m to 0.6 m height along the symmetry line of the model uncertainty for the time-averaged U and W velocities was 0.3%, and 5.5%, respectively ( Figure 4). The above mentioned LDA setup was used for all air velocity measurements in this study.

Stability Study
To establish a stable boundary layer airflow that represents a farmland terrain, the characteristics of wind profiles within the region of interest were assessed. The arrangement of the roughness elements was adjusted accordingly until reaching a desirable boundary layer flow following the VDI-guideline 3783/12 [40]. Air velocities at two vertical lines, one at the upstream edge of the scaled model (i.e., Line 1 in Figure 3) and the other at the downstream end of the region of interest (i.e., Line 2 in Figure 3), from 0.003 m to 0.6 m height along the symmetry line of the model

Stability Study
To establish a stable boundary layer airflow that represents a farmland terrain, the characteristics of wind profiles within the region of interest were assessed. The arrangement of the roughness elements was adjusted accordingly until reaching a desirable boundary layer flow following the VDI-guideline 3783/12 [40]. Air velocities at two vertical lines, one at the upstream edge of the scaled model (i.e., Line 1 in Figure 3) and the other at the downstream end of the region of interest (i.e., Line 2 in Figure 3), from 0.003 m to 0.6 m height along the symmetry line of the model were measured. Measurements were carried out using the 2D LDA without the presence of the scaled model. Free stream wind speed in the wind tunnel was 8 m·s −1 . Parameters that define the roughness class of the boundary layer, e.g., the wind profile power exponent, roughness length, and turbulence intensity were examined.

Measurement of Airflows Downwind the Building
In order to investigate airflow characteristics and to predict air pollutant dispersion properties downwind the building, air velocities at three vertical planes (Planes 1-3, Figure 3 In this study, X, Y and Z denoted the axes of the coordinate system, which were aligned with the streamwise, spanwise and vertical wind directions, respectively. The origin of the coordinate system was located at the centre of the downstream edge of the scaled model at the wind tunnel floor. The airflow measurement positions and the definition of the coordinate system are depicted in Figure 3. According to the wind profile measurement results (described in Section 3.1), a stable and fully developed turbulent flow that represents a farmland terrain boundary layer could be obtained when the free stream wind tunnel speed (U inlet ) was at least 8 m·s −1 . Therefore, all airflow measurements downwind the scaled model were performed at U inlet of 8 m·s −1 .

Data Analysis
The mean and the standard deviation of the air velocity in streamwise and vertical directions at each measurement position were provided by BSA Flow Software (Dantec Dynamics, Skovlunde, Denmark). Air velocity and air turbulence characteristics were processed from 40,000 samples at each position to ensure statistically reliable results.
The velocity magnitude calculated from two-dimensional velocity components was defined as: where V 2D is the velocity magnitude, m·s −1 ; U and W are mean air velocities in streamwise and vertical directions, respectively, m·s −1 . Turbulence intensity and turbulent kinetic energy of the airflow were calculated by Equations (2) and (3), respectively.
where TI is the turbulence intensity, %; TKE is the turbulent kinetic energy, m 2 ·s −2 ; V 2D is the velocity magnitude, m·s −1 ; σ U and σ W are standard deviations of the instantaneous air velocity in streamwise and vertical directions, respectively, m·s −1 .
It is noted that in this study V 2D , TI, and TKE were calculated based on a two-component velocity analysis.

Reynolds Number Independence Study
Dimensionless streamwise (U/U inlet ) and vertical (W/U inlet ) velocity profiles in the vertical direction (Z) at Line 1 in Figure 3 are compared in Figure 5 for different free stream wind speeds at the wind tunnel inlet (U inlet ), in order to identify the critical U inlet above which the flow can be considered Reynolds number independent [41]. The results showed that under U inlet of 4 m·s −1 and 6 m·s −1 conditions, the values of both U/U inlet and W/U inlet were different to the other cases. When U inlet increased from 6 m·s −1 to 8 m·s −1 , the average relative change of dimensionless air velocities was 0.391 and 0.006 for U/U inlet and W/U inlet , respectively. In contrast, when U inlet increased from 8 m·s −1 to 10 m·s −1 , the average relative change was only 0.003 and 0.001 for U/U inlet and W/U inlet , respectively. It indicated that when U inlet exceeded 8 m·s −1 , the dimensionless wind profile did not change considerably with further increase of the wind tunnel wind speed. Therefore, the Reynolds number independence was reached and a fully-developed turbulent flow was obtained at U inlet of 8 m·s −1 .

Reynolds Number Independence Study
Dimensionless streamwise (U/Uinlet) and vertical (W/Uinlet) velocity profiles in the vertical direction (Z) at Line 1 in Figure 3 are compared in Figure 5 for different free stream wind speeds at the wind tunnel inlet (Uinlet), in order to identify the critical Uinlet above which the flow can be considered Reynolds number independent [41]. The results showed that under Uinlet of 4 m·s −1 and 6 m·s −1 conditions, the values of both U/Uinlet and W/Uinlet were different to the other cases. When Uinlet increased from 6 m·s −1 to 8 m·s −1 , the average relative change of dimensionless air velocities was 0.391 and 0.006 for U/Uinlet and W/Uinlet, respectively. In contrast, when Uinlet increased from 8 m·s −1 to 10 m·s −1 , the average relative change was only 0.003 and 0.001 for U/Uinlet and W/Uinlet, respectively. It indicated that when Uinlet exceeded 8 m·s −1 , the dimensionless wind profile did not change considerably with further increase of the wind tunnel wind speed. Therefore, the Reynolds number independence was reached and a fully-developed turbulent flow was obtained at Uinlet of 8 m·s −1 .

Stability Study and Wind Profile Properties
To ensure a stable airflow within the region of interest, wind profiles at the upwind edge of the model (Line1) and the end of downwind airflow measurement region (Line 2) were measured without the scaled model. The vertical profiles of the streamwise air velocity component for Line 1 and Line 2 are compared in Figure 6. The results showed that the two profiles agreed well with each other with an average relative difference of 3.5%, indicating that a stable airflow was achieved within the region of interest.

Stability Study and Wind Profile Properties
To ensure a stable airflow within the region of interest, wind profiles at the upwind edge of the model (Line1) and the end of downwind airflow measurement region (Line 2) were measured without the scaled model. The vertical profiles of the streamwise air velocity component for Line 1 and Line 2 are compared in Figure 6. The results showed that the two profiles agreed well with each other with an average relative difference of 3.5%, indicating that a stable airflow was achieved within the region of interest. Figure 6 showed that the mean streamwise air velocity profile of the incident flow (i.e., at the upwind edge of the model) in the vertical direction followed a power law with the exponent of 0.14 and the coefficient of determination R 2 of 0.98:   Figure 6 showed that the mean streamwise air velocity profile of the incident flow (i.e., at the upwind edge of the model) in the vertical direction followed a power law with the exponent of 0.14 and the coefficient of determination R 2 of 0.98: where U and Uref are the mean streamwise air velocity at height Z and a reference height Zref, respectively, m·s −1 . To achieve the best fit of the model for the air velocity profile, Uref = 5.55 m·s −1 and Zref = 0.4 m were chosen in this study.
The streamwise air velocity at the building height (UB), calculated by Equation (4), was 4.77 m·s −1 . This resulted in a building Reynolds number (ReB) of 53,759, which was calculated by Equations (5) and (6).
where ν is the kinematic viscosity of the air, m 2 ·s −1 ; D is the characteristic length of the scaled pig building model, calculated as the hydraulic diameter of the cross-section of the model, m; WB and ZB are the width and the height of the scaled model, respectively, m. The average ambient temperature and relative humidity during the whole experimental period were 22.4 °C and 34.2%, respectively. Therefore, ν took the value of 1.53 × 10 −5 m 2 ·s −1 at temperature of 22.4 °C.
By fitting the air velocity profile of the incident flow to a logarithmic law, it gave the friction velocity (u*) of 0.22 m·s −1 and the full-scale roughness length (z0) of 6.4 × 10 −3 m, in which the von Karman constant took the value of 0.4. This gave the roughness Reynolds number (Rez0) of 92, which was calculated according to Equation (7).
It is stated in the VDI-guideline [40] that a moderately rough (corresponds to grassland or farmland) turbulent boundary layer should meet the following requirements: the profile power law exponent falls within the region of 0.12-0.18, and the roughness length within 0.005-0.1 m. In this study, both the power exponent of 0.14 and z0 of 6.4 × 10 −3 m satisfied the aforementioned criteria, hence the generated boundary layer can be considered to represent the airflow over farmland terrain. The streamwise air velocity at the building height (U B ), calculated by Equation (4), was 4.77 m·s −1 . This resulted in a building Reynolds number (Re B ) of 53,759, which was calculated by Equations (5) and (6).
where ν is the kinematic viscosity of the air, m 2 ·s −1 ; D is the characteristic length of the scaled pig building model, calculated as the hydraulic diameter of the cross-section of the model, m; W B and Z B are the width and the height of the scaled model, respectively, m. The average ambient temperature and relative humidity during the whole experimental period were 22.4 • C and 34.2%, respectively. Therefore, ν took the value of 1.53 × 10 −5 m 2 ·s −1 at temperature of 22.4 • C. By fitting the air velocity profile of the incident flow to a logarithmic law, it gave the friction velocity (u*) of 0.22 m·s −1 and the full-scale roughness length (z 0 ) of 6.4 × 10 −3 m, in which the von Karman constant took the value of 0.4. This gave the roughness Reynolds number (Re z0 ) of 92, which was calculated according to Equation (7).
It is stated in the VDI-guideline [40] that a moderately rough (corresponds to grassland or farmland) turbulent boundary layer should meet the following requirements: the profile power law exponent falls within the region of 0.12-0.18, and the roughness length within 0.005-0.1 m. In this study, both the power exponent of 0.14 and z 0 of 6.4 × 10 −3 m satisfied the aforementioned criteria, hence the generated boundary layer can be considered to represent the airflow over farmland terrain.
Additionally, both the building Reynolds number (Re B ) of 53,759 and the roughness Reynolds number (Re z0 ) of 92 were considerably higher than the reported critical Re B of 4000 [42] and critical Re z0 of 2.5 [43], respectively for a Reynolds number independent flow. It further indicated that the airflow generated in our wind tunnel at U inlet of 8 m·s −1 was fully developed turbulent.
Therefore, all subsequent measurements were performed at U inlet of 8 m·s −1 . Figure 7a shows contours of mean streamwise air velocity normalised by the air velocity at the building height (U B ) at Plane 1 for three roof slopes of the exercise yard of the scaled pig building model. Considerably lower air velocities were observed right behind the building with negative values in the centre (represented in dark blue in Figure 7a). The negative U velocities indicated reverse flows through Plane 1, which occurred slightly below the building height. The region of the reverse flow was affected by the roof slope in the spanwise wind direction (Y). It became narrower in the upper part but wider in the lower part for a steeper roof. Beyond the reverse flow region, the air flowed towards downstream with a reduced air velocity. Figure 7b compares the dimensionless mean vertical air velocity contours at Plane 1 for roof slopes of 5 • , 15 • , and 25 • . It was found that at both sides of the building, W velocities were negative, which indicated a downward flow direction. In contrast, within the region right behind the building upward airflows were observed. This indicated a complex three-dimensional air movement downwind the building caused by the air passing through and around the building. It is noted that the velocity contours were not strictly symmetric even though the scaled model was geometrically symmetric and the wind direction was perpendicular to the model sidewalls. The most possible reason for the asymmetrical airflow contours was the disturbance of the LDA probe to the flow field when it was positioned close to the building. The cross-sectional area of the LDA probe was 0.003 m 2 , accounting for 5% blockage to the airflow at the position right behind the building. This was not intended but was unavoidable, restricted by the focal length of the LDA. The potential influence of the air velocity sensor to the airflow was also reported by Sauer et al. [25], who observed an asymmetrical air velocity pattern at a vertical plane downstream from four aligned swine building models measured using a 3D hot-film anemometer. The results raised a potential interesting topic that is to quantify the disturbance of the LDA probe to the flow field, which could be done with CFD simulations.

Mean Air Velocities Downwind the Building
Appl. Sci. 2020, 10, x FOR PEER REVIEW 10 of 19 Additionally, both the building Reynolds number (ReB) of 53,759 and the roughness Reynolds number (Rez0) of 92 were considerably higher than the reported critical ReB of 4000 [42] and critical Rez0 of 2.5 [43], respectively for a Reynolds number independent flow. It further indicated that the airflow generated in our wind tunnel at Uinlet of 8 m·s −1 was fully developed turbulent.
Therefore, all subsequent measurements were performed at Uinlet of 8 m·s −1 . Figure 7a shows contours of mean streamwise air velocity normalised by the air velocity at the building height (UB) at Plane 1 for three roof slopes of the exercise yard of the scaled pig building model. Considerably lower air velocities were observed right behind the building with negative values in the centre (represented in dark blue in Figure 7a). The negative U velocities indicated reverse flows through Plane 1, which occurred slightly below the building height. The region of the reverse flow was affected by the roof slope in the spanwise wind direction (Y). It became narrower in the upper part but wider in the lower part for a steeper roof. Beyond the reverse flow region, the air flowed towards downstream with a reduced air velocity. Figure 7b compares the dimensionless mean vertical air velocity contours at Plane 1 for roof slopes of 5°, 15°, and 25°. It was found that at both sides of the building, W velocities were negative, which indicated a downward flow direction. In contrast, within the region right behind the building upward airflows were observed. This indicated a complex three-dimensional air movement downwind the building caused by the air passing through and around the building. It is noted that the velocity contours were not strictly symmetric even though the scaled model was geometrically symmetric and the wind direction was perpendicular to the model sidewalls. The most possible reason for the asymmetrical airflow contours was the disturbance of the LDA probe to the flow field when it was positioned close to the building. The cross-sectional area of the LDA probe was 0.003 m 2 , accounting for 5% blockage to the airflow at the position right behind the building. This was not intended but was unavoidable, restricted by the focal length of the LDA. The potential influence of the air velocity sensor to the airflow was also reported by Sauer et al. [25], who observed an asymmetrical air velocity pattern at a vertical plane downstream from four aligned swine building models measured using a 3D hot-film anemometer. The results raised a potential interesting topic that is to quantify the disturbance of the LDA probe to the flow field, which could be done with CFD simulations.   Figure 8 depicts contours of dimensionless streamwise and vertical air velocities at Plane 2 for the three roof slopes. Lower air velocities at the position behind the building than around the building were still observed. This indicated that the presence of the building affected the airflow field at a distance of 1.0 m (i.e., 7.7Z B , where Z B is the building height, corresponding to 50 m in full scale) downwind the building. In contrast to Plane 1, no reverse flows occurred at Plane 2 for all three cases. It was found that U velocities right behind the building at Plane 2 were obviously higher than those at Plane 1. However, Plane 1 presented higher air velocities at the height above 0.16 m than Plane 2 (Figures 7a and 8a), which was caused by the wind shear effect that accelerated the air when the air flowed over the top of the building. Compared with Plane 1, Plane 2 presented a more homogeneous vertical air velocity pattern with slightly upward airflows occurring near the centre (Figure 8b). Only small differences among different roof slope cases for both U and W velocities at Plane 2 were observed. It demonstrated that the slope of the leeward roof did not considerably affect the airflow at Plane 2.

Mean Air Velocities Downwind the Building
Appl. Sci. 2020, 10, x FOR PEER REVIEW 11 of 19 Figure 8 depicts contours of dimensionless streamwise and vertical air velocities at Plane 2 for the three roof slopes. Lower air velocities at the position behind the building than around the building were still observed. This indicated that the presence of the building affected the airflow field at a distance of 1.0 m (i.e., 7.7ZB, where ZB is the building height, corresponding to 50 m in full scale) downwind the building. In contrast to Plane 1, no reverse flows occurred at Plane 2 for all three cases. It was found that U velocities right behind the building at Plane 2 were obviously higher than those at Plane 1. However, Plane 1 presented higher air velocities at the height above 0.16 m than Plane 2 ( Figures. 7a and 8a), which was caused by the wind shear effect that accelerated the air when the air flowed over the top of the building. Compared with Plane 1, Plane 2 presented a more homogeneous vertical air velocity pattern with slightly upward airflows occurring near the centre (Figure 8b). Only small differences among different roof slope cases for both U and W velocities at Plane 2 were observed. It demonstrated that the slope of the leeward roof did not considerably affect the airflow at Plane 2. Airflow streamlines and contours of the air velocity magnitude (V2D) normalised by UB at Plane 3 for the three roof slopes are compared in Figure 9. The axes X and Z indicate the downwind distance to the building and the height from the wind tunnel floor, respectively. There was an elliptic-shaped low air velocity region with V2D/UB of smaller than 0.4 behind the building with a height of 0.13 m (i.e., 1ZB) and a downstream distance of 0.5 m (i.e., 3.8ZB and corresponding to 25 m in full-scale) for all three roof slopes. This region is expected to have a high concentration of air pollutants. This is because airborne pollutants are mainly transported by mean airflow, which results in the concentration of pollutants being inversely related to the air velocity downwind the building if no recirculation regions are presented [44]. If applying additional air pollutant treatment technologies, for example using pollutant traps or sprinklers to collect or wash high-concentrated air pollutants in this low air velocity region might effectively mitigate air contaminants and release their burden to the surrounding environment and residents. Additionally, within this low air velocity region, a wake zone with reversed airflow was observed at Plane 3, implying an anti-clockwise air recirculation in which air pollutants would accumulate. The size and shape of the wake was affected by the roof slope. For the roof slope of 5°, 15°, and 25°, the wake height (in Z direction) was 0.103 m, 0.099 m, and 0.084 m, respectively. Accordingly, the wake length (in X direction) was 0.254 m, 0.261 m, and 0.290 m, respectively. It showed that the larger the roof slope was, the lower and longer the wake became, and vice versa. As a result, the accumulated air pollutants would disperse farther for a steeper roof slope. The wake length observed by Tominaga et Airflow streamlines and contours of the air velocity magnitude (V 2D ) normalised by U B at Plane 3 for the three roof slopes are compared in Figure 9. The axes X and Z indicate the downwind distance to the building and the height from the wind tunnel floor, respectively. There was an elliptic-shaped low air velocity region with V 2D /U B of smaller than 0.4 behind the building with a height of 0.13 m (i.e., 1Z B ) and a downstream distance of 0.5 m (i.e., 3.8Z B and corresponding to 25 m in full-scale) for all three roof slopes. This region is expected to have a high concentration of air pollutants. This is because airborne pollutants are mainly transported by mean airflow, which results in the concentration of pollutants being inversely related to the air velocity downwind the building if no recirculation regions are presented [44]. If applying additional air pollutant treatment technologies, for example using pollutant traps or sprinklers to collect or wash high-concentrated air pollutants in this low air velocity region might effectively mitigate air contaminants and release their burden to the surrounding environment and residents. Additionally, within this low air velocity region, a wake zone with reversed airflow was observed at Plane 3, implying an anti-clockwise air recirculation in which air pollutants would accumulate. The size and shape of the wake was affected by the roof slope. For the roof slope of 5 • , 15 • , and 25 • , the wake height (in Z direction) was 0.103 m, 0.099 m, and 0.084 m, respectively. Accordingly, the wake length (in X direction) was 0.254 m, 0.261 m, and 0.290 m, respectively. It showed that the larger the roof slope was, the lower and longer the wake became, and vice versa. As a result, the accumulated air pollutants would disperse farther for a steeper roof slope. The wake length observed by Tominaga et al. [23] was 2.5Z B and 2.8Z B for a gable-roof building (without openings) with the roof slope of 16.7 • and 26.6 • , respectively, which is slightly greater than 2.0Z B and 2.2Z B for the roof slope of 15 • and 25 • respectively obtained in our study. This is because the sealed building structure could result in a much lower pressure field behind the building and thus a larger air recirculation region. In contrast to the results in this paper where the wakes were found from the floor level until above the roof, in a wind tunnel study of Ikeguchi and Okushima [18] using naturally ventilated dairy building models without outdoor exercise yards, the wakes were only observed above the leeward roof. One possible reason could be the large sidewall openings used in their study, which resulted in a cross ventilation through the building. In contrast, in our study, a small slot inlet opening with a high sidewall led to an up-jet airflow pattern inside the housing area. Additionally, the exercise yard with an upwind sidewall further complicated the airflow and resulted in the reversed flow behind the building. Therefore, our results are expected to be applicable to naturally ventilated livestock buildings with small inlet openings and downwind outdoor yards in a perpendicular wind direction.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 12 of 19 al. [23] was 2.5ZB and 2.8ZB for a gable-roof building (without openings) with the roof slope of 16.7° and 26.6°, respectively, which is slightly greater than 2.0ZB and 2.2ZB for the roof slope of 15° and 25° respectively obtained in our study. This is because the sealed building structure could result in a much lower pressure field behind the building and thus a larger air recirculation region. In contrast to the results in this paper where the wakes were found from the floor level until above the roof, in a wind tunnel study of Ikeguchi and Okushima [18] using naturally ventilated dairy building models without outdoor exercise yards, the wakes were only observed above the leeward roof. One possible reason could be the large sidewall openings used in their study, which resulted in a cross ventilation through the building. In contrast, in our study, a small slot inlet opening with a high sidewall led to an up-jet airflow pattern inside the housing area. Additionally, the exercise yard with an upwind sidewall further complicated the airflow and resulted in the reversed flow behind the building. Therefore, our results are expected to be applicable to naturally ventilated livestock buildings with small inlet openings and downwind outdoor yards in a perpendicular wind direction.  Figure 10 shows the variations of the dimensionless air velocities with distance to the building at heights (Z) of 0.02 m, 0.06 m, and 0.12 m, which corresponded to 1 m, 3 m, and 6 m in full-scale. At Z of 0.02 m, the roof slope of 25° resulted in lower U/UB and W/UB velocities than roof slopes of 5° and 15° until a distance of around 0.7 m from the building. At Z of 0.06 m and 0.12 m, there were not many differences for U/UB velocities among the three roof slopes. However, for W/UB velocities, the 5° roof slope presented the highest values among three cases until a distance of 0.6 m at Z of 0.06 m and a distance of 0.8 m at Z of 0.12 m, respectively, indicating a more upwards airflow direction. At a distance beyond 0.8 m (i.e., 6.2ZB) from the building, the roof slope had almost no effect on the air velocities. This is in line with the previous observations at Plane 2 shown in Figure 8 where no considerable differences could be detected among three roof slope cases.  At a distance beyond 0.8 m (i.e., 6.2Z B ) from the building, the roof slope had almost no effect on the air velocities. This is in line with the previous observations at Plane 2 shown in Figure 8 where no considerable differences could be detected among three roof slope cases.

Air Turbulence Downwind the Building
Spatial distributions of the turbulent kinetic energy (TKE) normalised by the square of the streamwise air velocity at the building height (UB) at Plane 1 and Plane 2 for roof slopes of 5°, 15°, and 25° are compared in Figure 11. It was found that TKE had different patterns from air velocities, where the region of low air velocities generally presented higher TKE values. At Plane 1, the highest TKE occurred at the centre of the building around the building height. This might be caused by the air collision of the reverse flow occurred in the wake and the bulk air flowed over the building, which increased the air turbulence in that region. At Plane 1, the high TKE region was the smallest when the roof slope was 5°. At Plane 2, higher TKE was observed behind and above the building than in other areas. It indicated that the influence of the building on the downwind air turbulence reached Plane 2. Only minor TKE pattern differences among different roof slopes were observed at Plane 2, which implied that the roof slope had almost no effects on TKE at Plane 2.

Air Turbulence Downwind the Building
Spatial distributions of the turbulent kinetic energy (TKE) normalised by the square of the streamwise air velocity at the building height (U B ) at Plane 1 and Plane 2 for roof slopes of 5 • , 15 • , and 25 • are compared in Figure 11. It was found that TKE had different patterns from air velocities, where the region of low air velocities generally presented higher TKE values. At Plane 1, the highest TKE occurred at the centre of the building around the building height. This might be caused by the air collision of the reverse flow occurred in the wake and the bulk air flowed over the building, which increased the air turbulence in that region. At Plane 1, the high TKE region was the smallest when the roof slope was 5 • . At Plane 2, higher TKE was observed behind and above the building than in other areas. It indicated that the influence of the building on the downwind air turbulence reached Plane 2. Only minor TKE pattern differences among different roof slopes were observed at Plane 2, which implied that the roof slope had almost no effects on TKE at Plane 2. Black dash lines indicate the profile of the scaled model. Figure 12 illustrates the normalised TKE distributions at Plane 3 for three roof slopes. The highest TKE was observed in the region in the vicinity of the building height. It decayed gradually as the distance to the building increased. The high TKE region that is represented in red and yellow in Figure 12 extended up to approximately 0.5 m (i.e., 3.8ZB) downwind from the building. It was located right above the region in which the air velocity was low (as depicted in Figure 9). This could be attributed to the flow reattachment on the leeward roof that led to the production of TKE. Ntinas et al.
[45] also found a high turbulence production in the air recirculation region behind the scaled model. It is known that the transport of airborne pollutants is affected by both the mean air velocity and the air turbulence [46]. Therefore, in that low air velocity region, the aerial pollutants would be likely transported via the energetic turbulent eddies/turbulence diffusion.   Figure 12 illustrates the normalised TKE distributions at Plane 3 for three roof slopes. The highest TKE was observed in the region in the vicinity of the building height. It decayed gradually as the distance to the building increased. The high TKE region that is represented in red and yellow in Figure 12 extended up to approximately 0.5 m (i.e., 3.8Z B ) downwind from the building. It was located right above the region in which the air velocity was low (as depicted in Figure 9). This could be attributed to the flow reattachment on the leeward roof that led to the production of TKE. Ntinas et al.
[45] also found a high turbulence production in the air recirculation region behind the scaled model. It is known that the transport of airborne pollutants is affected by both the mean air velocity and the air turbulence [46]. Therefore, in that low air velocity region, the aerial pollutants would be likely transported via the energetic turbulent eddies/turbulence diffusion. Black dash lines indicate the profile of the scaled model. Figure 12 illustrates the normalised TKE distributions at Plane 3 for three roof slopes. The highest TKE was observed in the region in the vicinity of the building height. It decayed gradually as the distance to the building increased. The high TKE region that is represented in red and yellow in Figure 12 extended up to approximately 0.5 m (i.e., 3.8ZB) downwind from the building. It was located right above the region in which the air velocity was low (as depicted in Figure 9). This could be attributed to the flow reattachment on the leeward roof that led to the production of TKE. Ntinas et al.
[45] also found a high turbulence production in the air recirculation region behind the scaled model. It is known that the transport of airborne pollutants is affected by both the mean air velocity and the air turbulence [46]. Therefore, in that low air velocity region, the aerial pollutants would be likely transported via the energetic turbulent eddies/turbulence diffusion.  The influences of the roof slope on the air turbulence, denoted by the turbulence intensity (TI) and dimensionless turbulent kinetic energy (TKE/U 2 B ), with respect to the downwind distance from the building (X) at heights (Z) of 0.02 m, 0.06 m, and 0.12 m are described in Figure 13. A considerably higher TI was found at X of 0.2 m-0.4 m for Z of 0.02 m and 0.06 m, where the air vortex was formed in the wake as depicted in Figure 9. The highest TI was observed in the most downward roof structure (i.e., roof slope of 25 • ), which was attributed to an enhanced disturbance of the roof to the airflow. At Z of 0.12 m, the TI had lower values compared with Z of 0.02 m and 0.06 m, and decreased slowly with distance from the building. The influence of the roof slope on the TI was negligible when the distance from the building exceeded 0.5 m (i.e., 3.8Z B ) (Figure 13a). The highest TKE/U 2 B occurred at Z of 0.12 m and decayed gradually away from the building. Even though the roof slope had a big impact on the TI, no notable differences were observed among different roof slope cases for TKE/U 2 B , indicating that the variations in the leeward roof slope might have little influence on the turbulent kinetic energy of the airflow downwind the building. The same TKE/U 2 B values for different cases implied that the velocity fluctuations, as explained by Equation (3), were not affected by the roof slope. However, the velocity magnitude (V 2D ), which was calculated from the U and W velocities (Equation (1)) were affected by the roof slope, as seen in Figure 10. According to Equation (2), TI was calculated as the ratio of the square root of TKE to V 2D . Therefore, different TI values were observed among different roof slope cases although the same TKE/U 2 B occurred. The influences of the roof slope on the air turbulence, denoted by the turbulence intensity (TI) and dimensionless turbulent kinetic energy (TKE/U B 2 ), with respect to the downwind distance from the building (X) at heights (Z) of 0.02 m, 0.06 m, and 0.12 m are described in Figure 13. A considerably higher TI was found at X of 0.2 m-0.4 m for Z of 0.02 m and 0.06 m, where the air vortex was formed in the wake as depicted in Figure 9. The highest TI was observed in the most downward roof structure (i.e., roof slope of 25°), which was attributed to an enhanced disturbance of the roof to the airflow. At Z of 0.12 m, the TI had lower values compared with Z of 0.02 m and 0.06 m, and decreased slowly with distance from the building. The influence of the roof slope on the TI was negligible when the distance from the building exceeded 0.5 m (i.e., 3.8ZB) (Figure 13a). The highest TKE/U B 2 occurred at Z of 0.12 m and decayed gradually away from the building. Even though the roof slope had a big impact on the TI, no notable differences were observed among different roof slope cases for TKE/U B 2 , indicating that the variations in the leeward roof slope might have little influence on the turbulent kinetic energy of the airflow downwind the building. The same TKE/U B 2 values for different cases implied that the velocity fluctuations, as explained by Equation (3), were not affected by the roof slope. However, the velocity magnitude (V2D), which was calculated from the U and W velocities (Equation (1)) were affected by the roof slope, as seen in Figure 10. According to Equation (2), TI was calculated as the ratio of the square root of TKE to V2D. Therefore, different TI values were observed among different roof slope cases although the same TKE/U B 2 occurred.