Study using a CFD approach of the efficiency of a roof ventilation closure system in a multi-tunnel greenhouse for nighttime microclimate optimization

Submitted on December 1 1, 2019 and accepted on August 13, 2020. 1 The present study was funded by the Sistema General de Regalías and the Gobernación de Cundinamarca through the project “Fortalecimiento de la competitividad del sector floricultor colombiano mediante el uso de ciencia, tecnología e innovación aplicadas en el departamento de Cundinamarca” code: BPIN2013000100217. 2 Universidad Jor ge Tadeo Lozano, Facultad de Ciencias Naturales e Ingeniería, Departamento de Ciencias Biológicas y Ambientales, Bogotá, Colombia. edwina.villagranm@utadeo.edu.co 3 Universidad Jor ge Tadeo Lozano, Facultad de Ciencias Naturales e Ingeniería, Departamento de Ciencias Básicas y Modelado, Bogotá, Colombia. carlos.bojaca@utadeo.edu.co Corresponding autor: edwina.villagranm@utadeo.edu.co Study using a CFD approach of the ef ficiency of a roof ventilation closure system in a multi-tunnel greenhouse for nighttime microclimate optimization 1


INTRODUCTION
In Colombia, greenhouse horticultural production has an approximate area of 10000 ha, dedicated almost exclusively to cut flowers and tomato crops. These production systems are characterized by the use of low cost, low-tech, naturally ventilated greenhouses with limited possibilities to control of the microclimate generated inside the cultivation area.
The climate control in naturally ventilated greenhouses is mainly carried out by passive means, a method that does not require any extra equipment or external energy to regulate the internal climate. The ventilation is carried out by opening or closing of the side and roof vents of the greenhouse. Generally, the roof ventilations are fixed, that is, they remain open and in constant exchange of air with the external environment. During the day period, the microclimatic conditions of the greenhouse are determined by the efficiency of the natural ventilation of the greenhouse. This efficiency depends on the characteristics of the structure and its surroundings, which have already been widely studied worldwide (Bournet et al., 2006;Rojano et al., 2014;Espinoza et al., 2017;McCartney & Lefsrud, 2018;Villagrán et al., 2019;Villagrán & Bojacá, 2019b;Villagrán & Bojacá, 2019c;Akrami, et al., 2020a). In Colombia, the ventilation efficiency of the greenhouses is low because does not comply with the minimum recommended value of 40 air renewals h -1 ; that is, the minimum value from which a passive greenhouse has the capacity to regulate the thermal and humidity excesses of the air (Montero, 2006;Villagrán & Bojacá, 2019d). During the night hours, the climate control is limited to the total closure of the side vents, aiming at keeping most of the energy gained during the day and maintaining as high as possible the temperature of the cultivation area. In this way, the objective is to prevent the thermal inversion of the greenhouse, a phenomenon in which the temperature of the greenhouse is lower than that of the surrounding outdoor air (Villagrán & Bojacá, 2019a). One of the causes of this phenomenon is the limited capacity of the greenhouses to store the heat captured during the day, since this is lost in a very fast process during the early hours of the night, especially through the roof ventilation (Teitel et al., 2008). These low temperatures limit the potential production of the crop due to the reduction of their development rates  and in more critical situations they can cause irreversible crop losses when extreme events such as frost occur.
The microclimatic study of a greenhouse can be developed through in-situ experimentation or scale models, approximations through analytical models of mass and energy balances or through numerical simulation through methodologies such as computational fluid dynamics (CFD) (Akrami et al., 2020b;Chen, 2009). The latter methodology has been widely used in various branches of engineering over the past four decades to analyze problems involving fluid flow and its interactions with heat and mass transfer. (Mohammadrezaei et al., 2017). In the agricultural sector, this approach has been applied to the aerodynamic study of different types of greenhouses (Bartzanas et al., 2004;Benni et al., 2016;Tong & Christopher, 2018).
Several of these works have included the use of the discrete-order radiation model (DO) that allows modeling the energy flow by radiation in a greenhouse . The CFD technique coupled with the DO model has been used to analyze the implementation of passive alternatives for night climate control in naturally ventilated greenhouses, among others. Montero et al. (2005) studied the nocturnal thermal distribution and energy flows of a greenhouse in northern Patagonia under three atmospheric conditions and two types of coverings. The use of the simulation model made it possible to compare the different passive energy saving strategies, showing the positive effect of increased air temperature through the use of double plastic cover on the roof. Montero et al. (2013) carried out 2D simulations for a passive greenhouse, considering alternative scenarios to evaluate the use of external and internal curtains. The authors reported that on clear and dry nights the phenomenon of thermal inversion was presented under the three simulated scenarios. Espinal-Montes et al. (2015) developed a work based on a numerical model CFD 3D with the objective to evaluate the nocturnal thermal behavior of a passive Mexican greenhouse, the main results of this study reported the generation of the phenomenon of thermal inversion with average gradients of -0.8 and -3.1 °C for the configuration of open vents and for the configuration totally closed respectively. Other works related to the night-time condition (Bournet et al., 2006;Majdoubi et al., 2016;Majdoubi et al., 2017) analyzed the nocturnal heat loss according to the radiative and convective energy losses of the roof and the effect of the use of a thermal cover. In the Colombian context, a study developed by Villagrán & Bojacá, (2019e), by means of 24 twodimensional numerical simulations analyzed the effect of a thermal screen and a porous screen on the thermal behavior of a gothic tunnel greenhouse. The use of these screens generated higher thermal values for cloudy and clear sky conditions compared to the standard configuration where these types of screens are not used.
In order to explore night-time climate control alternatives that improve the thermal behavior of naturally ventilated greenhouses, the objective of this paper is to evaluate a system of inflatable air ducts that close the fixed roof openings during the night in a multi-tunnel type greenhouse located in the savannah of Bogotá. This evaluation is carried out using a two-dimensional CFD simulation model that allows the analysis of thermal and air flow behaviors, which were later validated by the construction of a full-scale prototype of the greenhouse.

Numerical simulation
The CFD numerical simulation methodology allows the calculation of the air flow pattern and the thermal distribution field inside a greenhouse, behavior that is governed mainly by the Navier-Stokes equations. In this work a commercial software of finite volumes was used to solve the two-dimensional equations in stationary state and that can be represented by means of the following equation of diffusion transport by convection: (1) where is the velocity vector (m s -1 ), p is the fluid density (kg m -3 ), Γ represents the diffusion coefficient (m 2 s -1 ), is the nabla operator, φ represents the concentration of the transported amount in a dimensional form (the momentum conservation equation, the scalars of the mass and the energy conservation equations) and S indicates the source term (Piscia et al., 2015). To evaluate the turbulent nature of the airflow the standard model k-ε was adopted, this model has been widely used and validated in similar studies, the results obtained have proven to be reliable and on the other hand this model demands lower computational requirements The turbulent nature of the flow is simulated through the application of one of the Reynolds averaged turbulence models (RANS) available in the simulation software. The model selected was the k-ε standard which has been used and validated in many studies like this research. This model allows reliable simulation results to be obtained with a lower computational work cost than would be required if another turbulence model were used (Katsoulas et al., 2006;Bournet et al., 2017). Rev. Ceres, Viçosa, v. 67, n.5, p. 345-356, sep/oct, 2020 The air buoyancy produced by density changes generated by temperature variations was coupled to the simulation using the energy equation and the Boussinesq approximation. The radiation model selected for the nighttime climate simulation was the discrete order (DO) model. This model allows the simulation of radiative processes in environments that include semi-transparent walls or roofs, such as greenhouses. This DO model is also one of the most widely used to model the radiation flow that occurs under night-time climate conditions inside a greenhouse without heating. This flow is characterized by a phenomenon of thermic radiation from the floor surface that dissipates towards the walls, roof and external environment of the greenhouse (Montero et al., 2013;Piscia et al., 2015). The DO radiation model is represented in the numerical model by the following equation: ( 2) where is the intensity of the radiation at a wavelength; are the vectors that indicate the position and direction, respectively; is the direction of dispersion vector; , are the dispersion coefficients and the spectral absorption, respectively; n is the refractive index; represents the divergence operator; σ is the Stefan-Boltzmann constant (5.669×10 -8 W m -2 K -4 ) and Φ, T and Ω are the phase function, the local temperature (°C) and the solid angle, respectively.

Mesh generation
In this pre-processing stage, a large computational domain should be built including the cross-sectional characteristic section of the greenhouse studied and its surroundings. To this purpose, we followed the recommendations given by Kim et al. (2017), stating that the limits of the air inlet, air outlet and upper limit should be at a minimum distances of 7H, 15H and 6H from the greenhouse, respectively, where H is the maximum greenhouse height (Figure 1). A two-dimensional simulation approach was adopted, since these models result in quick and precise solutions for transport phenomena studies, especially for greenhouses equipped with side and rooftop vents . In addition, these twodimensional approaches, when the ventilation openings are perpendicular to the flow of external air as is the case in this research, allow a solution to be obtained for the thermal distribution inside the greenhouse that is fairly close to that offered by a three-dimensional model (Benni et al., 2016;Mesmoudi et al., 2017) This model made it possible to determine the thermal behavior and air flow patterns generated in the greenhouse under night weather conditions. The greenhouse model consisted of five 11 m wide and 55 m long spans, covered with standard commercial polyethylene and with minimum heights under the gutter and under the roof of 4.56 and 7.26 m, respectively. The greenhouse included double fixed overhead vents with a 1 m opening on each span and rolling windows on all its sides with a maximum opening of 3.6 m ( Figure 2).
Once the appropriate size of the computational domain was defined, the meshing process and its respective test of independence to the size of the grid was carried out. This step is essential to obtain reliable numerical results in the shortest time possible (Tadj et al., 2017). For the grid sensitivity analysis, a total of 9 simulations were performed where the change variable corresponded to the size of the numerical grid. These simulations were carried out under the standard configuration of the greenhouse, configuration that does not include the system of closing the roof ventilation areas. Likewise, for each simulation, the initial temperature conditions in the greenhouse's external environment were defined and kept constant at 13.1°C and an external wind speed of 0.4 m s -1 . The response variable analyzed was the temperature behavior inside the greenhouse, as it was done in the study developed by He et al. (2017).
The selection of the grid is made by analyzing the numerical results of the 9 simulations developed. For this it is necessary to obtain the numerical data of the temperature behavior inside the greenhouse and to observe the variation of the average value between the different sizes of grid evaluated (Figure 3). In this case it could be observed how after the size of grid number 6 the temperature behavior is almost constant and its variation is less than 0.05 °C. Therefore, this grid size guarantees a numerical solution independent of the grid size and at the same time a solution with a lower computational cost. According to these results the computational domain was divided into an unstructured grid with a total of 648,147 square elements (Figure 4).
The selected grid presented a total of 6 element sizes. The smallest grid elements were 0.03 m long by 0.03 m high all inside the greenhouse and a grid with larger elements of 0.5 m and 0.5 m near the limits of the computational domain ( Figure 5).
The quality parameters evaluated were the cell size and the cell-to-cell size variation, finding that 98.8% of the cells in the mesh were within the high-quality range (0.95-1). Orthogonal quality was also evaluated, where the minimum value obtained was 0.93, a result that was classified within the high-quality range (Flores-Velázquez et al., 2011). A semi-implicit method was adopted for the equations linked to the pressure to solve the pressure-moment equations by means of the SIMPLE algorithm. The convergence criteria were set at 10 -6 for the equations of momentum, turbulence, energy and continuity and at 10 -8 for radiation .

Boundary conditions
The left border of the greenhouse was established as the entrance of the wind and we evaluated wind speeds below 1.5 m s -1 (0.21 -1.4 m s -1 ). This wind speed range is the predominant one in the study area for the night hours. These values were established from the climatic characterization carried out with historical data recorded during a period of fifteen years as reported in Table 1 together with the average temperature values.
The upper surface of the computational domain was established as an edge boundary with symmetrical conditions in order not to influence the air flow dynamics due to friction losses, the right and left surfaces were established as air inlet and outlet boundary limits and their velocity inlet and pressure outlet boundary conditions respectively. The cover, floor and walls of the greenhouse were established with wall type edge conditions, as well as the physical and thermal properties established in Table 2 were adopted, taken from Mesmoudi et al. (2017). The lower part of the computational domain corresponding to the greenhouse floor was modeled as an opaque solid and additionally it was established a limit condition of heat transfer to the interior of the greenhouse, procedure similar to the one established in the work developed by Piscia et al. (2012).
In addition, the optical properties of the cover material required by the DO radiation model were included, which provides a solution to the phenomenon of long wave thermal radiation. Therefore, it was considered for a transparent polyethylene are reflection coefficient (ρ = 0.11), absorption coefficient (α = 0.69) and transmission coefficient (τ = 0.19) (Montero et al., 2005).

Considered scenarios
The climatic evaluations were carried out assuming three generic meteorological conditions from which the simulated scenarios were defined from a series of combinations of meteorological conditions, night time and operation of the inflatable duct system (Table 3). The meteorological conditions ranged from clear sky to overcast nights.

Model validation
The validation methodology considered the compilation of experimental data to determine the thermal behavior of a multi-tunnel greenhouse covered with polyethylene of passive type to real scale, established in the Centro de Bio-Sistemas of the Universidad de Bogotá Jorge Tadeo Lozano (4°53'2.56 "N, 74°0'47.55" W, 2655 masl). The total area covered was 3025 m 2 resulting from 55 meters of cross section by 55 meters of longitudinal section, the cross section of the greenhouse was oriented in a North-South direction.
The greenhouse is equipped with an automated system of polyethylene ducts that are inflated with air between 16:00 and 07:00 hours the next day, allowing the closure of the upper ventilation areas. The ducts were located along the roof openings of the greenhouse and were connected internally by a duct that ran transversely throughout the greenhouse. This duct allowed the inflation of the entire system through a fan located at one end of the transver- Rev. Ceres, Viçosa, v. 67, n.5, p. 345-356, sep/oct, 2020 sal duct in one of the side walls of the greenhouse at a height of 6 meters.
For validation purposes, one of the side spans of the greenhouse was isolated from the other four by a plastic partition installed from the ground to the gutter ( Figure  6). This allows both spans to be exposed to the same condition of being laterally connected to each other, making it possible to compare the results obtained for the two configurations evaluated. The above was considered due to the logistical difficulty that prevents the experiment from being carried out under the two configurations evaluated in two contiguous greenhouses with similar characteristics.
The field data recording period used for the validation was from September 16 to September 21, 2018. During this period, the inflatable ducts system was only activated for the roof vents of the isolated span while the roof vents of the other four spans remained open. Once the simulations were completed, we performed the post-process step by  extracting the temperature and velocity data of the indoor air at the points which coincided with the internal sampling points. The experimental and numerical simulation data were compared by calculating the mean quadratic error (RMSE) and mean absolute error (MAE). (2) (3) where Vmi and Vsi, represent the observed and simulated values, respectively, and m the number of observations compared.

CFD model validation
The simulation conditions used to represent the configuration of the experimental greenhouse for validation purposes considered the atmospheric conditions of scenario C at hour 19:00, that is, an outside wind speed of 1.4 ± 0.31 m s -1 and an average temperature of 12.8 ± 0.24°C. The air movement for the evaluated configuration shows two different behaviors (Figure 7a), in the first span the convective movement generated from the soil towards the area of the roof is observed with an average wind speed of 0.2 m s -1 . For the remaining spans, a movement similar to those described above is observed for simulated cases without the presence of the duct system and the average wind speed for this condition was 1.1 m s -1 . Figure 7b shows the results of the thermal behavior for the 19:00 hour. Under these conditions an average temperature of 13.9°C was obtained for the first span. For the remaining spans the average temperature was 12.1°C, which confirms the advantages of the airtightness offered by the closure duct system.
The comparison between simulated and measured temperatures for the 19:00 hour is shown in Figure 8. In general, it is observed that for temperature the trend of the measured and simulated data is similar in the transversal axis of the evaluated greenhouse, although the simulated data show a smaller magnitude than the data obtained experimentally. The RMSE and MAE for the temperature were 0.47 and 0.44°C, results that are like those obtained in studies such as the one developed by Espinal-Montes et al. (2015). On the other hand, for wind speed it was found that the simulated and measured values are in a close range, with a MAE and RMSE value of 0.02 and 0.025 m s -1 respectively. This leads to the conclusion that the CFD model makes satisfactory predictions, therefore, the outcome of the simulations presented in this study are valid and adequately represent the thermal behavior of the studied greenhouse.  Rev. Ceres, Viçosa, v. 67, n.5, p. 345-356, sep/oct, 2020

Air flow
During night conditions different behaviors were observed depending on the activation or not of the closing system of ducts. When the ducts are in operation, a convective type movement can be observed in each of the greenhouse spans where rotating movements are formed from the heat source, which in this case is the soil, towards the cold zone represented by the roof surface (Figures 9b, d, f). This behavior is characteristic for closed greenhouses that lack active heating systems (Mesmoudi et al., 2017), because the movements of air flow are induced by free convection through the phenomenon of buoyancy that is mainly generated by the effect of the energetic loss of the soil surface inside the greenhouse.
The velocity values oscillated between 0.15 and 0.52 m s -1 , where the maximum value obtained was for the simulated case of scenario A. The atmospheric condition of the sky, as a supposed blackbody for this scenario, presented a capacity to extract energy from the interior of the greenhouse more accelerated as compared to scenarios B and C.
The scenarios where the duct system was turned off showed an air movement driven by the wind and thermal components of the natural ventilation phenomenon (Figures 9a, c, e). The streams of air flow generated from the outside to the interior of the greenhouse, through the fixed overhead vents, move through the cross section of the greenhouse to go out through the overhead vents again. The simulated values for wind speed ranged  between 0.15 and 0.42 m s -1 . These results are influenced by the thermal differentials between the interior and exterior of the greenhouse (∆T = T internal average -T external average) and the external wind speed.

Thermal behavior Scenario A -Dry clear sky night
Figures 10, c, e show the thermal behaviors for the simulation of nocturnal conditions with clear nights without using the roof closing ducts. The average internal temperatures for each case were 9.48, 7.41 and 6.59 °C for the 19:00, 00:00 and 05:00 hours, respectively. The ∆T generated under this condition for the considered hours were -3.61, -2.7 and -2.62 ºC respectively, which indicates that under these atmospheric conditions and without any type of heating system, the thermal inversion occurs, as has already been demonstrated in previous studies (Mesmoudi et al., 2012;Montero et al., 2013). This phenomenon is generated mainly by an accelerated loss of heat through the roof and the flow of heat extracted from the interior of the greenhouse through the fixed overhead openings (Teitel et al., 2008).
Figures 10b, d, f show the thermal behaviors for 19:00, 00:00 and 05:00 hours where the temperature inside the greenhouse reached average values of 11.4, 8.9 and 8.2 °C respectively, with ∆T of 1.7, -1.2 and -1.0 °C between the inside and the outside of the greenhouse. Under the atmospheric conditions of the scenario A, the ducts system does not have the total capacity to maintain the thermal conditions inside the greenhouse in values similar to the outside conditions, although it should be noted that the use of the ducts limits the heat losses by free convection flows through the roof openings and for this reason the thermal inversion is lower. The recommendation in this situation for passive greenhouses is to ventilate through the side vents, seeking to generate an adequate ventilation rate that allows to level the energetic levels of the outside air with the inside of the greenhouse (Espinal-Montes et al., 2015).

Scenario B -Humid clear sky night
The results under clear and humid night conditions indicated thermal behaviors more favorable than those found for scenario A as a consequence of the sky Figure 11: Contours of temperature distribution for the scenario B, (a) hour 19:00 without duct system, (b) hour 19:00 with duct system, (c) hour 00:00 without duct system, (d) hour 00:00 with duct system, (e) hour 05:00 without duct system and (f) hour 05:00 with duct system. Figure 10: Contours of the temperature distribution for the scenario A, (a) hour 19:00 without duct system, (b) hour 19:00 with duct system, (c) hour 00:00 without duct system, (d) hour 00:00 with duct system, (e) hour 05:00 without duct system and (f) hour 05:00 with duct system.
Rev. Ceres, Viçosa, v. 67, n.5, p. 345-356, sep/oct, 2020 temperature increase and the corresponding reduction in the energy transfer value of the interior of the greenhouse towards the external environment (Iglesias et al., 2009).
The average temperature inside the greenhouse obtained in the simulations with the use of the duct system were 13.4, 10.5 and 9.4 °C for the hours 19:00, 00:00 and 05:00, respectively (Figures 11b, d, f); with these results, the average ∆T were 0.3, 0.4 and 0.2 °C. On the other hand, the average values of the interior temperature of the greenhouse obtained in the simulations without using the duct system were 12.7, 9.9 and 8.5 °C for the hours 19:00, 00:00 and 05:00 (Figures 11a, c, e), resulting in ∆T of -0.4, -0.2 and -0.7 ºC. Under these atmospheric conditions, the use of the duct system is an alternative that allows to increase the value of the thermal tightness of the greenhouse, obtaining average greenhouse temperatures significantly higher than those of the exterior. These ∆T would be sufficient to avoid or limit the potential damage to crops under some weather conditions of the Bogota plateau where the outside temperature is close to 0 °C during certain times of the year.

Scenario C -Humid overcast night
Under this condition, which in turn is the least critical for night-time conditions, it should be mentioned that the simulated cases are the ones with the highest ∆T positive generated among all the scenarios evaluated. The thermal behaviors for the simulations without the presence of the duct system are shown in Figures 12a, c, e. The average temperatures were 14.1, 10.9 and 9.5 °C for the hours 19:00, 00:00 and 05:00, resulting in ∆T of 1.0, 0.7 and 0.3 ºC. For the simulated cases with the presence of the duct system, the average internal greenhouse temperatures were 15, 11.7 and 10.4 °C for the 19:00, 00:00 and 05:00 hours, respectively, which generated ∆T of 1.9, 1.5 and 1.2 °C (Figures 12b, e, f).
These higher ∆T values are clearly a positive effect of the duct system and their restriction to let the hot air from inside the greenhouse escape to the outside through the roof openings. For this scenario, it could be observed that the heat released through the soil was sufficient to counteract the energy losses through the roof, which allows to maintain a higher energy level inside the greenhouse.

CONCLUSIONS
According to the results, it can be assured that the use of numerical CFD models are appropriate to evaluate microclimate optimization systems in passive greenhouses under nocturnal conditions. In this investigation, the results indicated that the use of duct systems in roof vents limit or prevent the occurrence of the phenomenon of thermal inversion for atmospheric conditions such as clear and humid nights. Likewise, in atmospheric conditions with humidity and presence of cloudiness, the roof closure duct system helps to maintain a higher internal energy level compared to situations where the duct system is not present or not in operation.
For the critical atmospheric condition of cloudless sky and low humidity the duct closure system reduced the negative thermal gradient. Even so, this system is not enough to prevent the greenhouse thermal inversion, so the recommendation for passive greenhouses is to ventilate through the side vents when this condition occurs.
This study confirms that the CFD technique can continue to be used as an alternative to evaluate night-time conditions applied to other greenhouse designs or other climate control alternatives for this condition.

CONFLICT OF INTEREST
The authors declare that there is no conflict of interesting carrying this research and publishing this manuscript. Figure 12: Contours of temperature distribution for the scenario C, (a) hour 19:00 without duct system, (b) hour 19:00 with duct system, (c) hour 00:00 without duct system, (d) hour 00:00 with duct system, (e) hour 05:00 without duct system and (f) hour 05:00 with duct system