Study of the Effects of Vent Configuration on Mono-Span Greenhouse Ventilation Using Computational Fluid Dynamics

Abstract

The rise in the human population, its density and scarcity of resources require cost-effective solutions for sustainable energy and water resources. Smart and sustainable agriculture is one important factor for future green cities to tackle climate change as a cost-effective solution to save energy and water. However, greenhouses (GH) require consistent ventilation due to their internal temperatures, and this can be an energy-intensive operation. Therefore, it is necessary to analyse the potential factors involved. In this study, the effect of vent configuration of a mono-span greenhouse with roof and side vents at low wind speeds was investigated using computational fluid dynamics (CFD). The validated simulations were then performed on different models to analyse the effects of the vents’ locations on the ventilation requirements. The side vents were found to contribute most to the ventilation. The position of the side vent was found to affect the convection loop in the greenhouse and the air velocity at the plant level. The humidity was shown to be highest under the windward side vent. The roof vent was found to affect the temperature and air velocity in the roof of the greenhouse but had very little effect on the distributions at the plant level.

1. Introduction

A greenhouse (GH) is an indoor production facility with a transparent cover, used to grow plants under ideal conditions, as it allows good control of the external conditions in which the plants are grown, such as temperature, air velocity and flow pattern, ventilation rate, humidity and CO₂ concentration. Traditionally, GHs are used in cold climates to protect crops from the wind and increase the internal temperature, to grow crops that would usually not be able to grow in cold climates. However, recently, smart GHs are also being used in warm regions to provide shading and remove excess heat and humidity, so that plants can be grown more effectively. This is especially important in the Middle East and North African (MENA) region, with scant water resources [1,2]. Furthermore, high temperatures with intense solar energy and low humidity cause open field agriculture to be less economical [2], while available intense solar energy provided situations for implementing sustainable technologies [3,4,5,6,7]. By using GHs, the water demands for irrigation can be reduced by up to 90% for countries with water stresses [2]. Linked to this, many countries in the MENA region have deployed desalination technologies for irrigation [8]. Since the temperature inside the GH reaches to a certain level, ventilation, shading and also the cooling technologies are important to provide suitable temperature for the plants’ growth [9]. While water demands can be controlled in the GHs, the ventilation needs to be analysed experimentally and/or computationally to enable the design of cost-effective smart ventilation technologies. Computational Fluid Dynamics (CFD) can simulate the GH micro-climate status and investigate the arrangements of ventilators based on the air velocities. Generally, after development, a CFD model is verified against experimental results, to ensure that CFD results are accurate and simulate physical GHs in the same conditions. Experimental approaches have also been deployed in previous studies to analyse the ventilation requirements on real GHs [10,11] using a tracer gas technique [12]. Some others analysed the distribution of temperature and velocity contours in the GHs by deploying thermocouples placed at certain predefined points to measure the air velocity and temperature [13]. Bartzanas et al. [14,15] utilised a 3D sonic anemometer to observe the air velocities and temperature in a GH for validation of their CFD model. Conjointly, the CFD methods were broadly used to analyse the ventilation requirements for GHs [16,17,18,19,20,21,22]. Although there are several other types of GHs, such as Parral [23], Saw-tooth [24], Canarian [25], Multi-span [26,27], Venlo [20], Open-roof [28,29], and Tropical [30,31] GHs, mono-span types of GHs [20,32,33,34,35] are largely being used in agriculture due to its simple geometry and popularity [36]. This study aims to investigate a mono-span GH with different ventilation requirements under the Egyptian environmental conditions, which are high temperatures, low wind velocities and low air humidity, with data obtained from the literature in Borg Elarab [37]. The effects of the design of the GH, as well as the environmental conditions, on the air-flow pattern, velocity, and the internal temperature and humidity of the GH are analysed using CFD. This helps decision-makers design smart ventilation systems to provide optimal micro-climate conditions for plant growth, which increases the crop yield [38]. These micro-climate conditions generally include the GH temperature, wind speed and relative humidity [24]. Temperature and wind speed were analysed using CFD models [39]; hitherto the relative humidity was simulated only by Kim et al. [40]. Their study analysed the humidity distributions using a fog cooling system and several dehumidifiers which limit the final results, while in our case study, a general pattern of the relative humidity is presented without using any cooling/heating systems. This could help to analyse the micro-climate conditions prior to be used as a sustainable bio-system with natural ventilation [41], without the use of energy resources for forced ventilation or adding cooling systems [42]. Since the studies showed that the ambient conditions, such as wind direction can influence the ventilation [19,43], providing an optimal solution for the ventilation is not possible, but analysing influential parameters can affect the design and implementing of vents and cooling mechanisms for having a sustainable horticultural system.

2. Materials and Methods

The basic theory of air-flow through a vent is described below, with the full analysis in Boulard and Baille [44]. Air-flow through an opening is caused by a pressure difference, which can be described using Bernoulli’s Equation for an incompressible gas with respect to the vertical height [44] (Equation (1)), and the discharge coefficient of the vent (Equation (2)).

$$\Delta P\left(y\right)=\frac{1}{2}\rho \zeta v^2\left(y\right)$$

(1)

$$A_1=\frac{1}{{\zeta }^{0.5}}$$

(2)

The distribution of the air speed can then be deduced from Equations (1) and (2) [44].

$$v\left(y\right)=\frac{\left|\Delta P\left(y\right)\right|}{\Delta P\left(y\right)}{A}_1{\left(\frac{2}{\rho }\left|\Delta P\left(y\right)\right|\right)}^{0.5}$$

(3)

The air-flow rate can then be calculated by integrating the velocity in Equation (3) over the height of the opening [44].

$$G=L{{\displaystyle \int }}_0^{H}v\left(y\right)dy$$

(4)

The buoyancy effect must also be described by the Equations. The vertical distribution of the pressure difference can be deduced if the vertical temperature distribution is assumed to be homogeneous, and the Boussinesq approximation is used [45]. In this case, the pressure distribution is given by Equation (5), and if the air velocity is assumed to be constant over the surface of the vent, the volumetric flow rate is given by Equation (6).

$$\Delta P\left(y\right)=\Delta P_0-\rho g\left(\frac{\Delta T}{T}\right)y$$

(5)

$$G=\frac{S}{2}{A}_1{\left[2g\frac{\Delta T}{T}\frac{H}{4}\right]}^{0.5}$$

(6)

Boulard and Baille [44] also studied the Equations when the buoyancy and wind effects were combined, for the models that they analysed.

The CFD simulations were performed using ANSYS Fluent 19.0, ANSYS Inc., USA [46], which is a commercial code commonly used in different industries [47,48,49,50,51,52,53], and is also being used to analyse the ventilation requirements for optimal agricultural production [54,55,56], mainly to reduce the energy consumption [55,57] based on the temperature/air humidity analysis [58]. ANSYS Fluent uses the Finite Volume method, in which the domain of interest is divided into numerous small, non-overlapping control volumes or cells to develop the mesh. The continuum Equations of conservation of mass, momentum and energy describing the flow are discretised on this mesh by integrating over the volume of each individual cell and applying Gauss’ theorem to convert spatial derivatives into fluxes across the cell faces. The effect of thermal buoyancy is included through the Boussinesq approximation in which the buoyancy is modelled by a temperature-dependent body force in the momentum Equation [59]. Since the airflow is turbulent, the mean flow variables must be solved together with a turbulence model which describes the effect of the turbulence on the mean flow. The standard k-ε turbulence model was used in the simulations, as this is the most common and reliable model used in the literature [46], and most studies have found good convergence and accurate results using the model [20,21,36,60]. The standard k-ε model is a semi-empirical model based on Reynold’s Equation and the Navier–Stokes Equations [61], which solves two Equations for the turbulent kinetic energy (k) and the dissipation rate of this energy (ε) [46]. Equations (7) and (8) [46] are the developed transport Equations used in Fluent to solve for k and ε respectively.

$$\frac{\partial }{\partial t}\left(\rho k\right)+\frac{\partial }{\partial x_i}\left(\rho k{u}_i\right)=\frac{\partial }{\partial x_j}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_k}\right)\frac{\partial k}{\partial x_j}\right]+G_k+G_b-\rho ϵ-Y_M+S_k$$

(7)

$$\frac{\partial }{\partial t}\left(\rho ϵ\right)+\frac{\partial }{\partial x_i}\left(\rho ϵu_i\right)=\frac{\partial }{\partial x_j}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_{ϵ}}\right)\frac{\partial ϵ}{\partial x_j}\right]+C_{1ϵ}\frac{ϵ}{k}\left(G_k+C_{3ϵ}{G}_b\right)-C_{2ϵ}\rho \frac{{ϵ}^2}{k}+S_{ϵ}$$

(8)

The eddy viscosity is computed as displayed in Equation (9) by combining k and ε.

$${\mu }_t=\rho C_{\mu }\frac{k^2}{ϵ}$$

(9)

The constants, along with $C_{\mu }$ in Equation (9), are given default values in the software, which are as follows:

$$C_{1ϵ}=1.44,C_{2ϵ}=1.92,C_{\mu }=0.09,{\sigma }_k=1.0,{\sigma }_{ϵ}=1.3$$

When the force due to gravity is applied to the model, the Boussinesq approximation simplifies the model by ignoring density differences unless the term in the Equation is multiplied by the force due to gravity, which gives rise to buoyancy forces due to temperature differences [46]. Due to the use of the Boussinesq approximation, the Equations for $G_k$ and $G_b$ differ from the standard Equations used in the standard k-ε model as shown in Equations (10) and (12), whilst Equation (11) defines the modulus of the mean rate-of-strain tensor, and Equation (13) defines the coefficient of thermal expansion, which is a property of the fluid [46].

$$G_k={\mu }_t{S}^2$$

(10)

$$S=\sqrt{2S_{ij}{S}_{ij}}$$

(11)

$$G_b=\beta g_i\frac{{\mu }_t}{P{r}_t}\frac{\delta T}{\delta x_i}$$

(12)

$$\beta =\frac{-1}{\rho }{\left(\frac{\delta \rho }{\delta T}\right)}_{\rho }$$

(13)

A full derivation of the Equations, along with the assumptions made in the standard k-ε model, can be found in Mohammadi and Pironneau [61]. The heat transfer Equation is used in the model to solve the transfer of heat between different regions of space. Heat is transferred via conduction, convection and radiation. The energy Equation can be found in Equation (14):

$$\frac{\partial }{\partial t}\left(\rho E\right)+\nabla ·\left(\overrightarrow{v}\left(\rho E+p\right)\right)=\nabla ·\left(k_{eff}\nabla T-{\displaystyle \sum }_j{h}_j\overrightarrow{J_i}+\left({\stackrel{=}{\tau }}_{eff}·\overrightarrow{v}\right)\right)+S_h$$

(14)

For an incompressible flow as used in the model, the Equations used for the energy E and the sensible enthalpy h are defined in Equations (15) and (16).

$$E=h-\frac{p}{\rho }+\frac{v^2}{2}$$

(15)

$$h={\displaystyle \sum _j}{Y}_j{h}_j$$

(16)

In Equation (16), $Y_j$ is the mass fraction of the species j, and $h_j$ is the enthalpy of the species, defined by Equation (17), where T_ref is 298.15 K and $c_{p,j}$ is the heat capacity of the species at constant pressure.

$$h_j={\displaystyle \underset{T_{ref}}{\overset{T}{\int }}}{c}_{p,j}dT$$

(17)

The species transport Equation can be used to study the humidity contours in the GH, as the user can specify mass fractions of different species at different points in the model. The species transport Equation solves conservation Equations describing convection, diffusion and reactions for each species in the model [46]. The general form of the conservation Equation used by Fluent to predict the mass fraction of each species at a given point in the model is given in Equation (18), and for a turbulent flow, the mass diffusion $\overrightarrow{J_i}$ is computed using Equation (19).

$$\frac{\partial }{\partial t}\left(p{Y}_i\right)+\nabla ·\left(\rho \overrightarrow{v}{Y}_i\right)=-\nabla ·\overrightarrow{J_i}+R_i+S_i$$

(18)

$$\overrightarrow{J_i}=-\left(\rho D_{i,m}+\frac{{\mu }_t}{S{c}_t}\right)\nabla Y_i-D_{T,i}\frac{\nabla T}{T}$$

(19)

In this study, it is assumed that no chemical reactions take place in the GH, so the species transport Equation relies solely on the user-defined species production. From the literature review, it was decided that the ventilation and temperatures of a mono-span GH would be analysed using CFD, as this was the most common GH used in agriculture [36]. The air surrounding the GH is also modelled, to reduce the effects of thes pressure gradients found in the model [16,39]. The full GH dimensions can be found in

Table 1. The dimensions of the greenhouse (GH) and external air used for the Computational Fluids Dynamics (CFD) geometry.

Dimension	Value
Length	8.25 m
Width	6.40 m
Ridge Height	3.75 m
Eaves Height	2.15 m
Side Ventilator Vertical Length	0.50 m
Side Ventilator Height (above ground)	a
Roof Ventilator Vertical Length	0.50 m
Roof Ventilator Distance from Ridge	b
External Upstream Air Length	2.00 m
External Downstream Air Length	8.00 m
External Air Height	5.00 m

, with a cross-section shown in Figure 1. It should be noted that for the initial 3D model, both the side and roof ventilators are continuous. Additionally, the external wind inlet is on the left of the geometry, so the ‘windward’ side of the GH is always the left-hand side of the GH shown on figures, and the ‘leeward’ side is always on the right. It should also be noted that the walls of the GH are cut from the surface to simplify the model, so no heat or mass transfer is solved for these areas. Finally, the ridge is the central point of the GH, and the eaves are the points at which the roof meets the wall. The GH model was imported into ANSYS Fluent, primarily as a 2D model. Initial analysis was performed on the 2D GH, as the literature [20,35,60] showed that for a wind perpendicular to the ridge of the GH, 2D and 3D studies gave similar results for the cross-section perpendicular to the ridge. The simulations were performed using a pressure-based solver, and steady-state analysis was performed so that the results produced were time-independent for constant wind speed, in line with most literature. Gravity was enabled, with the gravitational acceleration set to −9.81 ms⁻² in the y-direction (vertical direction).

In this study, the radiation from the sun was modelled using temperature boundary conditions. To replicate the effect of the plant transpiration, the floor of the GH was set to have an H₂O mass fraction output, which is the equivalent of relative humidity (RH), when the humidity model was included in the simulation. The fluid properties for both air and water were left as the default settings. However, as the Boussinesq approximation was used in the model, the Boussinesq density was selected, using the same value as the constant fluid density, and the thermal expansion coefficient for each fluid was also found. The thermal expansion coefficient of air was found to be 0.0034 K⁻¹ at 25 °C [62], and the thermal expansion coefficient of water was 0.000257 K⁻¹ at 25 °C [63]. The solids used in the simulations were glass and soil. The glass cover creates suitable microclimate conditions for the cultivation [64]. The material properties that affect the fluid flow and temperature Equations are the density, specific heat and thermal conductivity of the materials (see

Table 2. The material properties used in the simulation for glass and soil [46].

Material	Density (kg/m3)	Specific Heat (J/kgK)	Thermal Conductivity (W/mK)
Glass	2400	753	1.0
Soil	2200	871	0.5

). The glass was applied to all boundaries of the GH walls and roof, and the soil was applied to the GH floor. However, due to the setup of the model, these are only surface properties, as the wall thickness was not simulated in order to reduce computational requirements, so the effect of the solid properties on the model is small. The boundary conditions are based on experimental measurements by Molina-Aiz et al. [65] for the 2D (see

Table 3. The boundary conditions used in the 2D model [65].

Named Selection	Boundary Type	Boundary Conditions
GH Walls	Wall (glass)	T = 310 K
GH Roof	Wall (glass)	T = 310 K
GH Floor	Wall (soil)	T = 320 K
External Floor	Wall (soil)	T = 300 K
Inlet	Velocity-inlet	U = 0.2 m/s T = 300 K
Outlet	Pressure-outlet	N/A
External Roof	Symmetry	N/A

), and 3D models (see

Table 4. The boundary conditions used in the 3D model. Note: ‘end’ refers to the walls perpendicular to the GH ridge axis, whilst ‘side’ refers to the walls parallel to the GH ridge axis [65].

Named Selection	Boundary Type	Boundary Condition(s)
GH Side Walls	Wall (glass)	T = 310 K
GH End Walls	Wall (glass)	T = 310 K
GH Roof	Wall (glass)	T = 310 K
GH Floor	Wall (soil)	T = 320 K
External Floor	Wall (soil)	T = 300 K
External Walls	Wall (glass)	T = 310 K
Inlet	Velocity-inlet	U = 0.2 m/s T = 300 K
Outlet	Pressure-outlet	N/A
External Roof	Symmetry	N/A

).

A mesh convergence study was performed on the geometry to obtain the optimal mesh density for the developed CFD model. A higher cell density is required where strong gradients occur [36] to reduce computational requirements without losing accuracy, and therefore cell size controls were added near to the GH walls in the model. The simulation for each mesh was run until all residuals reached 10⁻³, and then the temperature and air velocity values and averages were taken over an area. It was found that the results converged at 835534 elements, with a maximum global face size of 10 mm, and a maximum size for each edge in the geometry of 5 mm (See

Table 5. The convergence study for the optimal number of elements.

Mesh	Number of Elements	Average Horizontal Velocity (m/s)	Average Horizontal Temperature (K)	Average Vertical Velocity (m/s)	Average Vertical Temperature (k)
1	8337	0.116	304.74	0.200	305.31
2	33,236	0.101	304.87	0.188	305.29
3	45,692	0.217	303.93	0.173	304.47
4	61,420	0.227	304.24	0.187	304.71
5	63,304	0.083	304.39	0.180	304.68
6	78,501	0.220	304.33	0.187	304.79
7	90,815	0.153	304.20	0.161	304.65
8	101,727	0.132	304.23	0.160	304.67
9	107,417	0.082	305.80	0.212	305.44
10	116,669	0.219	304.31	0.178	304.78
11	206,008	0.154	304.22	0.160	304.66
12	214,404	0.210	303.90	0.185	304.53
13	227,126	0.213	304.17	0.184	304.70
14	272,045	0.212	304.92	0.197	305.23
15	372,999	0.168	304.06	0.196	304.53
16	379,439	0.180	304.27	0.198	304.65
17	419,490	0.204	304.65	0.215	305.03
18	486,947	0.180	304.08	0.182	304.55
19	503,449	0.187	304.22	0.201	304.64
20	571,726	0.175	304.16	0.173	304.54
21	586,455	0.161	304.44	0.195	304.73
22	679,925	0.171	304.18	0.159	304.48
23	696,438	0.161	304.35	0.173	304.56
24	819,862	0.158	304.18	0.143	304.40
25	835,534	0.148	304.30	0.147	304.51
26	877,806	0.151	304.66	0.155	304.69

).

The 2D geometry was also used to produce a 3D model of the initial GH, by ‘extruding’ the geometry to create a 3D model with continuous roof and side vents in the same locations as the 2D model. The 3D model would be used to verify the flows in the 2D model, and to analyse the differences in the results of different planes in the GH (See Figure 2). The same setup and materials were created as those of the 2D model.

The initial 2D model was used as the baseline for the analysis, and each model simulation in the analysis performed had one parameter varying from the initial model, either in the geometry or the boundary conditions of the simulation. The positions of the side and roof vents, closing the side vents and changing the positions of the roof vents, and closing the roof vents and changing the positions of the side vents were all studied to analyse the effect on the air velocity, temperature and humidity. The side vents were studied for heights (at the bottom of the vent) between 0.25 m and 1.5 m above the ground in intervals of 0.25 m, and the roof vents were studied for distances (at the top of the vent) between 0.5 m and 3.0 m from the ridge of the GH in 0.5 m intervals. When one of the vents was closed, the interval distance was doubled for the open vent (i.e., 0.5 m interval for side vents, 1 m interval for roof vents). To analyse the boundary conditions, the wind speed, internal soil temperature and water vapour mass fraction were all individually adjusted to study the effect of each parameter on the initial model. The external wind speed was analysed for values of 0.1 m/s, 0.2 m/s, 0.4 m/s and 0.6 m/s to study the effect of the wind speed on the air-flow patterns and temperature contours in the GH. The soil temperature was analysed for 310 K, 315 K, 320 K and 325 K, as these were all found in experimental measurements of the floor temperature, to obtain different results and compare the effects of different temperatures on the model.

3. Validation

The air-flow patterns, velocities and temperature contours produced by the simulation are in very good agreement with the experimental work done by Sase et al. [35]. The air-flow pattern qualitatively replicates the pattern shown in the results of their work. The temperature contours are different, which is likely due to the different ventilator type used in the model. In terms of quantitative validation, the range of maximum temperature for the current study (320 K) is within the previous studies, 343 K [37] (with less than 7% deviation). In this study, the temperature of the region specified for the plants’ growth is around 300–306 K which is an acceptable maximum range for the cultivation (maximum of 308 K) [66]. The model also provides similar results to the CFD simulations performed by [20], which also attempted to model the experiment carried out by Sase et al. [35] using CFD methods. In this study, the vents themselves were modelled as ‘roll-up’ vents [14] instead of the pivoting vents used in Sase et al. [35]. This causes the velocity of the air through the windward side vent to increase as it passes through the vent, since the area where the air enters the GH through decreases, and therefore, to keep the momentum Equation, the velocity of the air passing through the vent must increase. This means that for an external wind velocity of 0.2 m/s, the average wind velocity entering the side vent is 0.43 m/s. The effect of the pivoting vent, as discussed by Bartzanas et al. [14], is that the air flowing towards the vent must first flow around the vent itself, which causes the air to slow down. The different air-flow pattern caused by this vent type can also be seen in Shyklar and Arbel [19] which shows similar trends as of this study.

4. Results and Discussion

The results of the initial 2D model (Figure 3) show that the air enters through the vents on the windward side of the GH. This air-flow rises over a loop of air, with some of the new air joining the loop, and the rest falling to the floor on the leeward side of the GH. This air then moves up the wall on the leeward side of the GH, with some exiting through the side vent, but most of the air follows the roof up to the roof vent, where it exits the GH. The temperature contour shows that where the loop occurs, the temperature is up to 6K higher than the temperature of the air entering the GH. This is most likely caused by the lower renewal rate of the air, which allows the temperature to rise. The water vapour contour in Figure 3 shows that the most humid area of the GH is under the windward side vent, which is also likely due to the reduced air renewal rate in this area.

The temperature contour is like that shown in the 2D initial model, with a higher temperature below the windward side vent than in the rest of the GH, as shown in Figure 4. The results of the 3D model show a similar air-flow pattern in the centre of the GH to that of the 2D model. An air current enters through the windward side vent, moves towards the floor and forms a clockwise rotating loop of air pattern below the windward side vent. The rest of the air moves along the floor, then up the leeward side wall of the GH, before exiting out of the leeward side and roof vents. In the 3D model, a slow counter-clockwise air loop also forms above the floor in the leeward side of the GH. A similar air-flow pattern forms closer to the walls of the GH, although the counter-clockwise loop of air does not form. However, the top view of the GH shows that the air-flow pattern is not completely independent of the side wall (see Figure 5), so its effect cannot be ignored, especially when the analysis is performed close to the wall. Figure 5 shows that the wall slows the air passing over it and causes a boundary layer profile to form, which in turn produces rotating air patterns close to the wall. These loops of air then affect the rest of the air-flow distribution.

As described previously, this is likely to be due to the low air renewal rate in this area. The temperature contour for the plane 0.50 m from the end wall of the GH also shows that the temperature does not vary much over the length of the GH. However, the end wall affects the temperature distribution (see Figure 5), as the wall is heated so that the air close to it is warmer than the external air, and there are symmetrical warmer contours form either side of the centre of the air-current, proving that the temperature is also not independent of the end wall.

The humidity contour is also like that of the 2D model (see Figure 4) for the 3D initial model in the centre and at the edge of the GH. However, the result shows that the humidity also builds up near the windward wall of the GH (see Figure 5).

The 2D and 3D initial models show similar air-flow patterns in the centre of the GH, as seen in Figure 3 and Figure 12, validating the assumption that the end walls of the GH have a negligible effect on the flow patterns in the centre of the GH. Furthermore, the air-flow pattern closer to the wall of the GH also has a similar pattern to that in the 2D model. However, the top view of the 3D model, seen in Figure 5, suggests that the flow is not independent of the end wall forces, and shows some air loops forming on the end walls that affect the rest of the flow. Comparison of the 2D and 3D temperature contours also shows that the air temperature distribution around the GH is very similar for the 2D and 3D models, with the warmer section of both Figure 3 and Figure 12 found under the windward vent where the air loop is created. This is confirmed by the results close to the end wall, which shows that the contour is still similar close to the wall of the GH. The top view (See Figure 5) shows a symmetrical increase in the temperature of up to 4 K just inside the windward vents, which is caused by the slower moving air in that region.

4.1. Results of the Side Vent Analysis

The results from the side vent geometry analysis show that for low wind speeds, the height of the vent influences the air-flow pattern inside the GH. When the vent is towards the floor of the GH, there is a fast-moving air current along the floor of the GH, as seen in Figure 6. This air causes a counter-clockwise rotating air loop to form above the air current, which fills almost all of the GH. As the side vent moves towards the centre of the wall, like in the initial model (see Figure 3), the counter-clockwise loop decreases in size, and a clockwise loop of air forms below the air entering the GH. When the side vents are at the top of the wall, the counter-clockwise loop no longer exists, and the air enters the windward side vent, follows the roof upwards to the ridge, then flows back down the leeward roof, where some air exits the leeward side vent and most of the air goes towards the floor and back to the windward wall, creating a large clockwise circular air-flow pattern (see Figure 7).

The dominance of the buoyancy effect is demonstrated by the air-flow patterns, as when air along the floor in every simulation is heated, it rises upwards instead of flowing along the floor. The buoyancy effect is the driving force behind the circular air-flows, as when the inlet flows are along the floor, the air rises up and forms a counter-clockwise circular flow, whereas when the vents are at the top of the wall, the stagnant air below the vent is heated and rises up into the inlet air-flow, which causes the clockwise flow of air around the GH.

The position of the side vents also affects the temperature distribution inside the GH, as in every model, the air below the vent has a higher temperature than that of the air entering the GH. For this model, using a vent position at the top of the wall causes temperatures in some locations to increase by 2 K to 4 K when compared to a vent positioned at the bottom of the wall. The temperature contours for a high and low side vent position are displayed in Figure 3 and Figure 4 and show that the temperature of the air in the GH for the horizontal distance from the windward vent, at a height of 1m for a low (0.25 m), central (1 m) and high (1.5 m) side vent position, and the temperature at plant level for the high side vent GH is at least 2 K higher than the temperature of the low side vent for every distance, whilst the temperature for the central side vent design is dependent on the distance from the vent, as the air current from the windward vent affects the temperature at this height. It should also be noted that while the low side vent provides better cooling, it has a higher velocity of air-flow at the plant level, as shown in Figure 6, which may cause damage to plants.

The water vapour mass fraction contour follows a similar pattern to that of the temperature contour. For low vent positions, the water vapour does not rise very far above the floor before it is taken by the air current from the windward side vent, and exits the leeward side vent, as shown in Figure 6. For higher side vents, the water vapour is taken by the air into the clockwise circling air pattern, and so the humidity in the air increases before exiting the leeward side vent, as seen in Figure 7. This also suggests that the ventilation rate is less efficient for higher side vent positions.

4.2. Results of the Roof Vent Analysis

The position of the roof vents does not have as great an impact on the variables studied as the position of the side vents. As the distance of the roof vent from the ridge increases, the clockwise circular air-flow pattern formed in the initial model increases in size to fill more of the GH. A larger percentage of the air from the roof vent joins the air current from the side vent, making the current larger and increasing the velocity of the air. Additionally, the fast-moving air tends to be at the same height as the roof vent, which means that for the low roof vents, such as in Figure 8, there is very little air movement in the roof of the GH.

The temperature contours are like that of the initial model, particularly at the plant level where the temperature increases by between 2 K and 6 K in all simulations, compared with the external air temperature. The effect of the roof vents on the temperature, as was with the air velocity, is that the temperature in the roof of the GH increases by up to 10 K as the roof vents are placed further away from the GH ridge, as seen in Figure 8. However, this is a small area, so the temperature increase would have a minimal effect on the conditions of the GH at the plant level.

The position of the roof vent has no effect on the water vapour contour, with every result being virtually identical to the initial model, as seen in Figure 3.

4.3. Results of Side Vent Analysis without Roof Vents

The air-flow pattern formed when the roof vents are removed, as seen in Figure 3, is still very similar to the air-flow in the initial model, which is shown in Figure 3. The air enters through the windward side vent of the GH and creates the same clockwise air loop below the vent. The difference between the original model is in the outflow of the air, as Figure 3 shows that in the initial model, most of the air exits the GH through the leeward roof vent. However, as Figure 9 shows, there is no roof vent in the model, so most of the air instead exits through the leeward side vent. This causes poor air renewal rate in the roof of the GH.

This is also shown in Figure 10, which plots the air velocities for the vertical height above the floor at the centre of the GH, for a roof vent placed close to the ridge (0.5 m), in the centre of the roof (2 m) and far away from the ridge (3 m). It shows that for the vent placed close to the ridge, the air velocity has two peaks, due to the two separate air currents developed in the initial model (Figure 3), and as the vent moves further from the ridge of the GH, the maximum velocity of the air occurs closer to the floor of the GH. It also shows that there is very little air movement above the air current caused by the roof vents. This can also be seen in Figure 7 where the temperature is analysed against the vertical height above the floor at the centre of the GH, for the initial model and the models without side vents. The position of the side vents is also analysed when the roof vents are removed, and for the high side vent design, the air-flow pattern is very similar.

This shows that the roof vent has very little effect on the air-flow for the high side vents as the clockwise circular air pattern still forms. However, for the low side vent position, the air enters the windward side vent, travels along the floor and out the leeward side vent, so the counter-clockwise air loop does not form, and there is very little air movement in the roof of the GH, as seen in Figure 11.

When the roof vents are further from the ridge (Figure 8) the temperature increases in the roof of the GH when the roof vents are removed, as the warm air is not efficiently ventilated. This means that the temperature in the roof increases by up to 6 K for the low side vent model, seen in Figure 11, when compared to the same GH design with roof vents included in the standard position, as in the initial model. However, as the high side vent design maintains the same clockwise air pattern without roof vents, the average temperature in this design is much lower, maintaining the same temperatures. The initial model shows a much lower temperature than the same model without the roof vent, particularly in the roof, as the temperature is 7 K higher for the closed roof vent model with the side vent in the same position as the original model. The high side vent position maintains the same air-flow pattern without the roof vent, which causes much more effective cooling than the other models analysed since the temperature stays at about 304 K until it reaches the very top of the roof. As with other results, the water vapour stays below the windward vent, where it builds up in the slower-moving air and where the temperature is also higher. Most of the water vapour is found below the windward vent, as the stronger air current flows over the corner and develops the rotating flow patterns. The water vapour is naturally removed by the air current and exits through the leeward side vent.

4.4. Results of the Roof Vent Analysis without Side Vents

The effect of the side vents on the ventilation of the GH is clearly seen by comparing Figure 12 with Figure 3. The air-flow pattern is different, and the average velocity of the air in the GH is lower, which causes a temperature increase of around 9K, and also causes the water vapour mass fraction to greatly increase. This temperature difference (see Figure 13) suggests that the ventilation rate is much lower, which is expected, as the side vents account for most of the air-flow into the GH in the initial model.

The effect of the position of the roof vent can also be studied. In the original roof vent position when the GH has no side vents, the air enters the windward roof vent and forms a clockwise circle of air that fills most of the GH, and some of the air forms a smaller, counter-clockwise circular air pattern by the leeward wall, as seen in Figure 12.

As the vent position moves away from the ridge of the GH, the counter-clockwise rotating air dissolves, and the clockwise loop fills the GH. When the vent is just above the wall, as far from the ridge as possible, the direction of the air-flow changes, as shown in Figure 14. Instead of moving along the roof and down the leeward wall, a fast-moving current of air moves directly from the windward vent to the centre of the floor and forms a counter-clockwise pattern by moving up the leeward wall, with some air exiting the leeward vent and the rest moving up the leeward roof to complete the circular flow. A smaller clockwise circular air-flow also forms from the air moving to the left once it reaches the floor.

The temperature contours show that the designs where the roof vents are further from the ridge are much more efficient at cooling the GH and maintaining a temperature closer to that of the outside temperature, which can be seen by comparing Figure 8 and Figure 9. The most likely cause of this is that as the vents are moved further apart, a lower percentage of the air that enters the GH can exit directly through the leeward vent, so the air renewal is more efficient, and the temperature is lower. The cold air-flow seen from the windward vent on Figure 14 is the fast-moving air described previously, which has a significant effect on the temperature of the rest of the GH. Comparing the temperature profile with the air-flow pattern for both Figure 8 and Figure 9 shows the effect of the flow on the temperature of the air, as both figures show that where the cooler external air enters the GH at relatively high velocities, the temperature contours show a lower temperature in that area.

The temperatures in the GH at different heights above the floor are compared in Figure 13, for the initial model; and for the model with a roof vent placed close to the ridge (0.5 m), slightly further from the ridge (1 m), in the centre of the roof (2 m), and far from the ridge (3 m); each without a side vent. The effect of the side vent can be seen by comparing the initial model with the same roof vent position (0.5 m) and no side vent, as the model without the side vent is between 6 K and 9 K higher for much of the GH. This also shows that the temperature decreases as the roof vent moves further away from the ridge, for the models with no side vents, as the model with roof vents furthest away from the ridge has a temperature distribution very similar to that of the initial model.

The pattern of the water vapour mass fraction contour is very similar for the different models analysed. The contours tend to the air-flow patterns, as the water vapour builds up in areas where the air velocity is lower, so the air renewal rate is also lower. The model shown in Figure 14 has a lower temperature than the model with the roof vents 2 m from the ridge but has a higher water vapour contour for the different air-flow patterns described above. This suggests that if the air can stagnate above the roof vents, as it does in Figure 14, then the percentage of water vapour in the air can increase, even though the temperature of the air is lower. This would increase the GH humidity level.

4.5. Results of the Sensitivity of the Wind Speed

The external wind speed does not have much effect on the air-flow patterns in the GH for this model, which is apparent when comparing Figure 15 with Figure 3. The air-flow is generally the same pattern as described in the initial model. This is expected to be the case, as for all the external wind speeds analysed, the buoyancy effect should still be the driving force, so changing the wind speed would have little impact on the air-flow distribution. However, the air velocity does increase as the external wind increases. The air velocity at the windward vent given as a ratio of the external wind speed was analysed in Figure 4 for the four simulations studied. Figure 16 shows the ratio decreasing until it reaches the value of one, as the external wind speed would increase until it is the dominant effect, at which point pressure differences caused by the buoyancy effect should have a negligible effect on the average wind velocity through the side vent, thus producing a wind velocity through the side vent that is equal to the external wind velocity. However, not enough analysis of higher wind speeds was performed to state with certainty that this is the case.

The temperature contours were also very similar to the original model. This also proves that the wind speed has very little effect on the ventilation of the GH at low wind speeds, when the buoyancy effect drives the ventilation mechanism. The water vapour level gradually decreases as the wind speed increases, the most likely cause of which is that the higher wind speeds can more efficiently ventilate the areas of the GH where the air is slower-moving. This causes the air renewal rate to increase in these areas, thus meaning less water vapour per unit volume of air.

4.6. Results of the Sensitivity of the Floor Heating Analysis

The analysis of the floor temperature on the air velocity and pattern showed that there was very little effect on either variable. The air-flow patterns that formed were almost identical to those of the initial model, as seen in Figure 3. This is likely due to the small changes in temperature that occur in the GH, which would not have much of an effect on the buoyancy forces acting on the air molecules.

The temperature and humidity contours are also qualitatively very similar in each of the models, for the same reasons as given above. However, the temperature does quantitatively increase with the increase in floor temperature. Figure 17 compares the temperatures at the vertical displacements above the floor, for the models with floor temperatures of 310 K, 315 K, 320 K and 325 K. The temperatures vary slightly for most of the GH, as increasing the floor temperature by 5 K increases the average temperature by around 1 K. However, Figure 18 shows that the temperature distributions around the GH are the same.

5. Conclusions

The initial 2D model designed in this study was validated by the work performed experimentally by Sase et al. [35] and using CFD methods by Mistrotis et al. [20]. The air-flow pattern and temperature contour were both found to be quantitatively very similar to both studies mentioned. Additionally, the humidity was found to be highest under the windward vent, where the air renewal rate was lower than in the rest of the GH. The initial 3D model showed a very similar flow pattern to the 2D model, where most studies found that for a wind direction perpendicular to the GH ridge, the variables were almost independent of the GH end walls. The humidity was also found to be higher in the same place than it was in the 2D model.

The geometrical analysis found that the position of the side vents caused a reversal of the direction of the air loop formed in the GH and found that the temperature was higher at the crop level when the side vents were located higher up the wall of the GH. The humidity was also found to increase with higher side vent positions. A similar conclusion was made by Bournet et al. [39,67] for multi-span GHs, which was that the high side vent provided a lower air renewal rate and greater heterogeneity of temperature at the plant level. The importance of the side vent was also analysed, as temperatures in the GH could increase by almost 10 K if no side vent was designed. The humidity was also much higher in the GHs with no side vents.

The roof vent position was found to have less of an effect on the ventilation, air-flow pattern and temperature at the plant level. However, it was shown to have a much greater effect on the ventilation in the roof, as the air became stagnant above the height of the roof vent, causing higher temperatures. The same conclusion was also made when the roof vents were removed, as the temperatures increased rapidly in the roof. The position of the roof vent had no effect on the humidity. When the side vents were removed, the roof vent position was more influential, as the average temperature decreased when the roof vent was moved further from the ridge.

The external wind speed was shown to have very little effect on the air-flow pattern, temperature and humidity for the wind speeds studied, as the buoyancy effect was the driving ventilation force in all simulations tested. The temperature was shown to increase with increased floor temperature, but this also had no impact on either the air-flow pattern or humidity.

Therefore, for the MENA countries, several side vents should be designed to open in different combinations when the temperature or the relative humidity reaches a certain amount or for different environmental conditions. The area which hosted the plants could be designed as a cavity to preserve their microclimate conditions, while the extra humidity could get condensed using dehumidification techniques. For the countries having hot climate conditions, one or multiple roof vents should be designed to release the extra heat for the intense conditions, when necessary.