3.3.1. Side Ventilation (S)

The distribution of the thermal behavior for each of the evaluated scenarios can be seen in Figure 11. In general terms, irrespective of the outside wind speed, it was found that the cool zones correspond to the ventilation regions where the air enters the greenhouse, while the regions of higher temperature are located in the region where the air flow exits from the interior of the structure. This behavior occurs mainly because the air that enters the greenhouse through the windward window of span 1, as it crosses the cross section, mixes with the warm air and increases its energy level by heat transfer [10,45].

On the other hand, the height of the greenhouse directly influences the magnitude and spatial distribution of the temperature; it was observed that as the height of the greenhouse increases, the magnitude of the temperature and the percentage of the area of the structure with high temperature values decrease. This is an effect generated by the higher renovation index and the higher level of thermal inertia obtained in higher greenhouses [1,35]. Likewise, it is possible to observe the effects of the external wind speed, where for the highest speed scenario S4 it was found that the temperature value inside the greenhouse was lower and also presented a greater homogeneous behavior, coinciding with what was reported Flores-Velázquez [81]. *Sustainability* **2021**, *13*, x FOR PEER REVIEW 17 of 29

**Figure 11.** Simulated thermal distribution patterns with ventilation configuration S and for the six greenhouse models (M1 to M6) under wind speeds S1 (0.2 ms−1) and S4 (1.6 ms−1). **Figure 11.** Simulated thermal distribution patterns with ventilation configuration S and for the six greenhouse models (M1 to M6) under wind speeds S1 (0.2 ms−<sup>1</sup> ) and S4 (1.6 ms−<sup>1</sup> ).

In numerical terms, the average air temperature inside the greenhouse ranged between maximum values of 37.5 ± 8.9 °C for S-M1S1 and a minimum of 23.3 ± 0.6 °C for S-M6S4. For the S1 scenario it can be observed that the temperature decreases 11.7 °C comparing greenhouses M1 and M6 respectively, likewise for the S4 scenario this temperature reduction between M1 and M6 is only 2.5 °C, therefore, the increase in greenhouse height In numerical terms, the average air temperature inside the greenhouse ranged between maximum values of 37.5 ± 8.9 ◦C for S-M1S1 and a minimum of 23.3 ± 0.6 ◦C for S-M6S4. For the S1 scenario it can be observed that the temperature decreases 11.7 ◦C comparing greenhouses M1 and M6 respectively, likewise for the S4 scenario this temperature reduction between M1 and M6 is only 2.5 ◦C, therefore, the increase in greenhouse height

will be more relevant in regions where calms or low wind speeds predominate. Although under these conditions it should also be analyzed up to what level it is convenient from

M4 and M5 models under this condition no less important reductions of 9.7 and 10.8 °C,

Another relevant criterion to be analyzed is the uniform distribution of temperature in the cross axis of the greenhouse. This is undoubtedly a factor that has received more attention in recent years, since it has a direct influence on the physiological behavior of the plants and on the final yield of the crops [58,82]. Therefore, for this study, the values obtained in the cross section of each of the greenhouses were plotted as a function of wind

In general, it is again observed the difference that exists between span 1 on the windward side (cool area) and span 9 on the leeward side (warm area) with relevant thermal differentials higher than 15 °C. Likewise, in some critical scenarios, the temperature in some areas of the cross section of the greenhouse exceeds 40 °C, a value that is quite inad-

equate for the growth and development of any vegetal species [83].

respectively, would be obtained.

speed (Figure 12).

will be more relevant in regions where calms or low wind speeds predominate. Although under these conditions it should also be analyzed up to what level it is convenient from the technical and economic point of view to increase the greenhouse height, since for the M4 and M5 models under this condition no less important reductions of 9.7 and 10.8 ◦C, respectively, would be obtained.

Another relevant criterion to be analyzed is the uniform distribution of temperature in the cross axis of the greenhouse. This is undoubtedly a factor that has received more attention in recent years, since it has a direct influence on the physiological behavior of the plants and on the final yield of the crops [58,82]. Therefore, for this study, the values obtained in the cross section of each of the greenhouses were plotted as a function of wind speed (Figure 12).

In general, it is again observed the difference that exists between span 1 on the windward side (cool area) and span 9 on the leeward side (warm area) with relevant thermal differentials higher than 15 ◦C. Likewise, in some critical scenarios, the temperature in some areas of the cross section of the greenhouse exceeds 40 ◦C, a value that is quite inadequate for the growth and development of any vegetal species [83]. *Sustainability* **2021**, *13*, x FOR PEER REVIEW 18 of 29

**Figure 12.** Spatial distribution of air temperature inside the greenhouse for the simulated scenarios under the S configuration.

#### under the S configuration. 3.3.2. Roof Ventilation (R)

already described.

3.3.2. Roof Ventilation (R) In qualitative terms, it is observed that under ventilation configuration R, lower temperatures are obtained with respect to configuration S (Figure 13). In general, it is ob-In qualitative terms, it is observed that under ventilation configuration R, lower temperatures are obtained with respect to configuration S (Figure 13). In general, it is observed that the highest temperature areas are located on the windward and leeward sides of the span, respectively. This behavior is related to the ventilation configuration and to the characteristics of the airflows of these regions near the greenhouse walls previously

**Figure 12.** Spatial distribution of air temperature inside the greenhouse for the simulated scenarios

served that the highest temperature areas are located on the windward and leeward sides of the span, respectively. This behavior is related to the ventilation configuration and to the characteristics of the airflows of these regions near the greenhouse walls previously

For this configuration it is also observed that there is a positive effect of wind speed

and greenhouse height on the magnitude and distribution of the temperature inside the structure. In the most critical case S-M1S1 it is observed that there are relatively important areas of high temperature in spans 1, 3 and 5, while in the S-M6S1 scenario these areas are significantly reduced. In S-M1S4 the high temperature area disappears in span 3 and the hot areas in spans 1 and 6 are reduced to an area very close to the side wall of each side and, finally, in S-M6S4 these high temperature areas disappear completely in the spans

of multi-tunnel greenhouses by Park et al. [84].

analyzed, a similar behavior was reported in a work on natural ventilation of four types of multi-tunnel greenhouses by Park et al. [84].

For this configuration it is also observed that there is a positive effect of wind speed and greenhouse height on the magnitude and distribution of the temperature inside the structure. In the most critical case S-M1S1 it is observed that there are relatively important areas of high temperature in spans 1, 3 and 5, while in the S-M6S1 scenario these areas are significantly reduced. In S-M1S4 the high temperature area disappears in span 3 and the hot areas in spans 1 and 6 are reduced to an area very close to the side wall of each side and, finally, in S-M6S4 these high temperature areas disappear completely in the spans already described. *Sustainability* **2021**, *13*, x FOR PEER REVIEW 19 of 29

**Figure 13.** Simulated thermal distribution patterns with ventilation configuration R and for the six greenhouse models (M1 to M6) under wind speeds S1 (0.2 ms−1) and S4 (1.6 ms−1). **Figure 13.** Simulated thermal distribution patterns with ventilation configuration R and for the six greenhouse models (M1 to M6) under wind speeds S1 (0.2 ms−<sup>1</sup> ) and S4 (1.6 ms−<sup>1</sup> ).

In numerical terms the average air temperature ranged between maximum value of 32.7 ± 3.1 °C for R-M1S1, value that is lower by 4.8 °C with respect to S-M1S1, the minimum value obtained was 25.0 ± 1.3 °C in R-M6S4 value that is higher by 1.7 °C with respect to S-M6S4. For the S1 scenario the temperature reduction obtained was 1.3, 2.2, 2.8, 3.6 and 4.3 °C in M2 to M6 compared to M1 respectively, while in S4 these reductions with respect to M1 were 0.5, 0.7, 1.1, 1.1, 1.4, 1.7 °C in M2 to M6, successively, results that continue to confirm the importance of roof ventilation in low wind speed conditions (S1). In numerical terms the average air temperature ranged between maximum value of 32.7 ± 3.1 ◦C for R-M1S1, value that is lower by 4.8 ◦C with respect to S-M1S1, the minimum value obtained was 25.0 ± 1.3 ◦C in R-M6S4 value that is higher by 1.7 ◦C with respect to S-M6S4. For the S1 scenario the temperature reduction obtained was 1.3, 2.2, 2.8, 3.6 and 4.3 ◦C in M2 to M6 compared to M1 respectively, while in S4 these reductions with respect to M1 were 0.5, 0.7, 1.1, 1.1, 1.4, 1.7 ◦C in M2 to M6, successively, results that continue to confirm the importance of roof ventilation in low wind speed conditions (S1).

For this case, the thermal behavior on the greenhouse cross axis in each scenario evaluated shows a totally different spatial temperature distribution than that observed in the previous ventilation configuration (Figure 14). Due to the air inlet and outlet flows through the roof vents and the airflow distribution patterns discussed in Figure 6, it is possible to observe the temperature variations occurring between span 1, 3 and 6 succes-For this case, the thermal behavior on the greenhouse cross axis in each scenario evaluated shows a totally different spatial temperature distribution than that observed in the previous ventilation configuration (Figure 14). Due to the air inlet and outlet flows through the roof vents and the airflow distribution patterns discussed in Figure 6, it is possible to observe the temperature variations occurring between span 1, 3 and 6 successively. It can

sively. It can also be seen how this temperature distribution tends to be more uniform and smaller in magnitude as the height of the greenhouse and the outside wind speed increase.

also be seen how this temperature distribution tends to be more uniform and smaller in magnitude as the height of the greenhouse and the outside wind speed increase. *Sustainability* **2021**, *13*, x FOR PEER REVIEW 20 of 29

**Figure 14.** Spatial distribution of air temperature inside the greenhouse for the simulated scenarios **Figure 14.** Spatial distribution of air temperature inside the greenhouse for the simulated scenarios under the R configuration.

#### under the R configuration. 3.3.3. Side Ventilation and Roof (RS)

3.3.3. Side Ventilation and Roof (RS) The spatial distribution of the temperature for this configuration allows us to observe qualitatively that the temperature value in each of the simulated scenarios has a lower magnitude compared to the S and R configurations (Figure 15). For these scenarios, it is The spatial distribution of the temperature for this configuration allows us to observe qualitatively that the temperature value in each of the simulated scenarios has a lower magnitude compared to the S and R configurations (Figure 15). For these scenarios, it is observed that the regions of lower temperature coincide with the ventilation areas where the air flow enters, both in the lateral sides and in the roof ventilation areas, which coincides with what was analyzed in the work by Villagran et al. [10].

observed that the regions of lower temperature coincide with the ventilation areas where the air flow enters, both in the lateral sides and in the roof ventilation areas, which coincides with what was analyzed in the work by Villagran et al. [10]. For scenario S1, it can be observed that the high temperature regions are in the spans located on the leeward side. This high temperature region gradually disappears as the height of the greenhouse increases. In the case of M1, the highest temperature region occurs in spans 4, 5 and 6, with a heat spot in the middle of the area of span 6. However, in M6, this heat patch disappears from the span 6 and only a small region of higher temperature is observed near the ventilation area on the leeward side of the greenhouse. Similar behavior is observed in scenario S4, although in this case the temperature values appear to be qualitatively lower than those obtained in S1.

Numerically, it was found that the temperature value inside the structure presented a maximum value of 29.4 ± 3.8 ◦C in RS-M1S1, being this value lower by 8.1 and 3.3 ◦C compared to S and R, respectively. Therefore, the RS ventilation configuration also allows obtaining a higher cooling efficiency under the same climatic conditions, which is a very relevant factor for the management of the microclimate in greenhouses located in the savanna of Bogota [1].

On the other hand, the minimum mean inside air temperature value was 23.2 ± 0.6 ◦C in RS-M6S4, a value that is lower by 1.7 and 0.1 ◦C with respect to those obtained in R and S, respectively. Therefore, again it can be observed that the greatest benefits in terms of thermal distribution are obtained for low wind speed conditions. It is also important to note that the temperature values obtained in most of the scenarios are within the range of 25 and 30 ◦C, ranges where species of commercial and food interest such as Tomatoes tend to develop adequately. As well as some ornamental species of commercial interest for the international market such as the Rose or Carnation [1,85].

Regarding the spatial distribution, it was again found that as the greenhouse height increases and the external wind speed increases, more stable and uniform temperature values are obtained in the cross section of the greenhouse (Figure 16). The thermal gradients for the same moment inside the greenhouse were reduced from 13.2 ◦C in M1S1 to 1.94 ◦C in M6S4, the latter value being a recommended limit to guarantee the homogeneity of a naturally ventilated greenhouse [39,86]. *Sustainability* **2021**, *13*, x FOR PEER REVIEW 21 of 29

**Figure 15.** Simulated thermal distribution patterns with the RS ventilation configuration and for the six greenhouse models (M1 to M6) under wind speeds S1 (0.2 ms−1) and S4 (1.6 ms−1). **Figure 15.** Simulated thermal distribution patterns with the RS ventilation configuration and for the six greenhouse models (M1 to M6) under wind speeds S1 (0.2 ms−<sup>1</sup> ) and S4 (1.6 ms−<sup>1</sup> ).

For scenario S1, it can be observed that the high temperature regions are in the spans located on the leeward side. This high temperature region gradually disappears as the height of the greenhouse increases. In the case of M1, the highest temperature region occurs in spans 4, 5 and 6, with a heat spot in the middle of the area of span 6. However, in M6, this heat patch disappears from the span 6 and only a small region of higher temperature is observed near the ventilation area on the leeward side of the greenhouse. Similar behavior is observed in scenario S4, although in this case the temperature values appear

compared to S and R, respectively. Therefore, the RS ventilation configuration also allows obtaining a higher cooling efficiency under the same climatic conditions, which is a very relevant factor for the management of the microclimate in greenhouses located in the sa-

On the other hand, the minimum mean inside air temperature value was 23.2 ± 0.6 °C in RS-M6S4, a value that is lower by 1.7 and 0.1 °C with respect to those obtained in R and S, respectively. Therefore, again it can be observed that the greatest benefits in terms of thermal distribution are obtained for low wind speed conditions. It is also important to note that the temperature values obtained in most of the scenarios are within the range of 25 and 30 °C, ranges where species of commercial and food interest such as Tomatoes tend to develop adequately. As well as some ornamental species of commercial interest for the

international market such as the Rose or Carnation [1,85].

to be qualitatively lower than those obtained in S1.

vanna of Bogota [1].

of a naturally ventilated greenhouse [39,86].

**Figure 16.** Spatial distribution of air temperature inside the greenhouse for the scenarios simulated **Figure 16.** Spatial distribution of air temperature inside the greenhouse for the scenarios simulated under the RS configuration.

#### under the RS configuration. *3.4. Response Surface Modeling*

*3.4. Response Surface Modeling* Finally, looking for a better and simpler interpretation of the results obtained by CFD simulation, these numerical results of the evaluated scenarios were adjusted to response Finally, looking for a better and simpler interpretation of the results obtained by CFD simulation, these numerical results of the evaluated scenarios were adjusted to response surface models. This methodology allows observing the behavior of the response variable in scenarios not simulated numerically, if these scenarios are within the limiting ranges of the simulated boundary conditions

Regarding the spatial distribution, it was again found that as the greenhouse height increases and the external wind speed increases, more stable and uniform temperature values are obtained in the cross section of the greenhouse (Figure 16). The thermal gradients for the same moment inside the greenhouse were reduced from 13.2 °C in M1S1 to 1.94 °C in M6S4, the latter value being a recommended limit to guarantee the homogeneity

surface models. This methodology allows observing the behavior of the response variable in scenarios not simulated numerically, if these scenarios are within the limiting ranges of the simulated boundary conditions The results for the fitted response surface models are summarized in Table 4. The second order model alternative gave the best results for all models. The lack of fit test showed that the models accurately fit the data exception made for the sidewalls scenario and with temperature as the response variable. The adjusted R-squared coefficient of de-The results for the fitted response surface models are summarized in Table 4. The second order model alternative gave the best results for all models. The lack of fit test showed that the models accurately fit the data exception made for the sidewalls scenario and with temperature as the response variable. The adjusted R-squared coefficient of determination indicated varying results whether the outcome variable was the internal air velocity or the temperature. Highest R-squared values were obtained for the models with wind speed as the dependent variable as compared to their temperature pairs under the same scenario.

termination indicated varying results whether the outcome variable was the internal air velocity or the temperature. Highest R-squared values were obtained for the models with wind speed as the dependent variable as compared to their temperature pairs under the same scenario. Despite the varying results for the adjusted R-squared coefficient, all terms included in the models were considered significant. Since most of the results indicated a good fit of the model to the calibration data, we evaluated the effect of varying levels of greenhouse height and external wind speed for each ventilation scenario through response surface plots. The response surface plots indicated a similar behavior for internal temperature and air velocity irrespective of the ventilation scenario (Figure 17). The individual effects of the predictors showed that increasing the greenhouse height effectively resulted in lower temperatures, while the internal air velocity was slightly affected by the greenhouse height. The major effect of the external wind speed on both, the temperature, and the internal air

velocity, is clearly depicted in Figure 17. Increasing the greenhouse height under increasing levels of external wind speed resulted in lower temperature values within the greenhouse.

**Table 4.** Response surface models and goodness of fit measures applied to the three ventilation scenarios and considering internal wind speed and temperature as dependent variables.


*IAV*: internal air velocity; *T*: temperature; *H*1: greenhouse height; *EWS*: external wind speed.

On the other hand, lower greenhouse heights combined with decreased external wind speeds resulted in lower internal air velocity. However, while all scenarios depicted the same trends, differences arise when looking at the scale of variation of the internal temperatures and air velocity. The sidewalls scenario showed that increasing the greenhouse height under increased external wind speeds have a more dramatic effect in terms of lowering the internal temperature of the greenhouse (Figure 17c), compared to roof or combined ventilation (Figure 17a,e).

The option to apply only the roof ventilation decreased the internal wind speed under almost null external wind speed no matter what the greenhouse height is. However, when the external wind speed increases along with the greenhouse height, internal air velocity also increases but the highest internal air velocity is reached with the lowest greenhouse height and the maximum external wind speed (Figure 17b). This results is opposed to the other scenarios where the top internal air velocity was reached when the maximum greenhouse height and external wind speed were evaluated (Figure 17d,f).

The optimal combination of greenhouse height and external wind speed that minimized the temperature and achieved an adequate internal wind speed seemed to be located around the maximum height and external wind speed for all scenarios. These results implied that higher values for both explanatory variables should be considered to find the lowest internal temperature and the highest internal air velocity. However, going beyond and increasing the region of analysis will yield unrealistic results due to greenhouse construction restrictions and maximum local wind speeds.

The results presented here are in line with those reported by Villagran et al. [1], and Bustamante et al. [87], indicating that increased greenhouse heights and the combination of roof and sidewalls ventilation yields improved internal climate conditions, particularly for greenhouse established in high altitude tropical regions.

It is also important to recommend that a future study regarding this greenhouse typology should focus on the structural design of the model evaluated and an economic analysis of the height alternatives proposed. This to be able to dimension the robustness of the structure and an associated cost per m<sup>2</sup> of covered area for each of the 6 models proposed, although in Colombia there are currently no regulations in force for the construction of greenhouses, a work of this magnitude can provide economic and structural tools for decision making. This information can be used for the generation of public policy management documents that can promote and regulate the use of greenhouses conceived with design criteria.

**Figure 17.** Effect of greenhouse height and external wind speed on internal air velocity and temperature for the roof (**a**,**b**), sidewalls (**c**,**d**), and roof + sidewalls (**e**,**f**) ventilation scenarios. **Figure 17.** Effect of greenhouse height and external wind speed on internal air velocity and temperature for the roof (**a**,**b**), sidewalls (**c**,**d**), and roof + sidewalls (**e**,**f**) ventilation scenarios.

## **4. Conclusions**

The numerical simulation tool is still an agile and accurate alternative to determine the renewal indexes and thermal distribution in passive greenhouses. It is also analyze unconstructed scenarios, facilitating decision making for farmers and providing an alternative analysis that can contribute to improve the technical and economic sustainability of protected agriculture.

The use of previously validated CFD models allows generating simulated data series of temperature and air velocity inside a greenhouse and associating them to response variables such as the height of the structure evaluated or others. This facilitated the association of the obtained data to other analysis methodologies such as regression with response surface models obtaining acceptable adjustments in the analyzed variables, therefore, it was possible to generate response surface graphs which facilitates the interpretation of the numerical results as a whole for each ventilation configuration analyzed (S, R, RS).

The natural ventilation of a greenhouse is highly dependent on the behavior of the wind speed of the study region. For the three ventilation configurations analyzed (S, R and RS), it was found that the renewal indexes were maximized as a function of the increase in the external wind speed. Although it should be noted that regardless of the ventilation configuration for the lowest wind speed S1 (0.2 ms−<sup>1</sup> ), none of the greenhouse models exceeded the recommended minimum value of 40 hourly renovations, although these IR values are less critical in the RC configuration. The optimization of the IR in the RS ventilation configuration allowed obtaining a higher cooling efficiency inside the greenhouse and a uniformity in the spatial distribution of the temperature. Maximum temperature values of 29.4 ± 3.8 ◦C were obtained in RS-M1S1, being this value lower by 8.1 and 3.3 ◦C with respect to S-M1S1 and R-M1S1 respectively.

The crop production in greenhouse in the future should be with the efficient use of resources, so that numerical simulation techniques will have the goal of adapting the climatic environment in the design and operation of protected agriculture. In addition to the greenhouse dimensions, the position and size of windows are among the other influential factors to be considered in future studies

**Author Contributions:** Conceptualization, E.V., M.A., J.F.-V.; Methodology, E.V., C.B., J.F.-V., M.A.; Software, E.V.; Validation, E.V., C.B.; Investigation, E.V., C.B.; Writing—original draft preparation, E.V., C.B., J.F.-V., M.A.; Writing—review and editing, E.V., C.B., J.F.-V., M.A. All authors have read and agreed to the published version of the manuscript.

**Funding:** The present study was funded by the Servicio Nacional de Aprendizaje (SENA), la Asociación Colombiana de Exportadores de Flores (Asocolflores) y al Centro de Innovación de la Floricultura Colombiana (Ceniflores) through the project "Generación de una herramienta de diseño u optimización de ventilación natural de los invernaderos dedicados a la producción de flores de corte en cuatro subregiones de la Sabana de Bogotá, mediante el uso de herramientas de simulación basadas en la técnica Dinámica de Fluidos Computacional (CFD)".

**Data Availability Statement:** Not applicable.

**Acknowledgments:** To the Universidad Jorge Tadeo Lozano for the support for the development of the doctoral research project called "*Development of a methodological proposal for the design, evaluation and microclimatic optimization of passive greenhouses through the use of numerical simulation tools*".

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

