**1. Introduction**

In Colombia, plasticulture and protected agriculture have promoted a method of agricultural production that generates higher crop yields per unit of available land area, compared to open field production [1]. The use of passive greenhouses, where the microclimate is managed by means of natural ventilation, has become more widespread. This type of greenhouse is also quite common in other regions of the world [2]. This cultivation method allows, among others, the intensification of horticultural or ornamental production [3], better management of water resources and fertilizer applications [4,5], partial or total control of the micro-climatic variables that affect crop growth and development [6]. In recent years, this production technology has also become a crop protection tool used by producers in the face of increasingly frequent and severe weather events due to climate change [3,7,8].

The use of greenhouses or protected agriculture structures with plastic roofing worldwide has been estimated at more than 3.5 million hectares, which have been established

**Citation:** Villagrán, E.;

Flores-Velazquez, J.; Akrami, M.; Bojacá, C. Influence of the Height in a Colombian Multi-Tunnel Greenhouse on Natural Ventilation and Thermal Behavior: Modeling Approach. *Sustainability* **2021**, *13*, 13631. https://doi.org/10.3390/ su132413631

Academic Editors: Muhammad Sultan, Yuguang Zhou, Walter Den and Uzair Sajjad

Received: 25 October 2021 Accepted: 6 December 2021 Published: 9 December 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

mainly in countries such as China, Korea, Spain, Mexico, France and Italy [9]. In the specific case of Colombia, statistics show that there are approximately 11,000 hectares destined for the ornamental and horticultural sectors. The main production areas are located mostly in regions with altitudes above 1500 m above sea level, where cold and humid mountain climate conditions predominate [10–12]. In Colombia, as in the other countries mentioned above, different types of greenhouses have been built and a high percentage of them are light structures with a low technological level [13,14].

These typologies of geometric designs of plastic roofing greenhouses already conceived, originated from some initial concerns of the producers such as land area, easily available structural and roofing materials, versatile structure designs for its construction, but mainly the final economic cost of the greenhouse or structure used. [15,16]. However, an adequate design process should include an analysis of the local climatic conditions so that, based on these conditions, the builders or decision makers can propose the appropriate architecture and geometry of the greenhouse that, together with the selected roofing material, will generate adequate microclimate conditions for the growth and development of the crops [17,18].

In terms of climatic information, there are already local or regional information tools or systems where it is possible to obtain data series of the main meteorological variables of a specific site. In the last three decades, significant progress has been made in improving the optical and thermal properties of polyethylene roofing material, seeking greater benefits for the development of the main horticultural and ornamental crops [19,20]. Therefore, in each country the producer has the possibility of selecting a type of plastic roofing according to his needs. This facilitates the selection of an optimal roofing material or the most suitable for the local climate conditions, so that inside the structures the behavior of temperature, relative humidity and vapor pressure deficit can be managed within the optimal range for agricultural production [21].

On the other hand, there is the architectural aspect and the geometrical dimensions of the greenhouse to be built. Although there are some design parameters established for naturally ventilated greenhouses, such as the width and length [22]. The information available on other parameters such as height is variable and many of the related studies recommend increasing the height of structures used in tropical and subtropical climate conditions [10,23]. However, the height level of the structure is not defined to a specific limit but varies between types of greenhouses and types of climatic conditions [17,24]. It would be ideal, therefore, to have at the availability of the greenhouse structure designers, a software or a design methodology that allows them to evaluate this characteristic and to be able to define the adequate dimensions. This design tool should be able to be implemented in different types of greenhouses and under different operation scenarios.

One of the methodologies that can offer quick solutions is computational fluid dynamics (CFD) simulation. This simulation tool has been implemented to analyze aspects related to the characteristics of multiple types of greenhouses [25–27]. In addition, when a CFD model is successfully validated experimentally, it is possible to analyze unbuilt scenarios, which promotes the efficient use of resources and avoids the construction of unsuitable greenhouse structures or undesirable microclimate conditions [28,29].

Regarding the use of CFD studies applied to the optimization of passive greenhouse structures, it should be mentioned that, for example, in a two-dimensional study for a Chinese solar greenhouse, it was determined that, as longer greenhouse sections were generated, higher temperatures were produced inside the structure [30]. Similarly, Villagrán et al. [10] determined that by increasing the roof ventilation surface with respect to the covered floor area (SVC/SC) from 2.5 to 20.1% in a traditional Colombian greenhouse, it was possible to obtain ventilation rates higher than 0.04 m<sup>3</sup> s <sup>−</sup><sup>1</sup> m−<sup>2</sup> for the prevailing conditions of the local climate.

Regarding the greenhouse width, it is known that in Almeria type structures [31], which include the arc type [32] and gothic tunnel type [33], as the more spans that are joined laterally, the greater the width of the greenhouse, the ventilation rate is reduced exponentially, which in turn generates higher temperature values inside the greenhouse [31–33]. For the height dimension, Boulard and Fatnassi [34] using a CFD-2D study on an 8-span arch-type greenhouse found that when the height above the greenhouse gutter increases from 3 to 5 m, the reduction of the thermal gradient between the interior and the exterior is reduced from 3.5 to 2.0 ◦C. This allowed the authors to conclude that the greenhouse height positively affects the natural ventilation phenomenon and allows to generate better thermal conditions inside the greenhouse [34]. joined laterally, the greater the width of the greenhouse, the ventilation rate is reduced exponentially, which in turn generates higher temperature values inside the greenhouse [31–33]. For the height dimension, Boulard and Fatnassi [34] using a CFD-2D study on an 8-span arch-type greenhouse found that when the height above the greenhouse gutter increases from 3 to 5 m, the reduction of the thermal gradient between the interior and the exterior is reduced from 3.5 to 2.0 °C. This allowed the authors to conclude that the greenhouse height positively affects the natural ventilation phenomenon and allows to generate better thermal conditions inside the greenhouse [34].

Regarding the greenhouse width, it is known that in Almeria type structures [31], which include the arc type [32] and gothic tunnel type [33], as the more spans that are

*Sustainability* **2021**, *13*, x FOR PEER REVIEW 3 of 29

Recently Fatnassi et al. [35] conducted a study in a tunnel type greenhouse with butterfly roof vents, where the thermal behavior of the greenhouse was evaluated numerically under four different gutter heights, 4, 5, 6 and 7 m, respectively. It was found that when the height is increased from 4 to 5 m, the temperature inside the greenhouse decreases by approximately 2 ◦C, while when it is increased from 6 to 7 m the temperature value only decreases by 0.2 ◦C. According to these results, the authors concluded that the effect of increasing greenhouse height on air temperature is clearly asymptotic, therefore, the increase in height should be analyzed as a climatic benefit versus economic cost ratio. Recently Fatnassi et al. [35] conducted a study in a tunnel type greenhouse with butterfly roof vents, where the thermal behavior of the greenhouse was evaluated numerically under four different gutter heights, 4, 5, 6 and 7 m, respectively. It was found that when the height is increased from 4 to 5 m, the temperature inside the greenhouse decreases by approximately 2 °C, while when it is increased from 6 to 7 m the temperature value only decreases by 0.2 °C. According to these results, the authors concluded that the effect of increasing greenhouse height on air temperature is clearly asymptotic, therefore, the increase in height should be analyzed as a climatic benefit versus economic cost ratio.

In this study, an experimentally validated CFD-2D model has been implemented under the climatic conditions of a region of the savanna of Bogota, Colombia. The CFD model was used to evaluate the effects on the air flow velocity and the interior temperature of a gothic-type greenhouse as the height level increases from 2.5 m to 5.0 m, under three ventilation configurations, side ventilation (R), roof vents (S) and combined ventilation (RS). Finally, the numerical data obtained in each of the developed simulation scenarios were grouped by ventilation configuration and fitted to response surface regression models. This was done to generate a more integrated and simpler method of analysis of the response variables analyzed and their relationship with the microclimate generated In this study, an experimentally validated CFD-2D model has been implemented under the climatic conditions of a region of the savanna of Bogota, Colombia. The CFD model was used to evaluate the effects on the air flow velocity and the interior temperature of a gothic-type greenhouse as the height level increases from 2.5 m to 5.0 m, under three ventilation configurations, side ventilation (R), roof vents (S) and combined ventilation (RS). Finally, the numerical data obtained in each of the developed simulation scenarios were grouped by ventilation configuration and fitted to response surface regression models. This was done to generate a more integrated and simpler method of analysis of the response variables analyzed and their relationship with the microclimate generated

#### **2. Materials and Methods 2. Materials and Methods**

#### *2.1. Description of Prototype Proposed for Analysis 2.1. Description of Prototype Proposed for Analysis*

The analysis proposed in this research was carried out on a gothic multi-tunnel greenhouse with a 200 µm thick polyethylene cover located in the Savannah of Bogota, Colombia. The greenhouse was composed of a total of six spans, each span had a width length of 9.33 m, for a total greenhouse width length of 55.98 m (Figure 1). The greenhouse had a gutter height of 4.0 m and a maximum roof height of 8.3 m. It was equipped with ventilation areas on the side walls, with an effective opening of 3.7 m, which equals a lateral ventilation surface with respect to the covered floor area of 13.21%. Likewise, the ventilation surface was complemented with a roof ventilation area of 1.5 m of total opening in each of the spans, therefore, the roof ventilation surface with respect to the covered surface is 16.1%. The analysis proposed in this research was carried out on a gothic multi-tunnel greenhouse with a 200 µm thick polyethylene cover located in the Savannah of Bogota, Colombia. The greenhouse was composed of a total of six spans, each span had a width length of 9.33 m, for a total greenhouse width length of 55.98 m (Figure 1). The greenhouse had a gutter height of 4.0 m and a maximum roof height of 8.3 m. It was equipped with ventilation areas on the side walls, with an effective opening of 3.7 m, which equals a lateral ventilation surface with respect to the covered floor area of 13.21%. Likewise, the ventilation surface was complemented with a roof ventilation area of 1.5 m of total opening in each of the spans, therefore, the roof ventilation surface with respect to the covered surface is 16.1%.

**Figure 1.** Geometric diagram of the cross section of the greenhouse studied.

#### *2.2. Computational Fluid Dynamics Modeling*

The numerical calculation of the physical processes involving fluid movement and related variables such as pressure, force, density, temperature and their associated changes in heat and mass transfer phenomena can be solved using computational fluid dynamics [36]. Currently, CFD is one of the most widely used analysis, design or redesign methodologies in different sectors of the economy, therefore, it is possible to find CFD studies approached from mechanical, environmental, aeronautical and automotive engineering [37]. On the other hand, its use in the agricultural sector has not been insignificant and in the last two decades it has been implemented as a tool for process analysis involving natural ventilation processes or energy optimization of greenhouse structures [22,38].

CFD analysis is a resolution technique that solves a set of nonlinear partial equations from algebraic discretization by means of numerical simulation using the finite volume method. A CFD simulation process is composed of three phases of development, preprocessing, processing and post-processing [39]. In the preprocessing phase, the computational domain, and the geometry of the structure on which the natural ventilation process will be analyzed are designed. In this phase the numerical meshing process applied to the whole computational domain including the greenhouse geometry is also generated, these processes must be performed based on criteria of good CFD simulation practices [22,40].

In the processing phase, the numerical solution of the evaluated problem is performed, when the air and energy flow in the model is calculated by solving the governing equations based on the physical laws of conservation of energy, mass and momentum. In this phase, the sub models considered as source terms of the governing equations are selected, such as turbulence, buoyancy or free convection, solar radiation, porous media and phenomena associated with mass transfer [41]. Finally, in the post-processing phase, it is possible to develop a qualitative and quantitative analysis of the numerical results obtained for each scenario evaluated [42].

It should be noted that the CFD-2D model used in this research has been previously successfully validated experimentally, verifying its high predictive capacity for air flow patterns and thermal distribution inside the greenhouse analyzed. The validation results can be verified in the works developed by Villagrán and Bojacá [16,43] and by Villagrán et al. [1]. Therefore, for this research work, details related to the validation process will not be discussed, although relevant aspects of each of the CFD simulation phases will be mentioned.

#### 2.2.1. Pre-Processing

The commercial software ANSYS ICEM CFD (v. 18.0, ANSYS Inc., Canonsburg, PA, USA, EE. UU.) was used, by means of which a large computational domain was created in a two-dimensional configuration. It included the geometry of the cross section of the analyzed greenhouse, as well as the process of meshing the computational domain was performed by means of this software. The two-dimensional CFD studies of natural ventilation are useful and have a capacity to predict the ventilation rate and the thermal distribution, inside a structure when the dominant external wind currents blow perpendicular to the ventilation areas [27,44].

It was determined that the dimensions of the computational domain, establishing as the reference for the maximum height of the structure (H), should have dimensions from the region of the airflow inlet to the windward sidewall of the structure of 15H. From the leeward sidewall to the airflow outlet boundary of 20H and a minimum height from ground level of 10H (Figure 2). These dimensions are established following the guidelines given in numerical studies of natural ventilation of greenhouses and are the dimensions that allow to adequately describe the physical phenomena associated with the analyzed problem within the atmospheric boundary layer [1,45].

*Sustainability* **2021**, *13*, x FOR PEER REVIEW 5 of 29

**Figure 2.** Overall dimensions of the computational domain. **Figure 2.** Overall dimensions of the computational domain. **Figure 2.** Overall dimensions of the computational domain.

In the meshing process, the computational domain volume was discretized into an unstructured grid of square or rectangular elements (Figure 3). The size of the numerical grid was defined by a mesh independence test as described in the work done by Villagran et al. [1], at the end of the analysis process, the selected numerical grid was composed of a total of 1,104,369 elements. The quality parameters evaluated were cell size and cell-tocell size variation; it was found that 95.2% of the cells of the mesh were within the highquality interval (0.95–1). The orthogonality criterion was also evaluated, where the minimum value obtained was 0.92, a value considered as high quality [10]. In the meshing process, the computational domain volume was discretized into an unstructured grid of square or rectangular elements (Figure 3). The size of the numerical grid was defined by a mesh independence test as described in the work done by Villagran et al. [1], at the end of the analysis process, the selected numerical grid was composed of a total of 1,104,369 elements. The quality parameters evaluated were cell size and cell-to-cell size variation; it was found that 95.2% of the cells of the mesh were within the high-quality interval (0.95–1). The orthogonality criterion was also evaluated, where the minimum value obtained was 0.92, a value considered as high quality [10]. In the meshing process, the computational domain volume was discretized into an unstructured grid of square or rectangular elements (Figure 3). The size of the numerical grid was defined by a mesh independence test as described in the work done by Villagran et al. [1], at the end of the analysis process, the selected numerical grid was composed of a total of 1,104,369 elements. The quality parameters evaluated were cell size and cell-tocell size variation; it was found that 95.2% of the cells of the mesh were within the highquality interval (0.95–1). The orthogonality criterion was also evaluated, where the minimum value obtained was 0.92, a value considered as high quality [10].

**Figure 3.** Detail of the cross-sectional grid of the greenhouse evaluated. **Figure 3.** Detail of the cross-sectional grid of the greenhouse evaluated.

**Figure 3.** Detail of the cross-sectional grid of the greenhouse evaluated. Once the geometry has been constructed, it is possible to define the boundary conditions to be established in the limits of the computational domain and in the geometry of the greenhouse (Figure 1). In this specific case the right side was defined as the air flow inlet boundary, for which a condition determined by a logarithmic wind speed profile was established, this logarithmic profile depends on the climatic characteristics and the type of soil in the study region, factors that have been previously determined in the work Once the geometry has been constructed, it is possible to define the boundary conditions to be established in the limits of the computational domain and in the geometry of the greenhouse (Figure 1). In this specific case the right side was defined as the air flow inlet boundary, for which a condition determined by a logarithmic wind speed profile was established, this logarithmic profile depends on the climatic characteristics and the type of soil in the study region, factors that have been previously determined in the work Once the geometry has been constructed, it is possible to define the boundary conditions to be established in the limits of the computational domain and in the geometry of the greenhouse (Figure 1). In this specific case the right side was defined as the air flow inlet boundary, for which a condition determined by a logarithmic wind speed profile was established, this logarithmic profile depends on the climatic characteristics and the type of soil in the study region, factors that have been previously determined in the work developed by Villagrán et al. [1].

developed by Villagrán et al. [1]. The left boundary on the other hand was determined as a pressure and airflow outflow boundary, while the upper boundary of the computational domain was set as a wall boundary condition with an imposed solar radiation flux. The wall-type boundary condition was also established for the lower boundary of the computational domain, the floor, the walls and the roof of the greenhouse, while the ventilation areas were set to the indoor boundary condition depending on the ventilation configuration analyzed. The physical and optical properties of the materials within the computational domain were also dedeveloped by Villagrán et al. [1]. The left boundary on the other hand was determined as a pressure and airflow outflow boundary, while the upper boundary of the computational domain was set as a wall boundary condition with an imposed solar radiation flux. The wall-type boundary condition was also established for the lower boundary of the computational domain, the floor, the walls and the roof of the greenhouse, while the ventilation areas were set to the indoor boundary condition depending on the ventilation configuration analyzed. The physical and optical properties of the materials within the computational domain were also defined with the values summarized in Table 1. The left boundary on the other hand was determined as a pressure and airflow outflow boundary, while the upper boundary of the computational domain was set as a wall boundary condition with an imposed solar radiation flux. The wall-type boundary condition was also established for the lower boundary of the computational domain, the floor, the walls and the roof of the greenhouse, while the ventilation areas were set to the indoor boundary condition depending on the ventilation configuration analyzed. The physical and optical properties of the materials within the computational domain were also defined with the values summarized in Table 1.

fined with the values summarized in Table 1.


**Table 1.** Physical and optical properties used in the CFD simulation. *Sustainability* **2021**, *13*, x FOR PEER REVIEW 6 of 29

Subsequently, six greenhouse models (M1–M6) were created with the only variations being the height under the gutter (H1) and the lateral ventilation dimension (SV), as described in Table 2. The meshing process for models M1 to M6 was performed on the same geometric model already discretized from model M4, simply decreasing the greenhouse height for models M1 to M3 and increasing it for models M5 and M6. Emissivity 0.95 0.9 0.7 Subsequently, six greenhouse models (M1–M6) were created with the only variations being the height under the gutter (H1) and the lateral ventilation dimension (SV), as described in Table 2. The meshing process for models M1 to M6 was performed on the same geometric model already discretized from model M4, simply decreasing the greenhouse height for models M1 to M3 and increasing it for models M5 and M6.

**Table 2.** Geometric details of the analyzed greenhouse models. **Table 2.** Geometric details of the analyzed greenhouse models.


the ventilation configurations used locally, which are: ventilation by the lateral sides (S), by the roof vents (R) and combined ventilation (RS). After this, the greenhouse model to be evaluated was defined from M1 to M6 and finally the climatic conditions to be used as starting conditions for each simulation were determined. These conditions were known from previous studies developed in the same study region and where it was established that the maximum daytime temperature averaged 22.3 °C, the maximum solar radiation level was 853 Wm-2 and wind conditions can vary between values of 0.2 and 1.6 ms−1 . Therefore, it was determined to perform the evaluations under four wind speeds S1 (0.2 ms−1), S2 (0.5 ms−1), S3 (1 ms−1) and S4 (1.6 ms−1). According to the above, 72 simulations were performed and coded as shown in Table 3. **Table 3.** Coding of the 72 simulated scenarios. **Scenarios** Finally, the simulation scenarios to be analyzed were defined, these were built from the ventilation configurations used locally, which are: ventilation by the lateral sides (S), by the roof vents (R) and combined ventilation (RS). After this, the greenhouse model to be evaluated was defined from M1 to M6 and finally the climatic conditions to be used as starting conditions for each simulation were determined. These conditions were known from previous studies developed in the same study region and where it was established that the maximum daytime temperature averaged 22.3 ◦C, the maximum solar radiation level was 853 Wm-2 and wind conditions can vary between values of 0.2 and 1.6 ms−<sup>1</sup> . Therefore, it was determined to perform the evaluations under four wind speeds S1 (0.2 ms−<sup>1</sup> ), S2 (0.5 ms−<sup>1</sup> ), S3 (1 ms−<sup>1</sup> ) and S4 (1.6 ms−<sup>1</sup> ). According to the above, 72 simulations were performed and coded as shown in Table 3.

Finally, the simulation scenarios to be analyzed were defined, these were built from

S-M1S1 S-M4S1 R-M1S1 R-M4S1 RS-M1S1 RS-M4S1



#### 2.2.2. Processing

For the numerical solution of the Reynolds-averaged Navier-Stokes (RANS) equations, the ANSYS FLUENT processing software was used. The air flow in the greenhouse was considered to be turbulent and with a density change associated with temperature changes that can be modeled by the Boussinesq approximation. The method of resolution selected for the CFD model was semi-implicit for the pressure and velocity bound equations using the SIMPLE algorithm. Finally, residual convergence criteria were established in 10−<sup>6</sup> for the energy equation and in 10−<sup>3</sup> for the other variables such as continuity, turbulence and momentum [41].

The influence of solar radiation intensity on the spatial behavior of temperature was considered using the discrete order (DO) radiation model, which allows modeling the radiation and calculating the convective exchanges that occur in the computational domain, [46].

It was also decided not to include any type of crop, because the purpose of the study is not to analyze the heat and mass transfer flows that occur between the indoor environment and any species of plants. On the contrary, the main objective is to make an analysis of the ventilation phenomenon under the most critical condition that can occur in a real scenario and that is in empty conditions and under the most extreme climatic conditions [37,47]. This is an approach that is still valid and continues to be implemented in a significant number of studies carried out in various countries [48–51]. In addition, as different types of horticultural, ornamental, aromatic and medicinal crops are grown in the region, the analysis developed in this study can be applied to any type of crop.

#### 2.2.3. Post-Processing

For this phase, where the main objective is to perform the qualitative and quantitative analysis of the evaluated scenarios, ANSYS CFD-Post software was used. Therefore, for each of the 72 simulations developed, two-dimensional plots of temperature and airflow patterns were generated and numerical values of air velocity and temperature in the greenhouse cross section at a height of 1.4 m above ground level were extracted. These data sets were analyzed for each simulation scenario and as a whole using response surface regression models.

#### 2.2.4. Response Surface Modeling

Response surface regression models were fitted to each of the three ventilation scenarios under consideration (roof (R), side (S), and roof—side (RS)). For each scenario, models were fitted to evaluate the response of varying levels of exterior wind speed and greenhouse height on the internal temperature and air velocity. The objectives of applying response surface models on top of the CFD simulated scenarios was threefold: (1) to establish an approximate relationship between the response variables and the exterior wind speed and greenhouse height that can be used for predicting temperature or internal wind speed values for a given combination of the predictors; (2) to determine the statistical significance of the explanatory variables through hypothesis testing; and (3) to determine the optimum combination of the explanatory variables that result in the maximum response over the region under consideration.

Response surface models are an extension of linear models and works similarly; however, extra arguments should be added to consider the response-surface portion of the model [52]. In the present study, second order models including first-order response, two-way interactions and quadratic terms for the explanatory variables were included in the model's formulation. Once calibrated, the models were tested for their goodness of fit through measures such as the coefficient of determination, lack of fit and pure error [53,54].

The output of the models was plotted as perspective plots to depict the combined effect of the experimental factors (external wind speed and greenhouse height) over the dependent variables (internal temperature and air velocity). The models were fitted using the rsm package (v. 2.10.2; [52]) included in the statistical software R (v. 4.0.4; [55]).

#### **3. Results and Discussion**

*3.1. Airflow Patterns*

3.1.1. Side Ventilation (S)

The airflow patterns under the ventilation configuration through the sidewall (S) areas of the greenhouse can be seen in Figure 4. In general terms, it is observed that the airflow pattern enters the greenhouse through the ventilation area on the windward side and moves horizontally in an airflow pattern that has a height very similar to the height under the gutter of each of the greenhouses evaluated (M1-M6), this horizontal flow leaves the greenhouse again through the ventilation area arranged on the leeward wall. *Sustainability* **2021**, *13*, x FOR PEER REVIEW 9 of 29

**Figure 4.** Air distribution patterns simulated with ventilation configuration S and for the six greenhouse models (M1 to M6) under wind speeds S1 (0.2 ms<sup>−</sup><sup>1</sup> ) and S4 (1.6 ms<sup>−</sup><sup>1</sup> ). **Figure 4.** Air distribution patterns simulated 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> ).

Another factor that showed a direct relationship with height was the uniformity of the air flow velocity behavior along the cross axis of the greenhouse (Figure 5). In general terms, it is observed that as the height of the structure increases, the velocity behavior curves show less oscillation between the points of high and low air velocity. This can also be verified with the standard deviation (sd) values. In S-M1S4 the sd value was ± 0.27 ms−<sup>1</sup> while for S-M6S4 it was ± 0.16 ms−1. This will undoubtedly promote a more homogeneous microclimate behavior inside the greenhouse, helping to limit the negative impacts on growth and development generated by non-uniform microclimates [58]. Another common feature is the air circulation loops that are generated in the upper part of each of the greenhouse spans, this occurs because due to the buoyancy phenomenon and the associated pressure components, the warm air inside the greenhouse moves vertically towards the roof area and its space in the lower part of each span is occupied by the fresh air that enters from the outside environment [3,56]. It should also be noted that qualitatively it is observed that as the height of the greenhouse increases, the air flow that moves over the nearby area where the crops are grown has a higher velocity; this same behavior is observed in the air flow that leaves the greenhouse through the leeward ventilation.

To quantitatively analyze the airflows, velocity data were extracted from a cross profile at a height of 1.4 m above ground level. In general, the airflows show a behavior for velocity that oscillates between a minimum of 0.25 <sup>±</sup> 0.05 ms−<sup>1</sup> for S-M1S1 and a maximum of 1.42 <sup>±</sup> 0.16 ms−<sup>1</sup> for S-M6S4. It was also observed that the air flows inside the greenhouse show a direct relationship with the wind speed outside, therefore, the higher the wind speed, the higher the air velocity inside the greenhouse. This occurs in greenhouses where this type of ventilation is analyzed and where porous insect-proof screens are not contemplated in the ventilation areas [29,57].

On the other hand, it is also observed that as the height of the greenhouse increases, the average air flow velocity inside the greenhouse increases for the wind speed scenarios analyzed (S1 to S4). In summary, the increase in velocity between the extreme simulation scenarios was 42.3, 20.9, 19.4 and 16.2% for M1S1, M1S2, M1S3, M1S4 compared to M6S1, M6S2, M6S3, M6S4, respectively. These velocity increases in the higher greenhouse models occur because the effect of air friction is reduced as the ratio (L/H1) between the greenhouse length (L) and the height above the gutter (H1) decreases, which coincides with what was previously reported in the study developed by Chu et al. [44].

Another factor that showed a direct relationship with height was the uniformity of the air flow velocity behavior along the cross axis of the greenhouse (Figure 5). In general terms, it is observed that as the height of the structure increases, the velocity behavior curves show less oscillation between the points of high and low air velocity. This can also be verified with the standard deviation (sd) values. In S-M1S4 the sd value was <sup>±</sup>0.27 ms−<sup>1</sup> while for S-M6S4 it was <sup>±</sup> 0.16 ms−<sup>1</sup> . This will undoubtedly promote a more homogeneous microclimate behavior inside the greenhouse, helping to limit the negative impacts on growth and development generated by non-uniform microclimates [58]. *Sustainability* **2021**, *13*, x FOR PEER REVIEW 10 of 29

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

#### the S configuration. 3.1.2. Roof Ventilation (R)

roof of this span.

spans 4 and 5.

3.1.2. Roof Ventilation (R) Figure 6 shows the simulated airflow patterns for this ventilation configuration. In this case, since the side windows are closed, the airflow patterns are forced in and out through the windows in the roof region [16]. For the low velocity scenarios S1 where the Figure 6 shows the simulated airflow patterns for this ventilation configuration. In this case, since the side windows are closed, the airflow patterns are forced in and out through the windows in the roof region [16]. For the low velocity scenarios S1 where the thermal component of natural ventilation dominates, a behavior where three flow patterns are generated is observed. In span 1, an air flow was identified that enters through the

thermal component of natural ventilation dominates, a behavior where three flow patterns are generated is observed. In span 1, an air flow was identified that enters through the ventilation area and becomes a convective movement loop between the floor and the

Subsequently, between span 2 and the middle of span 3 there is an airflow pattern that moves in the direction of the outside wind at ground level and a flow that moves in the opposite direction to the outside wind at the top of the roof of these spans. Finally, through the ventilation area of span 3 an airflow enters and mixes with another airflow entering from the ventilation area of span 4, flows that move towards the leeward wall and exits through the ventilation areas of spans 4 and 5. It should also be noted that the confluence of flows over the central area of span 3 generates a small low velocity loop near the floor region, being this an inadequate movement pattern that may cause the genera-

For the case of the S4 scenario wind speed, the natural ventilation of the greenhouse will depend on the thermal and wind components together [59,60]. In this case, we were able to identify an air flow pattern that enters through span 1 and forms a recirculation pattern between the roof of the span, the floor and the wall of the windward side, where the highest air flow velocity occurs at ground level, which coincides with that reported by Kwon et al. [61]. Likewise, part of the air flow entering through the window of span 1 is mixed with another air flow entering through the window of span 2, which generates an acceleration of the air flow above ground level; the same behavior is repeated with the air flows entering through spans 3 and 4; finally, the air flows leave the greenhouse through

tion of a heat spot over this region [29,33,35].

ventilation area and becomes a convective movement loop between the floor and the roof of this span.

Subsequently, between span 2 and the middle of span 3 there is an airflow pattern that moves in the direction of the outside wind at ground level and a flow that moves in the opposite direction to the outside wind at the top of the roof of these spans. Finally, through the ventilation area of span 3 an airflow enters and mixes with another airflow entering from the ventilation area of span 4, flows that move towards the leeward wall and exits through the ventilation areas of spans 4 and 5. It should also be noted that the confluence of flows over the central area of span 3 generates a small low velocity loop near the floor region, being this an inadequate movement pattern that may cause the generation of a heat spot over this region [29,33,35].

For the case of the S4 scenario wind speed, the natural ventilation of the greenhouse will depend on the thermal and wind components together [59,60]. In this case, we were able to identify an air flow pattern that enters through span 1 and forms a recirculation pattern between the roof of the span, the floor and the wall of the windward side, where the highest air flow velocity occurs at ground level, which coincides with that reported by Kwon et al. [61]. Likewise, part of the air flow entering through the window of span 1 is mixed with another air flow entering through the window of span 2, which generates an acceleration of the air flow above ground level; the same behavior is repeated with the air flows entering through spans 3 and 4; finally, the air flows leave the greenhouse through spans 4 and 5. *Sustainability* **2021**, *13*, x FOR PEER REVIEW 11 of 29

**Figure 6.** Simulated air distribution patterns with ventilation configuration R and for the six greenhouse models (M1 to M6) under wind speeds S1 (0.2 ms<sup>−</sup><sup>1</sup> ) and S4 (1.6 ms<sup>−</sup><sup>1</sup> ). **Figure 6.** Simulated air 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 air flow velocity inside the greenhouse over the area where the crops are grown ranged from a minimum of 0.26 ± 0.07 ms−1 for R-M1S1 to a maximum of 0.92 ± 0.36 ms−1 in R-M1S4. In this specific case no increase in air flow velocity is observed as the greenhouse height increases in the case of the low velocity scenario S1 it is In numerical terms, the air flow velocity inside the greenhouse over the area where the crops are grown ranged from a minimum of 0.26 <sup>±</sup> 0.07 ms−<sup>1</sup> for R-M1S1 to a maximum of 0.92 <sup>±</sup> 0.36 ms−<sup>1</sup> in R-M1S4. In this specific case no increase in air flow velocity is observed

observed that between M1 and M6 there is only an increase of 11.5% but this same value

the S3 scenarios the air flow velocity decreased up to 20.5% compared M1 with M5 and M6 and finally in S4 the air velocity also decreased up to 15% compared M1 with M6.

Regarding the behavior of the air flow in the cross axis of the greenhouse, it was found that for models M1, M2 and M3 there is a very similar distributed behavior in space, while for M4, M5 and M6 there is a small change in this spatial distribution, which may be caused by the change in greenhouse heights (Figure 7). Likewise, the inside air velocities increase as the wind speed increases, which has already been demonstrated in several

is obtained between M1 compared to M2, M3 and M4.

research works [22,33,49,62].

as the greenhouse height increases in the case of the low velocity scenario S1 it is observed that between M1 and M6 there is only an increase of 11.5% but this same value is obtained between M1 compared to M2, M3 and M4.

For the case of S2 the air velocity increased between M1 and M2 by 2.5% but compared to the other greenhouse models it decreased by up to 12.5% with M6 specifically. In the S3 scenarios the air flow velocity decreased up to 20.5% compared M1 with M5 and M6 and finally in S4 the air velocity also decreased up to 15% compared M1 with M6.

Regarding the behavior of the air flow in the cross axis of the greenhouse, it was found that for models M1, M2 and M3 there is a very similar distributed behavior in space, while for M4, M5 and M6 there is a small change in this spatial distribution, which may be caused by the change in greenhouse heights (Figure 7). Likewise, the inside air velocities increase as the wind speed increases, which has already been demonstrated in several research works [22,33,49,62]. *Sustainability* **2021**, *13*, x FOR PEER REVIEW 12 of 29

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

the R configuration. On the other hand, there is a quite heterogeneous behavior among each of the six spans, finding that the smallest standard deviation among these scenarios was ± 0.07 ms−1 in the low velocity scenarios S1, while for the high velocity scenarios S4 this standard deviation was up to ± 0.36 ms−1. The lowest air speed points were the regions near the On the other hand, there is a quite heterogeneous behavior among each of the six spans, finding that the smallest standard deviation among these scenarios was <sup>±</sup>0.07 ms−<sup>1</sup> in the low velocity scenarios S1, while for the high velocity scenarios S4 this standard deviation was up to <sup>±</sup> 0.36 ms−<sup>1</sup> . The lowest air speed points were the regions near the leeward and windward wall and the highest air speed points were the areas below spans 2 and 4.

#### leeward and windward wall and the highest air speed points were the areas below spans 3.1.3. Side and Roof Ventilation (RS)

2 and 4. 3.1.3. Side and Roof Ventilation (RS) The spatial distribution of the airflow patterns for the combined ventilation configuration can be seen in Figure 8. Qualitatively, more continuous and intense airflows can be The spatial distribution of the airflow patterns for the combined ventilation configuration can be seen in Figure 8. Qualitatively, more continuous and intense airflows can be observed over the entire cross-sectional area analyzed, both in the crop and canopy regions, which should positively impact the greenhouse renovation rates and directly the thermal performance of the greenhouse [62].

observed over the entire cross-sectional area analyzed, both in the crop and canopy regions, which should positively impact the greenhouse renovation rates and directly the

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

Numerically, the air flows presented a mean velocity that ranged from a minimum value of 0.31 ± 0.09 ms−1 in RS-M1S1 to a maximum of 1.53 ± 0.17 in RS-M6S4. For the case of S1, it is observed that the airflow velocity increases up to 29.1% in M5 and M6 compared to M1. Likewise, in S4 the increase in velocity was 15% in the M6 greenhouse with respect to the value obtained in M1, therefore, the greenhouse height has a significant impact on the airflow velocity. This is mainly due to the fact that, according to the analysis of natural ventilation carried out in previous studies, the greater the distance between the central axis of the side ventilation with respect to the central axis of the roof ventilation, the more dynamic the wind and thermal component is and therefore the better the air flow patterns are obtained [56,63]. On the other hand, the spatial behavior of the air flow velocity inside the greenhouse Numerically, the air flows presented a mean velocity that ranged from a minimum value of 0.31 <sup>±</sup> 0.09 ms−<sup>1</sup> in RS-M1S1 to a maximum of 1.53 ± 0.17 in RS-M6S4. For the case of S1, it is observed that the airflow velocity increases up to 29.1% in M5 and M6 compared to M1. Likewise, in S4 the increase in velocity was 15% in the M6 greenhouse with respect to the value obtained in M1, therefore, the greenhouse height has a significant impact on the airflow velocity. This is mainly due to the fact that, according to the analysis of natural ventilation carried out in previous studies, the greater the distance between the central axis of the side ventilation with respect to the central axis of the roof ventilation, the more dynamic the wind and thermal component is and therefore the better the air flow patterns are obtained [56,63].

cross section can be found for some of the simulated scenarios in Figure 9. In general terms, it can be observed that as the height of the greenhouse increases, the less oscillation between the points of higher and lower velocity exists in the spatial distribution curves. At the same time, it should be noted again that these changes in air flow velocity are mainly due to the dynamics of incoming and outgoing air flows and to the interaction of warm and cold air masses that mix and generate buoyancy flows [10,64]. Although these values are less critical in the M3 to M5 greenhouse scenarios under the RC combined ventilation configuration. On the other hand, the spatial behavior of the air flow velocity inside the greenhouse cross section can be found for some of the simulated scenarios in Figure 9. In general terms, it can be observed that as the height of the greenhouse increases, the less oscillation between the points of higher and lower velocity exists in the spatial distribution curves. At the same time, it should be noted again that these changes in air flow velocity are mainly due to the dynamics of incoming and outgoing air flows and to the interaction of warm and cold air masses that mix and generate buoyancy flows [10,64]. Although these values are less critical in the M3 to M5 greenhouse scenarios under the RC combined ventilation configuration.

Regarding air flow velocities, in general, under the three configurations evaluated, the values obtained are within the ranges reported in passive greenhouse ventilation studies [5,65]. Although for the low velocity condition S1, 100% of the airflows present velocities lower than the minimum recommended value of 0.5 ms−<sup>1</sup> for plant growth inside greenhouses [42,65].

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

#### the RS configuration. *3.2. Renewal Index (RI)*

Regarding air flow velocities, in general, under the three configurations evaluated, the values obtained are within the ranges reported in passive greenhouse ventilation studies [5,65]. Although for the low velocity condition S1, 100% of the airflows present velocities lower than the minimum recommended value of 0.5 ms−1 for plant growth inside greenhouses [42,65]. The renewal index (RI) was calculated using the method of integration of the volumetric flow of air leaving the ventilation areas of the greenhouse, the results obtained can be seen in Figure 10. The RI values ranged from a minimum of 5.89 Vol h−<sup>1</sup> obtained in S-M1S1 to a maximum of 72.3 Vol h−<sup>1</sup> obtained in S-M6S4. It is important to highlight that these are the contrasting scenarios, therefore, the minimum RI was obtained in the greenhouse with the lowest height and ventilation area, evaluated at the lowest wind speed (S1), while the maximum RI was obtained in the greenhouse with the highest height and ventilation area evaluated at the highest wind speed (S4).

*3.2. Renewal Index (RI)*  The renewal index (RI) was calculated using the method of integration of the volumetric flow of air leaving the ventilation areas of the greenhouse, the results obtained can be seen in Figure 10. The RI values ranged from a minimum of 5.89 Vol h−1 obtained in S-M1S1 to a maximum of 72.3 Vol h−1 obtained in S-M6S4. It is important to highlight that these are the contrasting scenarios, therefore, the minimum RI was obtained in the greenhouse with the lowest height and ventilation area, evaluated at the lowest wind speed For the lowest wind speed (S1), RI values between 5.89 Vol h−<sup>1</sup> and 19.1 Vol h−<sup>1</sup> were obtained in the side ventilation configuration S for M1 and M6, respectively. Likewise, for the ventilation configuration through the vents of the roof region R, RI values of 19.4 and 19.7 Vol h−<sup>1</sup> were obtained in M1 and M6, which are values compared to those obtained in S of 329% and 3.68%, respectively. Finally, for the RS combined ventilation configuration, an RI value was obtained for M1 of 26.5 Vol h−<sup>1</sup> which is a value 449% higher compared to the S ventilation configuration and 36.6% higher compared to the R ventilation configuration. For M6, a value of 36.8 Vol h−<sup>1</sup> was obtained, which is 92.7% and 86.8% higher than the S and R ventilation configurations, respectively.

(S1), while the maximum RI was obtained in the greenhouse with the highest height and ventilation area evaluated at the highest wind speed (S4). For the lowest wind speed (S1), RI values between 5.89 Vol h−1 and 19.1 Vol h−1 were obtained in the side ventilation configuration S for M1 and M6, respectively. Likewise, for the ventilation configuration through the vents of the roof region R, RI values of 19.4 and Regarding the S4 velocity scenario, RI values were obtained for the S ventilation configuration of 19.5 and 50.2 Vol h−<sup>1</sup> for M1 and M6, respectively. For the R ventilation configuration, values of 43.8 and 46.7 Vol h−<sup>1</sup> were obtained for M1 and M6, respectively, which represents an increase in M1 of 224.6% with respect to the S configuration; on the other hand, for M6 there was a reduction of 14.6% in RI with respect to S. In the case of the RS configuration, a value of 59.7 Vol h−<sup>1</sup> was obtained in M1, which is equivalent to a

19.7 Vol h−1 were obtained in M1 and M6, which are values compared to those obtained in S of 329% and 3.68%, respectively. Finally, for the RS combined ventilation configura-

pared to the S ventilation configuration and 36.6% higher compared to the R ventilation configuration. For M6, a value of 36.8 Vol h−1 was obtained, which is 92.7% and 86.8%

Regarding the S4 velocity scenario, RI values were obtained for the S ventilation configuration of 19.5 and 50.2 Vol h−1 for M1 and M6, respectively. For the R ventilation configuration, values of 43.8 and 46.7 Vol h−1 were obtained for M1 and M6, respectively,

higher than the S and R ventilation configurations, respectively.

higher value of RI by 306.1% with respect to the S configuration and 136.6% with respect to R. While in M6 an RI value of 72.3 Vol h−<sup>1</sup> was obtained, this value being 144 and 154.8% higher than those obtained in S and R, respectively.

To highlight the scenarios where RI values above or equal to the recommended minimum (RI <sup>≥</sup> 40 Vol h−<sup>1</sup> ) are achieved for naturally ventilated agricultural structures [11,66]. Under ventilation configuration S, this was only obtained under greenhouse models M4, M5 and M6 and under external wind speed conditions of 1.5 ms−<sup>1</sup> , with this same speed condition for ventilation configuration R adequate RI values are obtained in all greenhouse models (M1–M6).

While for the RS combined ventilation configuration, the RI values are adequate for all greenhouses (M1–M6) under wind speed scenarios higher than 1 ms−<sup>1</sup> , the same is true for the M4, M5 and M6 models under wind speeds higher than 0.5 ms−<sup>1</sup> . It is therefore relevant for greenhouse growers or greenhouse builders or decision makers to seriously analyze the prevailing wind speed conditions in the study region [62,67].

These results reaffirm some conclusions already determined in previous works related to the natural ventilation of greenhouses, where it was identified that the renewal indexes are dependent on the ventilation configuration used [68–71]. It was also identified that there is a linear relationship between the renewal index and wind speed [38,56,72,73], the ventilation configuration that allows to obtain the highest renewal index values is the combined configuration of roof and side vents [10,16,74,75].

On the other hand, in the side ventilation configuration, the increase of renewal indexes as a function of the increase of the side ventilation area is relevant and significant only in narrow greenhouses (width < 60 m) or with few attached spans (<6 spans) [76–79]. Finally, the IR in low external wind speed conditions are more stable and higher in greenhouse structures with relevant ventilation areas in the roof region [14,29,80]. *Sustainability* **2021**, *13*, x FOR PEER REVIEW 16 of 29

**Figure 10.** Calculated values of the simulated renewal indexes for each of the scenarios considered. **Figure 10.** Calculated values of the simulated renewal indexes for each of the scenarios considered.

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

*3.3. Spatial Distribution of Temperature*

what was reported Flores-Velázquez [81].

3.3.1. Side Ventilation (S)
