SIMULATION OF VENTILATION SYSTEMS IN A PROTECTED ENVIRONMENT USING COMPUTATIONAL FLUID DYNAMICS

Computational simulations of mass and energy flow help in implementing alternative cooling systems in protected environments. The aim of this study was to model and simulate the interaction between external and internal environments of a protected environment by means of computational fluid dynamics (CFD) techniques and validate micrometeorological variables for subsequent comparison between natural and indirect ventilation by ground heat exchangers. At the first phase, the micrometeorological variables global solar radiation (Qg), air temperature (Tair), and air relative humidity (RH) were monitored. The second phase consisted of the numerical modeling of finite volumes, with validation through recorded data, as well as simulation and comparison of two ventilation systems. The functional relationship between simulated and recorded meteorological elements presented a good linear association, with coefficients of determination of 0.97, 0.93, and 0.94 for Qg, Tair, and RH, respectively. Simulation of indirect ventilation system by ground heat exchangers presented a reduction of 4 °C in Tair and 15% in RH compared to that recorded inside the environment. The natural ventilation system allowed a reduction of 1 °C in Tair when compared to the protected environment.


INTRODUCTION
Protected cultivation is one of the technologies that has contributed to agriculture modernization since micrometeorological variables such as air temperature, air relative humidity, wind, and solar radiation are modified by plastic films, shading meshes, and climate conditioning systems.These systems promote a positive effect on plant production, development, and growth, in addition to reducing the effect of productive seasonality during the year (ROMANINI et al., 2010).
Monitoring of meteorological variables and physical characterization of the protected environment in general demands a high investment for the acquisition of instruments that allow recording and later analyzing the data that assist in decision-making.In addition, the adoption of abiotic factor control systems requires the actions of specialists aiming at minimizing possible interference in protected cultivation.Faced with an increasingly competitive market, this scenario is hardly available to producers and requires the insertion of methodologies and computational tools that can guarantee a greater interaction between the different actors in the agricultural sector, better use of available resources, and productive efficiency (ABDEL-GHANY & HELAL, 2011).
Thus, the use of mathematical models of computational fluid dynamics (CFD) can be used in a simulation environment, in a practical and applied way, to analyze the effect and management of the most diverse materials on meteorological variables inside the protected environment (ALI et al., 2014).
CFD technique allows simulating fluid interaction by means of a computer graphics interface, with surfaces defined by boundary condition through fundamental equations of conservation of matter, energy, and momentum (TORRE-GEA et al., 2011;ROMERO-GÓMEZ et al., 2010).BOURNET & BOULARD (2010) pointed out that, associated with mass and energy flow, this technique provides subsidies to estimate micrometeorological conditions in response to an air renewal system, or ambient cooling, in order to assess the configurations of ventilation systems under different climatic conditions.Ventilation in a protected environment consists of replacing the mass of heated air from the inside with the outside air.This allows evacuating much of the heat overload, reducing air temperature as water vapor and gases are dissipated.Ventilation systems adopted can present a natural or mechanical configuration.Natural ventilation requires central or lateral openings, or even the combination of both, being based on convection and advection energy transfer (BOURNET & BOULARD, 2010).
In addition to the direct forced ventilation methods, with the use of ventilation or exhaust fans, stands out the indirect air circulation system by ground heat exchangers, consisting of a bundle of tubes buried in the ground based on the model of MONGKON et al. (2013).The forced convection of air inside the tube performs the necessary thermal exchanges to keep the protected environments under desired temperature ranges (GARCIA, 2001).
In this context, the aim of this study was to model and simulate the interaction between external and internal environments of a protected environment by means of computational fluid dynamics (CFD) techniques and validate micrometeorological variables for subsequent comparison between natural and indirect ventilation by ground heat exchangers.

MATERIAL AND METHODS
CFD modeling data were recorded between June 21 and August 2, 2013, allowing the characterization of the day and time with greater availability of solar radiation (simulated day: 6/25/2013, 11:30 h) inside a protected environment and in the external environment.The measurements were made by automatic meteorological stations located at the Federal Rural University of Pernambuco, Recife -PE, Brazil (08°04'03" S, 34°55'00" W, and 4-m altitude).Local climate is characterized as megathermal (As`) according to Köppen classification, with winter precipitation and average annual temperature of 25.2 °C (PEREIRA et al., 2002).
The protected environment contained micropropagated sugarcane seedlings growing in trays set on benches.The internal floor was composed of sandy soil with the following dimensions and constructive characteristics: 7.0 m width, 21.0 m length, 3.0 m ceiling height, 4.5 m total height, arc-shaped ceiling covered with low-density polyethylene plastic (150 µm), being closed laterally with an anti-aphid screen with 3.2 µm porosity.A thermo-reflective screen (50% shading) was installed close and internally to the plastic cover.
In order to determine the boundary conditions and generation of simulations, the following internal and external meteorological variables were recorded in the protected environment: solar radiation (Qg, W m −2 ; MJ m −2 d −1 ), air temperature (Tair, °C), air relative humidity (RH, %), and wind velocity (Wv, m s −1 ).Protected and external environments were monitored by two automatic meteorological data logging stations (Datalogger Campbell ® , model CR1000).The first station was located at the geometric center of the protected environment, with sensors arranged at a height of 1.5 m from the ground, and the second station was installed 200 m from the protected environment (agrometeorological station).Both stations were programmed to register meteorological data every 5 seconds, generating averages every 15 minutes, 60 minutes, and 24 hours.
CFD was based on equations that govern fluid dynamics (continuity, momentum, and energy), obtained directly from a volume, or fixed element in a conservative system.
In the continuity equation, mass flows that enter in a fluid element must match exactly with those leaving (Equation 1).
(1) Momentum equation (conservative mode) admits that the sum of external forces acting on the fluid particle is equal to its linear momentum rate of change (Equation 2).
(2) Energy conservation (second law of thermodynamics) assumes that the energy variation rate of a fluid particle is equal to the addition of heat and work on the particle (Equation 3).
(3) CFD techniques apply these conservation laws along the discretized flow domain in order to determine the systematic changes of mass, momentum, and energy for the fluid that crosses the interfaces of each discrete region (TORRE-GEA et al., 2011).
To describe the solar radiation intensity in transparent materials, such as the plastic film of the cover used in the protected environment, the discrete ordinate method was used, with emphasis on the radiative transfer equation solution (Equation 4).
(4) where, r is the position vector, l is the direction, Kλ is the monochromatic absorption coefficient, Iλ is the monochromatic intensity, and Ibλ is the blackbody intensity given by the Planck function.
CFD modeling was divided into three stages: pre-processing, processing, and post-processing (NORTON et al., 2007).
Boundary conditions were applied according to critical situations analyzed during data recording period.Thus, day and time of greater radiation availability were selected for use in CFD models in order to generate the initial conditions of the protected environment and for comparison of ventilation systems and their configurations.
The geometric domain of calculation was discretized in a non-structured mesh of the surrounding protected environment to characterize the air boundary layer assuming 25 m wide and 14 m high, associated with the protected environment geometry (Table 1).Data validation was established by comparison between data registered in the protected environment with those simulated.A sequence of simulations was considered.In this case, the average differences between recorded and simulated values were determined, as well as the functional relationship between them, established by regression analysis for simulated day and time (6/25/2013, 11:30 h).
Mathematical modeling, through CFD, was used to distribute the airflows to the protected environment studied, as well as simulate, predict, and compare two ventilation systems: a natural ventilation system with an opening at the top of the coverage structure, and an indirect cooling by ground heat exchangers aiming at estimating its applicability in planning and managing micrometeorological variables.For all simulations, the maximum ventilation rates calculated by the software ANSYS FLUENT v. 15 were considered.
Natural ventilation system consisted of a 50-cm moving window arranged on the opposite side of the air boundary layer entrance, having as its boundary condition the air outflow.
Simulation with indirect ventilation system by ground heat exchangers was proposed based on the technique developed by MONGKON et al. (2013) and adapted to simulation conditions.The system consisted of tubes buried in the ground, being the air intake positioned on one side of the agricultural greenhouse, going through the buried tubes to a depth of 1 m.After cooled, the air was discharged at equidistant longitudinal points in the protected environment.
The 2D simulation of the transversal profile of the protected environment corresponded to one of these points of cooled air diffusion.The aluminum tube was configured with a diameter of 0.08 m, double flow discharge with input velocity into the tube of 1 m s −1 , and RH of 73%.The initial 10 cm of the tube were configured in the software as moving wall at an air velocity of 0.1 m s −1 and at the tube outlet was configured as air outflow.
Computational mesh geometry was generated with an orthogonality of 0.74, considered as adequate according to FLUENT (1998), who classifies the orthogonality with values from 0 to 1; the closer to zero this value, the lower the mesh quality is (Figure 2).The mesh was classified as semi-structured, composed of 37,533 knots and 36,176 elements; each element represented 10 cm 2 of the 350 m 2 of the total mesh.From the 2000 thousand programmed interactions, a large part of them converged between 250 and 300 interactions, with an average processing time of 5 minutes for each process.This corroborates with BOUHOUN-ALI et al. (2014), who tested different amounts of elements of computational mesh for a protected environment of 5.9 m in height and 9.6 m in width and observed that 36,900 elements ensured satisfactory results.2010), the advection effect acts close to lateral and ground meshes due to the effect of pressure variation to windward and leeward.Conduction and irradiation processes act from the plastic cover and materials present inside the protected environment, considering the associated effect of temperature variation on the vertical profile and containment of heated air mass (Figure 3a).
At the center of the protected environment occurs the stabilization of Tair values when compared to those close to the coverage, which is conditioned by the radiation transmission incident on the plastic film.On the other hand, in the lower part, boundaries present variations close to the lateral meshes, where the processes of convection and advection energy transfer are more active (Figure 3a).On average, the difference between Tair close to the coverage and that close to the ground was 1.2 °C.According to BAXEVANOU et al. (2010), the air temperature in the upper part of the protected environment is a function of radiation incidence and its angle; at the lower part, the air temperature is totally independent of the angle of incidence of Qg and dependent on convection and advection processes.
Figure 3b shows the vertical variation of global solar radiation (Qg, W m −2 ), presenting a variation in intensity and distribution along the transversal profile of the protected environment.This response was also observed by BAXEVANOU et al. (2010), who verified variations above 200 W m −2 attributed to air conditions inside the protected environment such as water content, which changes the absorption and refraction properties of the air and, consequently, interferes with solar energy transfer into the protected environment.
When analyzing the Qg distribution profile in the protected environment, a diffuse energy uniformity at the crop level is observed in the horizontal profile.According to FIDAROS et al. (2010), the knowledge of Qg distribution profile inside the protected environment assists in the management of plant positioning, increasing the use of radiant energy and photosynthetic efficiency of the crop.Figure 4a shows the variation of air relative humidity (RH, %) in the transversal profile of the protected environment.RH presents higher values at air inlet and outlet regions, close to the lateral meshes.In the lower parts from the center to the top of the structure, RH is reduced.This spatial variability is due to the ascending flow of air mass conditioned by the free convection concomitant to the increasing Tair in the vertical profile.BOURNET & BOURLARD (2010) observed from CFD simulations that for a Wv lower than 1.0 m s −1 , the thermosyphon effect (convection) is more active than the dynamic wind forces (advection), in which the hot air, less dense and with less water vapor, occupies the upper layers of the protected environment.The functional relationship between Qg measured inside the protected environment with values obtained for each numerical simulation by using computational fluid dynamics at the time of simulated day (9:30, 10:00, 11:30, 12:30, 13:00, 13:15, and 14:00) is shown in Figure 5a.The coefficient of determination was high (R 2 = 0.97), which indicates a good linear association.The angular coefficient was 0.97, i.e. the simulations presented an average error of 3% of the values obtained by the meteorological station inside the protected environment.
Tair values recorded inside the protected environment and those obtained by simulation presented a high coefficient of determination (R 2 = 0.93), indicating a good linear association (Figure 5b).The angular coefficient was 1.01; therefore, the simulations presented an error lower than 1% of the values obtained by the meteorological station inside the protected environment.
Variations between Tair measured and simulated was between 0.1 and 0.63 °C for the respective times studied, with an average of 0.21 °C.This fact can be explained by the absence of more detailed parameters referring to the soil-plant-atmosphere relationship.According to BOURNET & BOULARD (2010), the spatial heterogeneity of micrometeorological variables inside a protected environment considerably interferes with crop activities such as plant gas exchange, particularly transpiration and photosynthesis.These processes, in turn, can also alter the exchange of sensible and latent heat.MONTERO et al. (2013) used 2D modeling to assess the thermal inversion in protected environments under heating conditions with internal and external curtain and obtained thermal inversion values from 2.5 to 3.5 K for completely clear and cloudy nights.
Functional relationships between RH monitored and simulated (Figure 5c) presented a good association (R 2 = 0.94), with simulated results representing an average of 97% of those obtained in the monitoring and close to the angular coefficient (y = 0.9645x).VOLTAN et al. (2014) reported that RH is directly influenced by solar radiation and presents a spatial dependence related to the environment height, with an amplitude of 7.8% observed for the relative humidity.
In addition, KIM et al. (2008) observed that 2D models might present difficulties in simulating RH, generally underestimating them.This is due to the instability of water vapor movement, which presents a variable spatial movement.One solution to this problem would be the use of 3D geometry.However, ALI et al. ( 2014) pointed out the need for equipment with high processing capacity, but even so, the simulations would require much more time to obtain the graphical response interface, which could last for days.Figure 6 shows the simulated Wv vectors and their distributions in the protected environment profile.The natural ventilation system with an opening in the ceiling (Figures 6a1 and 6a2) presents an upward flow towards the opening, being this characteristic configured as an air outflow and the right side of the atmospheric boundary layer configured with the pressure outlet.The increased surface area of openings favors protected environments with a greater ceiling height, in which the thermal fluctuation has less abrupt variations when compared to protected environments with lower air volume (TORRE-GEA et al., 2014).
The ascending airflow in the natural ventilation system is considered by the initial condition in which the internal Tair in the simulated time is greater than the external Tair.Thus, the difference in temperature causes a difference in pressure between both environments, with the airflow Eng.Agríc., Jaboticabal, v.37, n.3, p.414-425, maio/jun. 2017 421 increasing towards the opening by the action of the thermosyphon effect.The opening in the coverage favored the air convection process, as well as the renewal of the humid air mass, and the internal RH can match that recorded externally (VAZQUEZ, 2013).
The maximum ventilation rate calculated by the software was 3.0 m 3 s −1 .KATSOULAS et al. ( 2006), considering anti-aphid screens and different opening formats to promote natural ventilation in a protected environment, found values ranging up to 5 m 3 s −1 for systems with lateral opening.OLIVEIRA et al. ( 2014) reported a maximum air velocity of 2.0 m 3 s −1 in a protected environment, evidencing that the greatest problems are related to low velocities and hence a less air renewal in the environment.ESPINAL-MONTES et al. (2015) reported that higher velocities in a protected environment are observed close to the openings and in the region opposite to them.
Vector distribution in the simulation of system by ground heat exchanger (GHE) is shown in Figure 6 (b1 and b2).The flow is ascending from the air intake (pipe on the right side) to the middle of the protected environment and, subsequently, is descending towards the air outflow (pipe on the left side).The air intake velocity was 0.75 m s −1 .ROMERO-GÓMEZ et al. (2010) observed in CFD simulations higher ventilation rates with wind velocities above 1.4 m s −1 .FIGURE 6. Vectors indicating air velocity (1) and flow distribution (2) in the natural ventilation system (a) and ground heat exchanger (b).
In the system with an open on the protected environment coverage, a reduction of 1 °C was observed in Tair value registered internally.Thus, the simulated Tair inside the environment was equal to the Tair recorded at the external meteorological station (Figure 7a1).Figure 7a2 shows that in GHE system, the difference between internal and simulated Tair was 4 °C, being attributed mainly to a low heat conduction of soil.MATHUR et al. (2015) reported that this system provides air-cooling inside the protected environment, where soil temperature, because it is lower than the ambient temperature during the day due to thermal conductivity, cools the tubing and, consequently, the air that passes inside the buried pipes.However, the authors warned that the continuous use of this system could cause thermal saturation in the soil as the day passes since the soil accumulates heat and the dissipation is slow.
GHE system was able to dehumidify 15% in relation to the recorded RH (Figure 7b2), and in the air intake (right side), the simulated RH was 0.23.The new air mass was mixed with that present in the environment.Thus, the simulated RH value was gradually changed due to the presence of water vapor inside the environment; at the air outflow, the simulated RH presented a maximum value (0.77), but still lower than the recorded RH (0.81).

CONCLUSIONS
The use of fluid dynamics with a computational graphical interface (2D) allowed simulating accurately the variability of global solar radiation, air temperature, and air relative humidity in the transversal profile of the protected environment, as well as its conditioning effects on the airflow.
The simulation of natural ventilation system presented a reduction in air temperature and the indirect ventilation system by ground heat exchangers was more efficient in attenuating air temperature and air relative humidity when compared to the values recorded inside the protected environment.

FIGURE 2 .
FIGURE 2. Protected environment representation and its surroundings (air boundary layer) by means of a semi-structured computational mesh.

FIGURE 7 .
FIGURE 7. Distribution of Tair and RH in the transversal profile of the protected environment in the natural ventilation system (a) and ground heat exchanger (b).

TABLE 1 .
Details of boundary conditions.

TABLE 2 .
Values of the initial variables used in boundary conditions for each simulation.