Simulation of the Thermal and Aerodynamic Behavior of an Established Screenhouse under Warm Tropical Climate Conditions: A Numerical Approach

Simulation of the Thermal and Aerodynamic Behavior of an Established Screenhouse under Warm Tropical Climate Conditions: A Numerical Approach

Edwin VillagranRoberto Ramirez Andrea Rodriguez Rommel Leon Pacheco Jorge Jaramillo 

Centro de Investigación Tibaitata, Corporación Colombiana de Investigación Agropecuaria - AGROSAVIA, Mosquera - Cundinamarca 250040, Colombia

Estación Experimental Enrique Jiménez Núñez, Instituto Nacional de Innovación y Transferencia en Tecnología Agropecuaria de Costa Rica – INTA., Cañas – Guanacaste 50601, Costa Rica

Centro de Investigación Caribia, Corporación Colombiana de Investigación Agropecuaria - AGROSAVIA, Sevilla – Magdalena 478020, Colombia

Centro de Investigación La Selva, Corporación Colombiana de Investigación Agropecuaria - AGROSAVIA, Rionegro - Antioquia 054040, Colombia

Corresponding Author Email:
15 January 2020
19 April 2020
18 June 2020
| Citation



In tropical countries agriculture protected with passive and low-cost structures is one of the main alternatives for intensifying agricultural production in a sustainable manner. This type of greenhouses has adequate efficiency in cold weather conditions meanwhile its use in hot weather conditions presents disadvantages due to the generation of an inadequate microclimate for the growth and development of certain species. This has generated an important interest for the use of screen houses (SH) for the horticultural and fruit production, and currently there are many studies on the behavior of microclimates in SH; however, these experiments were developed for climatic conditions in other latitudes. In this research, a study was developed using a computational fluid dynamics (CFD) 3D numerical simulation, with the aim of evaluating the thermal and aerodynamic behavior of an SH under two specific configurations (under rain (RC) and under dry conditions (DC)). The CFD model was validated by taking experimental temperature data inside the SH. The results showed that: i) the CFD model has an acceptable capacity to predict the behavior of temperature and air flows, ii) simulations can be performed under environmental conditions of day and night, and iii) the RC configuration affected the positive thermal behavior, which limited the presence of the thermal inversion phenomenon under nocturnal conditions, meanwhile under RC daytime conditions, it reduced the velocity of the air flows generating higher thermal gradients compared to DC.


computational fluid dynamics, temperature, screenhouse microclimate, wind speed

1. Introduction

Screen houses (SH) are a technological option offered by protected agriculture as an intermediate alternative between open field and greenhouse cultivation. With the implementation of these structures, the aim is to transform the land use from extensive to intensive or promote agricultural production in alternative and sustainable systems in order to generate the supply necessary to meet the demand for high-quality food throughout the year [1]. This type of structure is built on metal columns and support cables where a roof and side walls are installed, generally made of porous screens that are insect proof or shaded [2].

The adoption of this type of technology has generated a great boom since the end of the 90s and is currently a relevant component of farming systems undercover, which has gradually extended from the countries of the Mediterranean coast to regions in other latitudes, mainly with temperate or warm climates [3] and for different cultivation types and methods [4]. Commercially there is a great variety of screens that differ in material types, color and porosity. These characteristics affect their optical and aerodynamic properties; therefore, they have been strongly studied and modified seeking to improve the microclimatic conditions generated inside the SH [5, 6].

According to the manufacturing material of the porous screen used and its properties, various agricultural benefit objectives are sought such as (i) Shading for regions where solar radiation is excessive and with supra-optimal values [7]; (ii) reducing the vulnerability of crops to damage by weather events such as icy hail and wind gusts [8, 9]; (iii) cooling limitation in night-time conditions through the reduction in energy loss by radiation [10]; (iv) exclusion of insects and vectors that transmit viruses, allowing significant reductions in the application of pesticides [6, 11]; and (v) increase the efficient use of water, extending the growth period of the plants and delaying the ripening process of some horticultural products [12, 13]. In addition to the benefits mentioned above, this type of structure has become popular and widespread among farmers because they can potentially maximize the benefit of crops with a low-cost technological contribution compared to conventional greenhouses [9].

The knowledge of the microclimate in SH as well as in plastic greenhouses is essential to achieve adequate crop management [6]. The effects of different types of screens on the microclimate of plants have been studied since the beginning of the century [1, 11, 14, 15]. The use of screens mainly influences the radiation exchange and air flow dynamics, reducing its speed and modifying its turbulence characteristics [16], thus, affecting ventilation rates and heat exchanges, mass, and gases between the plants and their surrounding atmosphere. This usually translates into behaviors with high values of variables such as temperature and humidity that can cause physiological and environmental disorders conducive to the appearance of fungal diseases that affect the final crop yield [17].

The studies dedicated to the measurement, modeling and simulation of the microclimate distribution in conventional greenhouses have been extensive in the last three decades, obtaining results that have allowed to describe the distribution of temperature, humidity, CO2 concentration and the characteristics of airflow patterns, and develop management strategies to optimize the behavior of these variables [18-20]. On the other hand, the studies related in this field with SH are still scarce, although there are significant advances as summarized in the study developed by Tanny et al. [6]. Currently, there is a need to generate relevant information that allows researchers and farmers of horticultural products to obtain a deep understanding of the patterns and characteristics of the airflow in order to obtain a better design and positioning of screen houses [21] or study the aerodynamic effect of different types of screens on physical and biological processes in these systems [9].

One of the most used tools since the beginning of the century to characterize the microclimate distributed inside greenhouses and its interaction with the plant has been computational fluid dynamics (CFD). This tool models and simulates fluid flow and transfer of heat, mass, and momentum, obtaining great advances in the design and optimization of agricultural structures [22, 23]. The study of the microclimate in screenhouses can be approached through CFD numerical simulations, considering the roof material as a porous medium, which will allow evaluating a great variety of structures, screens and climatic environments in a relatively short period of time. Bartzanas et al. [1] developed a two-dimensional CFD study to assess the effect of a screen on radiation distribution, finding that the optical and spectral properties directly affect the distribution of solar radiation, and the degree of porosity of the screen reduces air velocity, affecting the thermal behavior inside the screenhouse. Other relevant studies using 3D CFD modeling were in charge of evaluating the behavior of air flows and the value of temperature in screenhouses used for tomato cultivation, reporting that these parameters are strongly affected by the degree of porosity of the screen [24]. Although these works have not been developed for the warm climate conditions of the Central American Caribbean

According to the above, the objective of this work was to determine through 3D CFD simulation the thermal and airflow patterns behavior of an insect proof screen house established in Guanacaste - Costa Rica. With the purpose of evaluating two configurations of the productive system used in different times of the year.

2. Materials and Methods

2.1 Experimental site and climatic conditions

The study area is in the coastal area of the canton of Abangares, province of Guanacaste in northwest Costa Rica (10º11' N, 85º10' W at an altitude of 10 m a.s.l.). This region has a warm tropical climate with a dry season, and according to the Köppen-Geiger climate classification, the area has an Aw climate [25]. The average multi-year average temperature is 27.7℃, with maximum and minimum averages of 36.9 and 21.1℃ (Figure 1a). The annual rainfall reaches a value of 1669.7 mm, distributed during the months of May to November (Figure 1a). The average wind speed oscillates in the year between 0.2 and 1.4 ms-1 (Figure 1b), with predominant directions between SE-SSE.

2.2 Description of the screenhouse

The development of the experimental study was carried out in a flat roof SH with a covered floor area of 1,496 m2, where the longitudinal section was in an east-west direction (E-W). The geometric characteristics of the structure were the following: width (X = 34 m), length (Z = 44 m) and height (Y = 5 m) (Figure 2a). The side walls and roof were covered with a porous insect-proof screen (Dimensions thread 16.1x10.2 and porosity ε = 0.33). Inside the screenhouse, small semicircular tunnels of 2.2 m of height and 1.2 m of width were built located along the longitudinal axis of the SH and on top of the cultivation beds, these tunnels were covered with polyethylene to be used during the rainy season, in order to avoid or reduce to the maximum the wetting of the foliage (Figure 2b).

Figure 1. Meteorological characteristics for the canton of Abangares, province of Guanacaste in northwest Costa Rica

Figure 2. Dimensions and interior detail of the screenhouse

2.3 Fundamental equations and physical models

The models explained in this section, describe the physical principles that govern the study analyzed, the models selected are those reported in the literature used for problems similar to this research and that have shown appropriate computational performance and numerical results adjusted to real behavior. The governing flow equations presented in Eq. (1). They represent as diffusion-convection equations of fluids for three conservation laws, which include the transport, impulse and energy equations of a compressible fluid and in a three-dimensional (3D) field in a steady state.

$\nabla(\rho \phi \vec{v})=\nabla(\Gamma \nabla \phi)+\mathrm{S} \phi$          (1)

where, ρ is the density of the fluid (kgm-3), ∇ is the nabla operator, ϕ represents the concentration of the transported quantity in a dimensional form (the momentum, the scalars mass and energy conservation equations), $\vec{v}$ is the speed vector (ms-1), Γ is the diffusion coefficient (m2s-1), and S represents the source term [26].

The turbulent nature of the air flow was simulated using the standard turbulence model k-ε, a model widely used and validated in studies focused on greenhouses, which has shown an adequate fit and accuracy with a low computational cost [27, 28]. Because wind speeds are lower in some areas inside the screenhouse, the effects of buoyancy influenced by the change in air density will be present [29, 30]. Therefore, they were modeled using the Boussinesq approximation, which is calculated using Eq. (2) and Eq. (3).

$\left(\rho-\rho_{0}\right) g=-\rho_{0} \beta\left(\mathrm{T}-\mathrm{T}_{0}\right) \mathrm{g}$          (2)

$\beta=-\left(\frac{1}{\rho}\right)\left(\frac{\partial \rho}{\partial \mathrm{T}}\right)_{\mathrm{p}}=\frac{1}{\rho} \frac{\mathrm{p}}{\mathrm{RT}^{2}}=\frac{1}{\mathrm{T}}$          (3)

where, g is the gravitational constant in (m s-2); β is the volumetric thermal expansion coefficient (°K-1); $\rho_{0}$ is the reference density in (kg m-3): R is the gas constant (J K-1 mol-1); p is the pressure on Pa, and $\mathrm{T}_{0}$ is the reference temperature (℃).

Likewise, the energy equation and the selected radiation model were considered, i.e. the one of discrete ordinates (DO) with angular discretization. The DO model has been widely used in greenhouse studies [31-34] and screenhouses [1]. This model allows calculating, by means of Eq. (4), the radiation and convective exchanges between the roof, the ceiling and the walls of a structure which, in the case of greenhouses, are treated as semi-transparent media. It is also possible to carry out the climate analysis in night conditions, simulating and solving the phenomenon of radiation from the floor of the greenhouse to the outside environment. For this purpose, the sky is considered as a black body with an equivalent temperature (TC) for two predominant scenarios of cloudy and wet nights and clear wet nights [35-37].

$\nabla .\left(I_{\lambda}\left( \begin{matrix} \Rightarrow \\ r \\\end{matrix},\begin{matrix} \Rightarrow \\ s \\\end{matrix} \right) \begin{matrix} \Rightarrow \\ s \\\end{matrix} \right)+\left(a_{\lambda}+\sigma_{s}\right) I_{\lambda}\left( \begin{matrix} \Rightarrow \\ r \\\end{matrix},\begin{matrix} \Rightarrow \\ s \\\end{matrix} \right)$

$=a_{\lambda} n^{2} \frac{\sigma T^{4}}{\pi}+\frac{\sigma_{S}}{4 \pi} \int_{0}^{4 \pi} I_{\lambda}{\left( \begin{matrix} \Rightarrow \\ r \\\end{matrix},\begin{matrix} \Rightarrow \\ s \\\end{matrix}' \right)} \Phi{\left( \begin{matrix} \Rightarrow \\ s \\\end{matrix}.\begin{matrix} \Rightarrow \\ s \\\end{matrix} '\right)} d \Omega^{\prime}$          (4)

where, $I_{\lambda}$ is the intensity of the radiation at a wavelength; $\begin{matrix}   \Rightarrow   \\   r  \\\end{matrix}$, $\begin{matrix}   \Rightarrow   \\   s  \\\end{matrix}$ they are the vectors that indicate the position and direction, respectively; $\begin{matrix}   \Rightarrow   \\   s  \\\end{matrix}$` is the direction vector of dispersion; $\sigma_{s}, a_{\lambda}$ are the coefficients of dispersion and spectral absorption; n is the refractive index; is the divergence operator; σ is the Stefan-Boltzmann constant (5.669×10−8Wm−2°K−4), Φ,T and Ω are the phase function, the local temperature (°K) and the solid angle, respectively.

The presence of insect screens was modeled using equations derived from the flow of a free and forced fluid through porous materials, taking into account their main characteristics of porosity and permeability [38, 39]. These equations can be derived using Eq. (5), which represents the Forchheimer equation.

$\frac{\partial p}{\partial x}=\frac{\mu}{K} u+\rho \frac{C f}{\sqrt{K}} u|u|$          (5)

where, u is the air velocity (ms-1); μ is the dynamic viscosity of the fluid (kgm-1s-1), K is the permeability of the medium (m2); Cf is the net inertial factor; ρ is the air density (kgm-3), and ∂x is the thickness of the porous material (m). The inertia factor Cf and the permeability of the screen K have been evaluated in different experimental studies from tests in the wind tunnel, the numerical results obtained in the experiments are adjusted to equations showing correlation with the porosity (ε) of the screen. The aerodynamic parameters for the insect-proof porous screens commonly used in protected agriculture are obtained by Eq. (6) and Eq. (7), which are the mathematical expressions that best fit the data obtained in the wind tunnel [39-42].

$C_{f}=0.00342 \varepsilon^{-2.5917}$          (6)

$K=2 \mathrm{X} 10^{-7} \varepsilon^{3.5331}$          (7)

Figure 3. Meshing of the computer domain

2.4 Computational domain and generation of the mesh

The construction of the computational domain, the meshing and the evaluation of the quality of the mesh were carried out following the existing guidelines of the good practices of CFD simulation, where the minimum criteria to be met for these three parameters that are directly related to the precision of the results and the required computational effort are established. The ANSYS ICEM CFD 18.2 preprocessing software was used to generate a large computational domain composed of the screenhouse (Figure 3b) and its surroundings, in order to guarantee an appropriate definition of the atmospheric boundary layer and avoid the generation of forced flows with velocities and unrealistic behaviors [43]. The dimensions of the computational domain were 184, 75 and 194 m for the X, Y and Z axes, respectively (Figure 3a). This size was determined following the recommendations of numerical studies of the wind environment around the buildings [44]. The computational domain was divided into an unstructured mesh of hexahedral elements composed of a total of 7,787,701 discretized volumes in space. This number of elements was obtained after verifying the independence of the numerical solutions from the airflow and the temperature behavior at a total of 7 different sized meshes where the one with the highest number of elements presented a value of 12,123,456 and the one with the lowest number of elements was 1,345,123. The independence test was performed following the procedure reported and used successfully by Villagran et al. [34]. The quality parameters evaluated in the mesh were the variation of cell-to-cell size, which showed that 92.3% of the cells in the mesh were within the high-quality range (0.9-1), and on the other hand, the criterion of orthogonality was evaluated, where the minimum value obtained was 0.92, results that are classified within the adequate quality range [45, 46].

2.5 Boundary conditions and convergence criteria

The CFD ANSYS FLUENT 18.2 processing software was used to perform the simulations under the conditions set forth in Table 1. It was run from a computer that was composed of an Intel® Xeon W-2155 processor with twenty cores at 3.30 GHz and 128 GB of RAM in a Windows 10 64-bit operating system. The semi-implicit solution method for the pressure-velocity equation (SIMPLE) was applied to solve the flow field of the simulated fluid. The convergence criteria of the model were established in 10-6 for all the equations considered [47]. With this computer equipment and the simulation criteria established to have an efficiency of results as well as agility and communication effort, the average simulation time was 53 hours for approximately 7,900 iterations.

The upper limit of the domain and the surfaces parallel to the flow were set with boundary conditions of symmetrical properties so as not to generate frictional losses of the air flow in contact with these surfaces. The simulations considered the atmospheric characteristics of the air and the physical and optical properties of the materials within the computational domain, are summarized in Table 1. At the lower limit and the walls of the greenhouse a non-slip wall boundary condition was fixed, at the left boundary, the entry condition for the average wind speed was imposed through a logarithmic profile [48]. The profile was linked to the main CFD module using the user-defined function and using the Eq. (8).

$v(y)=\frac{v^{*}}{K} \ln \left(\frac{y+y_{o}}{y_{o}}\right)$          (8)

where, $y_{o}$ is the roughness of the surface, that for this case, was set at 0.03 m according to the response standard of the local terrain, v* is the friction velocity v(y) is the average wind speed at height and above the ground level and K is the von Karman constant with a value of 0.42, the leeward limit was considered as an edge condition of pressure output type. The model is not considered without the presence of cultivation, since we want to obtain a solution independent of any type and size of plant. In addition, other boundary conditions imposed in the computational domain and the physical and optical properties of materials taken from works such as those of Flores Velasquez et al. [42] and Villagran et al. [49] are summarized in Table 1.

Table 1. Settings of the computational fluid dynamics (CFD) model simulations and boundary conditions

Boundary conditions

Entry domain

Velocity inlet logarithmic profile (Air velocity a 2 m height), and atmospheric pressure.

Output domain

Pressure outlet (Zero pressure and same condition of turbulence).

Treatment of porous medium

Screen porous jump, viscosity effect (α)=3.98 e-9 and drag coefficient (C2) 19185.

Heat source

Constant from the ground, Boussinesq hypothesis activated in the buoyancy effect of the turbulence model.

Physical and optical properties of the materials used






Density (ρ, kg m-3)





Thermal conductivity (k, W m-1 K-1)





Specific heat (Cp, J K-1 kg-1)





Coefficient of thermal expansion (K-1)








Scattering coefficient





Refractive index











Table 2. Initial boundary conditions for simulated configurations


Diurnal Period


Wind speed [ms-1

Wind direction [°]

Air Temperature [°C]

Solar radiation [W m-2]

Dry configuration (DC)





Rain configuration (RC)






Nocturnal Period


Wind speed [ms-1

Wind direction [°]

Air Temperature [°C]

Tc* [°C]

Dry configuration (DC)





Rain configuration (RC)





* Equivalent temperature of the sky.

Table 3. Initial boundary conditions to validate simulation



Wind speed [ms-1]

Wind direction [°]

Air Temperature [°C]

Solar radiation [W m-2]

Tc* [°C]

Diurnal Period







Nocturnal Period







* Equivalent temperature of the sky.

2.6 Measurements and experimental procedure

During the development of the experimental phase between July 01 and July 10, 2018, and in order to obtain data for the validation of the CFD model, ten-minute records of climatic variables inside and outside the SH were made. Outside, a conventional I-Metos weather station (Pessl Instruments Gmbh, Weiz, Austria) was used, located 50 m from the greenhouse and equipped with temperature sensors (range: -30℃ to 99℃, accuracy: ± 0.1℃), relative humidity (range: 10% to 95%, Accuracy: ± 1%), global solar radiation (range: 0 Wm-2 to 2,000 Wm-2, accuracy: ± 2%), wind speed (range: 0 ms-1 to 70 ms-1, precision: ± 5%), wind direction (range: 0° to 360°, resolution: 2° precision: ± 7°) and precipitation (range: 6.5 cm per measurement period; resolution: 0.01 cm; precision: ± 0.1%). The indoor air temperature of the screenhouse was registered by nine data-loggers type sensors HOBO® Pro RH-Temp H08-032-08 (Onset Computer Corp., Pocasset, USA) This measured the temperature in a range from −20℃ to 70℃ with accuracy of ±0.3℃, sensors that were located at a height of Y = 1.8 m above the ground level just at the center line of the screenhouse at X = 17 m and distributed uniformly on the longitudinal Z-axis = 40 m. additionally these devices were covered with a capsule that acted as a protective shield against direct solar radiation.

2.7 Simulated scenarios

The validated CFD numerical model was used as a simulation tool to determine the thermal and aerodynamic behavior of the screenhouse, evaluating two specific configurations, one with rain (RC) and the other one dry (DC), under diurnal and nocturnal climate conditions, establishing the initial conditions listed in Table 2.

2.8 Validation of the model developed

The validation of the CFD model was performed by comparing temperature data obtained experimentally in the SH and the data obtained by numerical simulation for two specific conditions, the initial boundary conditions were determined from the average values of the climatic variables obtained for the experimental period considering a specific time for day and night, respectively (Table 3). Validation is a necessary phase in order to adequately verify the results obtained from the numerical model and establish total independence of these parameter results such as the quality and size of the mesh [50].

Another way to evaluate the performance and accuracy of numerical models is through the calculation of some goodness-of-fit criteria that compare measured and simulated data. In this case, the mean absolute error (MAE) was calculated with Eq. (8), the mean square error (MSE) with Eq. (9) and finally the mean percentage error (MAPE) with Eq. (10).

$M A E=\frac{1}{n} \sum_{i=1}^{n}|X m i-X s i|$          (9)

$M S E=\frac{1}{n} \sum_{i=1}^{n}|X m i-X s i|^{2}$          (10)

$M A P E=\frac{\sum_{i=1}^{n} \frac{|X m i-X s i|}{|X m i|}}{n}$          (11)

where, Xmi is the value measured, Xsi is the simulated value and n the number of data compared. Once it is verified that the values of the goodness-of-fit criteria are close to 0, the model is validated and can be used to develop CFD simulations under the scenarios considered in this investigation.

3. Results and Discussion

3.1 Validation of the CFD model

The fit and performance of the CFD model were tested through a quantitative analysis. For the diurnal period it was found that the absolute differences between the values of the simulated and measured points ranged between 0.25°C and 1.08°C, meanwhile, for the night period, such differences were 0.17°C and 0.83°C. The Figure 4 shows the trend of the simulated and measured data under climatic conditions for the day and night period, it can be seen that the qualitative and quantitative behavior of the data sets are similar, which allows us to deduce at first that the CFD model makes adequate temperature predictions for SH studied.

On the other hand, for the goodness-of-fit criteria used to evaluate the numerical model, values of 0.70℃ and 0.55℃ were obtained for the MAE and MSE, respectively, and an MAPE value of 1.46% for the daytime condition, while for the night-time simulation conditions, values of 0.54℃ and 0.32℃ were obtained for the MAE and MSE, respectively, and an MAPE value of 1.32%. These values obtained for the temperature are in the same order of magnitude as those found by Ali et al. [51].

These experimental results allow us to conclude that the CFD numerical model has adequate capacity to predict the temperature behavior within the SH. Although no experimental measurement of airflow patterns was performed, it is known that the thermal behavior indoors is dependent on airflow patterns, therefore this model can be used as a tool to perform aerodynamic and thermal analysis within the SH structure.

3.2 Daytime period

3.2.1 Air flows

In Figure 5, the behavior of the wind speed inside the SH for the evaluated scenarios RC and DC was observed. For DC an air flow is observed with an average speed of 0.21 ms-1 and maximum and minimum values of 0.49 ms-1 and 0.06 ms-1, respectively (Figure 5a). The behavior of the flow shows a pattern with a higher air velocity in the roof area of the SH over the central length of the structure and that is directed towards the leeward wall. This behavior has already been described in previous studies [11, 42, 52]. For the DC scenario, two converging flow zones can be observed between the ground and the deck area, the zone located between the windward wall and X = 12 m, with a flow of low speeds and in the opposite direction to the outside air flow. On the contrary, the zone between X = 12 m and the leeward window shows higher velocities and a flow that has the same direction of the external flow for the upper average height of the SH, and an air movement for the height at the lower half in the inverse direction to the external flow. Likewise, the interaction area between the windward wall and the roof area of the SH has vectors of low intensity and speed, this is caused by the loss of impulse generated on the air flow by the presence of the insect-proof porous screen (Figure 5a). For the RC case, it was observed that the displacement of the air moves in a single convective cell, clearly differentiated in comparison with DC; this shows a movement in the same direction of the external air flow with average velocity values of 0.24 ms-1, just in the area above the small plastic tunnels and a reverse flow direction in the lower area of these with an average speed of 0.36 ms-1 (Figure 5b).

Figure 4. Comparison of simulated and measured temperature data

Figure 5. Simulated air velocity field inside the screenhouse (m s–1). (a) The configuration of DC, and (b) The configuration RC for the diurnal period

Figure 6. Normalized air velocity (Vint/Vext) inside the screenhouse during the diurnal period for the DC and RC configurations

In order to compare the airflow velocities inside the structure for both RC and DC, the normalized wind speed (VN) was calculated for the heights above the ground level Y = 1 m and Y = 2 m, respectively. This velocity represents the relationship between the interior and the exterior velocity of the air. In Figure 6 the VN curves for RC and DC in each of the heights evaluated along the width of the SH can be observed. For RC-1m a reduction of the air flow was found, which oscillated between 56% and 99.6% in comparison with the external air, meanwhile, the zone with lower velocities appears over the 5 m adjacent to the lateral leeward and windward walls. On the contrary, the highest velocities occurred on the area between 8 m and 30 m of the width of the SH, because the reduction of air flow is influenced by the strong pressure drops that are generated when the external air makes contact with the screen [52]. In the case of RC-2m, a higher reduction in airflow is observed, influenced mainly by the presence of the plastic covering the tunnels; in this case, the air flow values are below, between 80% and 98.8% respect to the external airspeed and its behavior is more homogeneous over the length evaluated in the width of the SH.

In the case of DC-1m and DC-2m, lower flow reductions are observed compared to the RC scenario. In this case, the reductions in airflow are between 26% and 88% in relation to the external wind speed; nonetheless, these values coincide with previous studies conducted by Flores-Velazquez et al. [53]. The behavior for both DC-1m and DC-2m is very similar, the greatest reductions in flow are observed between the windward wall and the zone that is 10 m adjacent to the wall; on the contrary, the lower reduction rates can be observed between 10 m and 29 m of the width of the SH, meanwhile, in the area between 29 m and the leeward wall, an increase in the air reduction indexes is observed again (Figure 6). This allows deducing that the presence of plastic tunnels inside the SH, generates greater reductions in air flow and spatial behaviors differentiated from this flow in RC compared with DC.

3.2.2 Thermal behavior

In Figure 7, the spatial behavior of the temperature inside the SH at a height of 2 m above the ground level is shown. For the DC scenario, an average temperature of 37.6 ± 0.2°C was obtained with maximum and minimum values of 38.3°C and 36.9°C, respectively. Qualitatively it was observed that the areas of higher temperature were generated near the windward side wall just where the lowest airflow velocity values are presented; on the contrary, the lower temperature zones were found in the areas near the front and leeward walls and a small area of the windward wall as well as over the area of greater air flow (Figure 7a). The vertical distribution of the temperature for this case showed a behavior directly related to the air movements as was shown by Teitel et al. [2], finding an area with values higher than 38°C, just in the area of interaction of the two convective airflow cells generated in DC, area that extends from the ground to the SH cover (Figure 7b).

In the RC scenario, the mean temperature value found was 35.5 ± 0.4℃ with maximum and minimum values of 36.5℃ and 34.1℃, respectively. The spatial distribution in the interior volume showed three areas adjacent to the windward wall as the zones of higher temperature, and these zones expanded heterogeneously across the width of the SH. On the other hand, the zones of low temperature were located near the leeward lateral wall and expanded on an area of the central part of the SH (Figure 7c). The distribution of the vertical temperature profile shows high-temperature zones just above the plastic tunnels and another zone with similar values, located on the area between the windward wall and the first plastic tunnel, an area that has low airspeeds and little air exchange (Figure 7d).

Figure 7. Simulated temperature profiles (°C) inside the screenhouse. (a) Top view of the DC configuration at 2 m of height, (b) Front view of the DC configuration, (c) Top view of the RC configuration at 2 m of height, and (d) Front view of the RC configuration for the diurnal period

Figure 8. Thermal gradient profile width of the screenhouse for DC and RC configurations during the diurnal period

The thermal gradient was calculated (∆T) for both DC and RC; this ∆T represents the difference between the air temperature outside and inside the SH. Figure 8 shows the ∆T for a height above the ground level of 1m and 2 m. In general, the average value of ∆T for RC is superior in 0.7℃ and 1.1℃ compared to the DC scenario for both 1m and 2m, respectively. Additionally, it was observed that the behavior of ∆T in RC-2m presents greater variability between nearby points obtaining values of ∆T between 1.7℃ and 2.5℃; this can be directly related to the presence of plastic tunnels (Figure 8).

3.3 Night period

3.3.1 Air flows

The distribution patterns of the air flow for RC and DC are presented in Figure 9. In the case of DC, two air movement cells can be observed; one moves in a clockwise direction from the central area of the SH towards the leeward wall with average speeds of 0.11 ms-1 with some zones of greater speed in the area adjacent to the roof and floor of the structure. On the other hand, the movement cell included between the windward wall presents a clockwise displacement in the upper part of the SH with average velocity values of 0.07 ms-1. This cell is complemented by another that shows a counter-clockwise displacement in the lower part of the structure, with average speeds of 0.10 ms-1 (Figure 9a).

This behavior differs from what was observed by Montero et al. [36] in greenhouses covered with impermeable plastic walls and may be influenced by air leaks through the porous material that occurs both on the front and side walls of the SH. For RC, a clearly differentiated air flow behavior was observed in two zones. A flow in the upper part of the screenhouse just above the plastic tunnels shows slight upward-downward currents that move from the lateral wall of the windward side of the leeward side wall, with an average velocity of this flow of approximately 0.13 ms-1. The other flow moves through the lower area of the plastic tunnels with two main characteristics, the first, a displacement contrary to the flow of the outside air, and the second, a more accelerated air velocity with approximate average values of 0.23 ms-1; this higher speed can be influenced by a greater air movement generated by free convection from a warm zone with lower air density (Figure 9b).

In Figure 10, the behavior of the wind speed (VN) can be observed for the two heights evaluated in RC and DC. The average reduction of the indoor air velocity compared to the outside air for DC was 68% and 81% for DC-1m and DC-2m, respectively. For this scenario, the air displacement patterns move in the direction of the flow of the outside air, except for DC-2m in the area between X = 4 m and X = 8 m on the width of the SH. In the RC scenario, an air movement is observed that moves in the opposite direction to the external airflow for the two evaluated heights, the VN values show a reduction of the airspeed in ranges of 31% and 88%, where the flow pattern that showed a more homogeneous velocity behavior was obtained for RC-1m, unlike that obtained in RC-2m, where high and low-velocity vectors are observed in relatively close points; this is clearly influenced by the presence of the plastic tunnels (Figure 10).

Figure 9. Simulated air velocity field inside the screenhouse (m s–1). (a) DC configuration; and (b) RC configuration during the nocturnal period

Figure 10. Normalized air velocity (Vint/Vext) inside the screenhouse for the DC and RC configurations during the nocturnal period

3.3.2 Thermal behavior

The spatial distribution of the temperature inside the SH for the night period can be seen in Figure 11. For the DC scenario an average temperature value of 20.4 ± 0.2°C, the spatial behavior of this variable was homogeneous inside the structure, the minimum and maximum values obtained under this condition were 20.2°C and 21.3°C, respectively. The areas of higher temperature were located near the front and side walls of the structure, and a small cell generated towards the center of the same, meanwhile the low-temperature zones were located in the central area of the SH between the coordinates of  X = 7 m and X = 14 m over the width dimension of the structure (Figure 11a). The vertical distribution of the temperature for DC shows a behavior where the soil in the central part of the structure is the zone of higher temperature with values close to 21.5°C, approximately in 20% of the volume evaluated. Additionally, low-temperature areas can be observed with mean values of 20.2°C in the central zone adjacent to the higher temperature area (Figure 11b).

Figure 11. Simulated temperature profiles (°C) inside the screenhouse. (a) Top view of the DC configuration at 2 m of height, (b) Front view of the DC configuration, (c) Top view of the RC configuration at 2 m of height; and (d) Front view of the RC configuration during the nocturnal period

Figure 12. Thermal gradient profile width of the screenhouse for the DC and RC configurations during the nocturnal period

The behavior for RC is presented in Figure 11 c and d. The average value of the temperature for this case was 23.9 ± 0.4°C; under these conditions a heterogeneity was observed in the spatial distribution of temperature where there are two zones clearly differentiated, one of them located between the central area of the screenhouse and the wall of the rear facade with average temperature values of 24.3°C, meanwhile, the area of low temperatures was located from the middle zone of SH to the front wall with values of 23.4°C (Figure 11c). The vertical distribution of the temperature shows the presence of two higher temperature zones with average values of 24.2°C, the first located between the soil and the cultivation beds, and a second located on the lower part of the plastic tunnels. On the other hand, the zones of lower temperature with average values of 23.3°C were located near the side walls and the roof of the screenhouse (Figure 11d).

Figure 12 shows the ∆T calculated for RC and DC for heights above the ground level of 1m and 2 m. One of the main differences observed is the numerical value of ΔT for each case, are on the one hand, a positive ∆T which occurs in the RC scenario for the two evaluated heights, with an average ∆T value of 0.7°C, and areas with ∆T values of 0.05°C, just over the areas surrounding the windward and leeward side walls, and areas with ΔT of 1.3°C in the central area of the SH. So it can be inferred that the presence of plastic tunnels may be influencing the thermal behavior, whereas for RC-2m, a behavior with greater variability between nearby points was found (Figure 12).

The opposite occurred in the DC scenario where the structure enters thermal inversion conditions; this phenomenon is characterized by lower indoor air temperatures compared to the values of the outside air temperature. Numerically, this could be checked when analyzing the ∆T values generated under this condition, finding an average value of ∆T of -0.5°C; however, these values are within the range of those reported in previous studies by Teitel et al. [21]. The maximum ΔT value was 0.3°C and occurred in the central zone of the SH and the minimum ∆T value was -0.9°C, which was located in DC-1m in an area on the central zone of the SH and moved in the direction of the lateral leeward window (Figure 12). The thermal inversion phenomenon occurs due to the cooling generated by infrared thermal radiation, poor ventilation and the presence of climatic conditions of low humidity and clear skies [36].

4. Conclusions

The research results are as follows:

(1) CFD 3D simulation proved to be an optimal, valid and accurate tool to determine the microclimatic behavior of a screenhouse during the day and night period under the climatic conditions of the study region.

(2) The presence of small tunnels inside the structure (RC) generates a negative effect on the speed and distribution of air patterns which translates into thermal conditions with values of ΔT up to 1.1°C, compared to the scenario where the tunnels are not used (DC).

(3) The presence of the small tunnels (RC) for the night period allows to improve the microclimate of the screenhouse limiting the phenomenon of thermal inversion characteristic of the DC scenario. 

(4) Any modification to the cultivation system under screen-house structures generates both positive and negative effects on the microclimate, therefore, these modifications cannot be made by following the farmers' criteria alone. It is recommended that in future studies, starting from the base that this investigation leaves with a validated CFD model, other variables of interest are included like, different commercial anti-insect meshes that have defined their aerodynamic properties, evaluations to shorter temporary scales that allow to simulate different meteorological conditions for the day and the night, evaluations with some type of crop or another geometric configuration for the screen-house, by the side of the experimentation it is recommendable to validate the flow patterns through sonic anemometry.


The authors wish to thank Corporación Colombiana de Investigación Agropecuaria (AGROSAVIA) and Instituto Nacional de Innovación y Transferencia en Tecnología Agropecuaria de Costa Rica (INTA) for their technical and administrative support in this study. The research was funded by The Regional Fund of Agricultural Research and Technological Development (FONTAGRO) as part of the project "Innovations for horticulture in protected environments in tropical zones: an option for sustainable intensification of family farming in the context of climate change in LAC".





screen houses

Computational fluid dynamics

Rain configuration


Dry configuration


Observed temperature data (°C)


Simulated temperature data (°C)



Gravitational acceleration, (m.s-2)

Thermal conductivity (W.m-1. K-1)


The normalized wind speed


Absolute mean error (°C)


Mean square error (°C)


Mean absolute percentage error (%)


The semi-implicit solution method for the pressure-velocity equation


source term


Air temperature (°C)


Temperature of the sky (°C)


components of speed (ms-1)


user defined function


Air speed (m.s-1)


the length of the roughness coefficient (m)

Greek symbols


the diffusion coefficient 


thermal gradient (°C)


thermal expansion coefficient (K-1)


turbulent kinetic energy dissipation rate



dynamic viscosity (kg.m-1.s-1)


turbulent viscosity (kg.m-1.s-1)


density (Kg.m-3)


Concentration of the transported quantity in a dimensional form


[1] Bartzanas, T., Katsoulas, N., Kittas, C. (2012). Solar radiation distribution in screenhouses: A CFD approach. In VII International Symposium on Light in Horticultural Systems, 956: 449-456.

[2] Teitel, M., Liang, H., Tanny, J., Garcia-Teruel, M., Levi, A., Ibanez, P.F., Alon, H. (2017). Effect of roof height on microclimate and plant characteristics in an insect-proof screenhouse with impermeable sidewalls. Biosystems Engineering, 162: 11-19.

[3] Tanny, J., Cohen, S. (2003). The effect of a small shade net on the properties of wind and selected boundary layer parameters above and within a citrus orchard. Biosystems Engineering, 84(1): 57-67.

[4] Shahak, Y., Gal, E., Offir, Y., Ben-Yakir, D. (2008, October). Photoselective shade netting integrated with greenhouse technologies for improved performance of vegetable and ornamental crops. In International Workshop on Greenhouse Environmental Control and Crop Production in Semi-Arid Regions, 797: 75-80.

[5] Manja, K., Aoun, M. (2019). The use of nets for tree fruit crops and their impact on the production: A review. Scientia Horticulturae, 246: 110-122.

[6] Tanny, J. (2013). Microclimate and evapotranspiration of crops covered by agricultural screens: A review. Biosystems Engineering, 114(1): 26-43.

[7] Möller, M., Cohen, S., Pirkner, M., Israeli, Y., Tanny, J. (2010). Transmission of short-wave radiation by agricultural screens. Biosystems Engineering, 107(4): 317-327.

[8] Ilić, Z.S., Milenković, L., Šunić, L., Fallik, E. (2015). Effect of coloured shade‐nets on plant leaf parameters and tomato fruit quality. Journal of the Science of Food and Agriculture, 95(13): 2660-2667.

[9] Mahmood, A., Hu, Y., Tanny, J., Asante, E.A. (2018). Effects of shading and insect-proof screens on crop microclimate and production: A review of recent advances. Scientia Horticulturae, 241: 241-251.

[10] Teitel, M., Peiper, U.M., Zvieli, Y. (1996). Shading screens for frost protection. Agricultural and Forest Meteorology, 81(3-4): 273-286.

[11] Tanny, J., Pirkner, M., Teitel, M., Cohen, S., Shahak, Y., Shapira, O., Israeli, Y. (2014). The effect of screen texture on air flow and radiation transmittance: laboratory and field experiments. Acta horticulturae, (1015): 45-51.

[12] Tanny, J., Haijun, L., Cohen, S. (2006). Airflow characteristics, energy balance and eddy covariance measurements in a banana screenhouse. Agricultural and Forest Meteorology, 139(1-2): 105-118.

[13] Pirkner, M., Tanny, J., Shapira, O., Teitel, M., Cohen, S., Shahak, Y., Israeli, Y. (2014). The effect of screen type on crop micro-climate, reference evapotranspiration and yield of a screenhouse banana plantation. Scientia Horticulturae, 180: 32-39.

[14] Cohen, S., Raveh, E., Li, Y., Grava, A., Goldschmidt, E.E. (2005). Physiological responses of leaves, tree growth and fruit yield of grapefruit trees under reflective shade screens. Scientia Horticulturae, 107(1): 25-35.

[15] Desmarais, G., Ratti, C., Raghavan, G.S.V. (1999). Heat transfer modelling of screenhouses. Solar Energy, 65(5): 271-284.

[16] Siqueira, M.B., Katul, G.G., Tanny, J. (2012). The effect of the screen on the mass, momentum, and energy exchange rates of a uniform crop situated in an extensive screenhouse. Boundary-layer Meteorology, 142(3): 339-363.

[17] Meneses, J.F., Baptista, F.J., Bailey, B.J. (2007). Comparison of humidity conditions in unheated tomato greenhouses with different natural ventilation management and implications for climate and Botrytis cinerea control. Acta Horticulturae, 801(801): 1013-1020.

[18] Teitel, M., Garcia-Teruel, M., Ibanez, P.F., Tanny, J., Laufer, S., Levi, A., Antler, A. (2015). Airflow characteristics and patterns in screenhouses covered with fine-mesh screens with either roof or roof and side ventilation. Biosystems Engineering, 131: 1-14.

[19] Villagran, E., Bojaca, C.R. (2019). CFD simulation of the increase of the roof ventilation area in a traditional colombian greenhouse: Effect on air flow patterns and thermal behavior. International Journal of Heat and Technology, 37(3): 881-892.

[20] Mesmoudi, K., Meguellati, K., Bournet, P.E. (2017). Thermal analysis of greenhouses installed under semi arid climate. International Journal of Heat and Technology, 35(3): 474-486. 

[21] Teitel, M., Garcia-Teruel, M., Alon, H., Gantz, S., Tanny, J., Esquira, I., Soger M., Levi, A., Schwartz, A., Antler, A. (2014). The effect of screenhouse height on air temperature. Acta Horticulturae, 1037: 517-523.

[22] Norton, T., Sun, D.W., Grant, J., Fallon, R., Dodd, V. (2007). Applications of computational fluid dynamics (CFD) in the modelling and design of ventilation systems in the agricultural industry: A review. Bioresource Technology, 98(12): 2386-2414.

[23] Bournet, P.E., Boulard, T. (2010). Effect of ventilator configuration on the distributed climate of greenhouses: A review of experimental and CFD studies. Computers and Electronics in Agriculture, 74(2): 195-217.

[24] Flores-Velazquez, J., Ojeda, W., Villarreal-Guerrero, F., Rojano, A. (2015). Effect of crops on natural ventilation in a screenhouse evaluated by CFD simulations. In International Symposium on New Technologies and Management for Greenhouses-GreenSys2015, 1170: 95-102.

[25] Peel, M.C., Finlayson, B.L., McMahon, T.A. (2007). Updated world map of the Köppen-Geiger climate classification. Hydrology and Earth System Sciences Discussions, 4(2): 439-473.

[26] Piscia, D., Montero, J.I., Bailey, B., Muñoz, P., Oliva, A. (2013). A new optimisation methodology used to study the effect of cover properties on night-time greenhouse climate. Biosystems Engineering, 116(2): 130-143.

[27] Drori, U., Dubovsky, V., Ziskind, G. (2005). Experimental verification of induced ventilation. Journal of Environmental Engineering, 131(5): 820-826.

[28] Villagrán, E.A., Bojacá, C.R. (2019). Effects of surrounding objects on the thermal performance of passively ventilated greenhouses. Journal of Agricultural Engineering, 50(1): 20-27.

[29] Villagrán, E.A., Bojacá, C.R. (2019). Determination of the thermal behavior of a Colombian hanging greenhouse applying CFD simulation. Revista Ciencias Técnicas Agropecuarias, 28(3). 

[30] Villagrán, E.A., Bojacá, C.R. (2019). Simulacion del microclima en un invernadero usado para la producción de rosas bajo condiciones de clima intertropicaL. Chilean Journal of Agricultural & Animal Sciences, 35(2): 137-150.

[31] Baxevanou, C., Fidaros, D., Bartzanas, T., Kittas, C. (2018). Yearly numerical evaluation of greenhouse cover materials. Computers and Electronics in Agriculture, 149: 54-70.

[32] Nebbali, R., Roy, J.C., Boulard, T. (2012). Dynamic simulation of the distributed radiative and convective climate within a cropped greenhouse. Renewable Energy, 43: 111-129.

[33] Yu, Y., Xu, X., Hao, W. (2018). Study on the wall optimization of solar greenhouse based on temperature field experiment and CFD simulation. International Journal of Heat and Technology, 36: 847-854.

[34] Villagrán, E.A., Romero, E.J.B., Bojacá, C.R. (2019). Transient CFD analysis of the natural ventilation of three types of greenhouses used for agricultural production in a tropical mountain climate. Biosystems Engineering, 188: 288-304.

[35] Iglesias, N., Montero, J.I., Muñoz, P., Antón, A. (2009). Estudio del clima nocturno y el empleo de doble cubierta de techo como alternativa pasiva para aumentar la temperatura nocturna de los invernaderos utilizando un modelo basado en la Mecánica de Fluidos Computacional (CFD). Horticultura Argentina, 28: 18-23. 

[36] Camacho, J.I.M., Muñoz, P., Guerrero, M.S., Cortés, E. M., Piscia, D. (2013). Shading screens for the improvement of the night time climate of unheated greenhouses. Spanish Journal of Agricultural Research, 1: 32-46.

[37] Villagrán, E.A., Bojacá, C.R. (2019). Numerical evaluation of passive strategies for nocturnal climate optimization in a greenhouse designed for rose production (Rosa spp.). Ornamental Horticulture, 25(4): 351-364.

[38] Campen, J.B. (2004). Greenhouse design applying CFD for Indonesian conditions. In International Conference on Sustainable Greenhouse Systems-Greensys2004, 691: 419-424.

[39] Valera, D.L., Álvarez, A.J., Molina, F.D. (2006). Aerodynamic analysis of several insect-proof screens used in greenhouses. Spanish Journal of Agricultural Research, 4(4): 273-279.

[40] Miguel, A.F., Van de Braak, N.J., Bot, G.P.A. (1997). Analysis of the airflow characteristics of greenhouse screening materials. Journal of Agricultural Engineering Research, 67(2): 105-112.

[41] Teitel, M. (2007). The effect of screened openings on greenhouse microclimate. Agricultural and Forest Meteorology, 143(3-4): 159-175.

[42] Flores-Velazquez, J., Montero, J.I. (2008). Computational fluid dynamics (CFD) study of large scale screenhouses. In International Workshop on Greenhouse Environmental Control and Crop Production in Semi-Arid Regions, 797: 117-122.

[43] Bournet, P.E., Khaoua, S.O., Boulard, T. (2007). Numerical prediction of the effect of vent arrangements on the ventilation and energy transfer in a multi-span glasshouse using a bi-band radiation model. Biosystems Engineering, 98(2): 224-234.

[44] Tominaga, Y., Mochida, A., Yoshie, R., Kataoka, H., Nozu, T., Yoshikawa, M., Shirasawa, T. (2008). AIJ guidelines for practical applications of CFD to pedestrian wind environment around buildings. Journal of Wind Engineering and Industrial Aerodynamics, 96(10-11): 1749-1761.

[45] ANSYS Fluent, V.18.0. Ansys Fluent Tutorial Guide., accessed on Nov. 21, 2019.

[46] Zhang, X., Wang, H., Zou, Z., Wang, S. (2016). CFD and weighted entropy based simulation and optimisation of Chinese Solar Greenhouse temperature distribution. Biosystems Engineering, 142: 12-26.

[47] Baxevanou, C., Fidaros, D., Bartzanas, T., Kittas, C. (2018). Yearly numerical evaluation of greenhouse cover materials. Computers and Electronics in Agriculture, 149: 54-70.

[48] Richards, P.J., Hoxey, R.P. (1993). Appropriate boundary conditions for computational wind engineering models using the k-ϵ turbulence model. Journal of Wind Engineering and Industrial Aerodynamics, 46: 145-153.

[49] Villagrán, E.A., Bojacá, C. (2019). Study of natural ventilation in a Gothic multi-tunnel greenhouse designed to produce rose (Rosa spp.) in the high-Andean tropic. Ornamental Horticulture, 25(2): 133-143.

[50] Ramponi, R., Blocken, B. (2012). CFD simulation of cross-ventilation for a generic isolated building: impact of computational parameters. Building and Environment, 53: 34-48.

[51] Ali, H.B., Bournet, P.E., Cannavo, P., Chantoiseau, E. (2018). Development of a CFD crop submodel for simulating microclimate and transpiration of ornamental plants grown in a greenhouse under water restriction. Computers and Electronics in Agriculture, 149: 26-40.

[52] Tanny, J., Teitel, M., Barak, M., Esquira, Y., Amir, R. (2007). The effect of height on screenhouse microclimate. In International Symposium on High Technology for Greenhouse System Management: Greensys2007, 801: 107-114.

[53] Flores-Velazquez, J., Villarreal Guerrero, F., Lopez, I.L., Montero, J.I., Piscia, D. (2012, July). 3-Dimensional thermal analysis of a screenhouse with plane and multispan roof by using computational fluid dynamics (CFD). In Ist International Symposium on CFD Applications in Agriculture, 1008: 151-158.