Received 10 March 2016; accepted 3 May 2016; published 6 May 2016

1. Introduction

The use of canary plastic greenhouses type has been spreading along the Atlantic coast of Morocco, especially in souss valley [1] . In the climate conditions prevailing in this region, efficient climate control is the major challenge that growers have to face. Moreover, the resulting microclimate in this region is far from satisfactory for the crop during a large part of the year.

In order to reduce the need for pesticide application, the greenhouse roofs and sidewall ventilation openings were equipped with fine insect screens which act as mechanical barriers to insect entry inside the greenhouse. However, the use of insect screens reduces the greenhouse air renewal and tends to increase the inside air temperature and humidity, and therefore nearly precludes production.

Consequently, the performance and control of air exchange rate are crucial factors for crop growth during the hot and cold periods of year and for improved production yield and quality.

The ventilation process is the driving force for the greenhouse airflow and climate distribution. But up to now, very few studies were performed on canary greenhouse, equipped with very fine mesh screen, ventilation.

Demrati [2] has used the energy balance method to predict the greenhouse ventilation rate of a banana crop canary greenhouse. Some authors used computer programs to investigate the influence of insect screens and predict the inside microclimate distribution, Fatnassi [3] [4] have characterized and developed a three dimensional CFD simulation model of the inside climate distribution of 0.5 ha canary type greenhouse, equipped with insect screens (6 × 6) on side and roof ventilation openings.

Also, the effect of wind speed and direction and the ventilation openings design on natural ventilation of a Spanish “Parral” greenhouse equipped with insect proof was analyzed by Campen [5] and Molina-Aiz [6] .

Demrati [7] has developed and tested models to enable banana crop transpiration flux to be determined, in a real-canary greenhouse conditions, the transpiration model of the banana crop under cover was validated with respect to the estimated rates of heat and mass transfer implied from indoor and outdoor climatic measurements.

Result obtained by Majdoubi [1] showed that the relatively “bad” ventilation performance of canary greenhouse is not a result of the low value of the greenhouse wind-related ventilation efficiency coefficient (C_w) but of the low value of the discharge coefficient (C_d), caused by a high pressure drop along the air circulation generated both by the use of fine insect screens with small openings and obstruction due to the crop row orientation which were perpendicular to the prevailing wind direction.

As a previous CFD study Majdoubi [8] analyzed the air circulation and the distributed climate during daytime in a 1-ha Canary type tomato greenhouse in the coastal area of southern Morocco. Analyses show that even with low outside wind speed, the outside wind governs in turn the inside airflow direction, inducing a strong windwise air current above the canopy and a very slow reverse flow inside the crop canopy. The weak air exchange within the canopy governs the climate at this level, with a major increase in air temperature and a more moderate increase in specific humidity.

Very few studies have been conducted on nocturnal canary greenhouse climate distribution. Therefore, the main objectives of the present research were to investigate and predict in detailsthe nocturnal distributed climate, as the airflow, temperature and humidity fields, of a real large scale Moroccan Canary greenhouse (1.125 ha area) equipped with insect screens (20/10) on the roof and sidewalls ventilation openings. Simulation conducted by using a commercial CFD package (CFD2000®/Storm), a realistic model based on energy, mass and momentum exchanges was considered. As the regime in greenhouse is turbulent, the turbulent transfers were described by a k-ε model. Likewise, the dynamic influences of insect screens and tomato crop on airflow movement were modeled by means of the concept of porous medium with the Boussinesq assumption. Atmospheric radiations contribution was included in the model by customising the roof cover temperature deducted from its energy balance by means of view factor law. Also, the CFD code was customized in order to simulate in each element of the crop cover the sensible and latent heat exchanges between tomato crop and greenhouse.

The experimental data of greenhouse temperature and humidity distribution were measured and were used to validate this model, the external climate conditions being considered as boundary conditions for the simulation.

2. Computational Fluid Dynamics Theory

The velocity field U, and the associated temperature field T or water vapor content w can be deduced at any time from the resolution of the mass, momentum and energy balances equations, CFD2000® manual [9] :

(1)

where is the studied variable, i.e. the 3 components of the speed vector, T the temperature or a given mass component such as the air humidity content w, is the diffusion coefficient of the quantity, the source term and the speed component. The governing equations are discretized following the procedure descripted by Patankar [10] . This consists of integrating the governing equations over a control volume. Theses equations are numerically solved using a finite volume code CFD2000 using the PISO algorithm developed by Issa [11] . This code also allows for the modeling of the turbulent constraints by means of the standard k-e turbulence model [12] .

2.1. Flow through a Porous Medium

In our CFD model, the insect screens and the crop cover were simulated as porous medium and the air flows governed by the Darcy-Forchheimer Equation (2):

(2)

where U is the air speed, m is the dynamic viscosity of the fluid, K is the permeability of the porous medium and is the non-linear momentum loss coefficient. The values of the aerodynamic proprieties (K and C_F) of the porous medium (screen) were calculated using relations from literature, which correlate these proprieties with the porosity [13] .

(3)

(4)

where is screen porosity, it can be deduced from the dimensions of the thread [13] :

(5)

where d, L and are mesh’s wire diameter, length and width respectively.

For the crop, this sink of momentum is proportional to the leaf density and it may be expressed by unit volume of the cover by the commonly used formula [14] :

(6)

where LAI is the leaf area index and C_D is the drag coefficient of the vegetal cover. For a mature greenhouse tomato crop, Haxaire [15] found, using wind tunnel facilities.

For the crop and the range of air speed observed into the cover, the term in U of Equation (2) can be neglected in front of the quadratic term and the non-linear momentum loss coefficient C_F and the permeability K of the medium can be deduced from the combination of Equations (2) and (6):

(7)

For our simulations, the tomato crop cover was assimilated to a unique 2.6 m high parallelepiped block of porous medium with the same length and width than the greenhouse and with a leaf area index, LAI, equal to 3, the drag coefficient C_D being equal to 0.32 [15] .

2.2. Meshes and Boundary Conditions

The limits of the computational domain include the greenhouse and the free space, high (38 m), windward (30 m), leeward the near greenhouse (surrounding greenhouse 90 m) and on the lateral sides (2 × 30 m) of the greenhouse. The computational grid used Cartesian body fitted coordinates and after several trials with different densities, the calculations were based on a 192 by 44 by 112 grid with finer resolutions were imposed near soil, walls and roof, due to stronger thermal gradients [8] . The boundary conditions prescribed a null pressure gradient in the air, at the limit of the computational domain. The outside air speed was perpendicular to the West- East sidewalls ventilation openings with a measured value equal 1.05 m∙s⁻¹. The driving force of natural convection is the wind force and buoyancy force arising from small temperature differences within the flow according to the boussinesq hypothesis. On average 10 days were necessary before converging to a solution for each trial using a 2.5 GHz frequency computer with 512 MB random access memory (RAM) (Figure 1).

The averages and standard deviations of the climate boundary conditions are summarized in Table 1.

2.3. Simulation of Dynamic, Thermal and Hydrous Effects of the Crop Cover

The radiative flux, reaching each mesh of the crop cover was assimilated to a “volumic heat source boundary condition” and partitioned into convective sensible and latent heat flux (water vapor) depending on the heat and water exchanges between air and the virtual solid matrix representing the crop cover and characterized by its surface temperature:

(8)

The sensible heat flux can be expressed with respect to the difference of temperature between inside air and the canopy:

(9)

and the latent heat flux deduced from a similar relation :

(10)

where and r are respectively the specific heat of air at constant pressure (J∙kg⁻¹∙K⁻¹) and the air density (kg∙m⁻³), is the greenhouse air temperature (K), is the vegetation temperature (K), the aerodynamic resistance between the leaf and the air within each mesh (s∙m⁻¹), is the latent heat of vaporisation of water (J∙kg⁻¹), is the saturated water content of the air at vegetation temperature (kg∙kg⁻¹), is the air water vapor content (kg∙kg⁻¹), is the leaf stomatal resistance (s∙m⁻¹) and the Lewis number.

In order to take into account these new exchanges, the CFD code was customized by means of a “source model” [15] :

(11)

where the different terms of this equation were identified with the terms of the sensible and latent heat transfer equations between the crop cover and the air within each mesh, i.e. for the temperature :

and

and for air humidity:

and

The aerodynamic resistance was deduced from the air speed within each mesh:

(12)

where is the characteristic length of the leaf (m), the interior air speed (m∙s⁻¹), is the air thermal conductivity (W∙m⁻¹∙K⁻¹) and is the air viscosity.

The radiative flux was considered as not limiting and the tomato leaf stomatal resistance deduced from air temperature and saturation deficit [16] :

(13)

where: and.

Tomato crop temperature and humidity can be customized and calculated in simulation model according to the following equations [8] :

(14)

(15)

With LAI is the crop stand leaf area index and R(z) is the global radiation inside the greenhouse in W∙m⁻².

2.4. Simulation of the Contribution of Solar and Atmospheric Thermal Radiation

Atmospheric radiations contributions were included by customising in the model the roof cover temperature deducted from its energy balance by means the view factor law [8] .

(16)

with

and

knowing the temperatures of the greenhouse system and his optical and thermal properties we can then expressed the plastic cover temperature from inside and outside air mean temperatures. The temperature is introduced in the model by using the option “user defined” available in the menu “Wall”.

3. Materials and Methods

3.1. The Greenhouse

The experimental greenhouse is a large scale commercial canary plastic greenhouse covered with a 200 µm thickness single layer polyethylene plastic. Its surface occupies 1.125 ha (90 m length and 125 m width) with a height of 5 m at the gutter and 5.5 m at the ridge and was surrounded by similar greenhouses. The orientation of its spans and tomato crop rows was North-South, i.e., perpendicular to the direction of the prevailing sea wind.

Greenhouse ventilation was performed by means of seventeen roof ventilation openings (0.6 × 125 m² each, total of 1.275 m²) covered with insect screens (20 meshes∙cm⁻¹ in width, 10 meshes∙cm⁻¹ in length with a wire diameter of 0.28 mm). The sidewalls ventilation openings were equipped with similar insect screens. The maximal opening areas reached 875 m² for the West-East sides and 630 m² for the North-South sides. During the experiment, the roof and sidewalls openings were maintained unchanged and equal to 1275 m² and 1505 m² respectively (Figure 2).

The greenhouse was occupied by a tomato crop (Solanum Lycopersicum, cv. Gabriella) planted with a plant density of 1.8 plant∙m⁻² and north-south oriented rows, i.e. perpendicular to the prevailing see breeze direction coming from west. Irrigation was performed by means of a drip irrigation system.

3.2. Experimental Measurements

The following parameters were systematically recorded:

・ The inside net radiation R_net was monitored by means a net radiometer (Q-7, Campbell Scientific Ltd, UK), situated between the top of the crop canopy and the polyethylene film of the roof.

・ The inside conductive heat flux exchange at soil surface FS which was measured by means of a conductive flux meter (HFT3, Campbell Scientific Ltd, UK) situated 1 mm below the soil surface.

・ The inside and outside air temperature T and relative humidity HR were measured by means of thermo- hygrometers probes (HMP45 AC, Vaisala, Etoile Internationale, Paris).

・ The inside soil surface, roof cover, tomato leaf temperatures were measured by means of thermocouples (Copper-Constantan) which were stuck on the plastic or positioned 1mm below the soil surface together with fine thermocouples which were inserted in the principal vein of the terminal leaflet on the lower face of a leaf.

・ The outside wind speed U and direction Dv were measured by means of a cup anemometer and a wind vane (W200P, Campbell Scientific Ltd, UK) located 8.5 m high on a mast situated 3 m over the greenhouse ridge (A 100R. Campbell Scientific Ltd, UK).

・ The outside global radiation Rgo was measured by means a pyranometer (SP-LITE, Kipp & Zonen, Campbell Scientific Ltd, UK).

・ All these measurements have been sampled every 5 seconds and the data were then averaged and stored every fifteen minutes into two data loggers (Models 21 X and CR23, Campbell Scientific, Inc., Logan, UT).

・ The parameters described in Table 1 were systematically recorded, in order to determine the simulation boundary conditions and to validate the simulation model.

4. Model Validation

4.1. Distributed Inside Climate

Figure 3 shows CFD predicted and measured air temperature profiles from west to east in the middle of the greenhouse at 1 and 4 m high.

It can be seen that the temperature difference between numerical results and those measured in the experimental greenhouse was very weaker. This difference ranges between 0.1 and 0.6˚C. The average difference between the computational and experimental data was 0.4˚C. The same figure shows that the greenhouse air

temperature is almost homogenous along the west-east direction contrarily to result found by majdoubi [8] during daytime.

The specific air humidity measurements inside experimental greenhouse and those simulated by using CFD2000 models are shown in Figure 4, at 1 and 4 m high above soil surface.

It can be seen that values for the measured and predicted results were very close. The differences between specific air humidity measured and those simulated ranges between 0 and 0.5 ´ 10⁻³ kg∙kg⁻¹. The average difference value reaches 0.13 ´ 10⁻³ kg∙kg⁻¹. Also we can observe that air specific humidity have the same value along the west-east direction.

4.2. Verification Respect to the Global Energy Balance

Estimation of the exchange rate G:

Integrating air speed over a complete cross-section of the greenhouse ventilation opening allows for estimation of air flow rate G, as follows [6] :

(17)

which

represents the opening length (m), U is the horizontal air velocity (m∙s⁻¹), is the opening height (m) and z is the vertical distance.

In our case for an external air velocity equal to 1.05 m∙s⁻¹, airflow rate was estimated as G = 93.5 m³∙s⁻¹ which allows 5.6 air renewals per hour [17] (Table 2).

Estimation of the net radiation R_net:

The energy balance of our greenhouse can be written as [2] :

(18)

With

- F_S the soil surface flux

- the sensible heat flux extracted by ventilation.

- the latent heat flux extracted by ventilation.

- the sensible heat flux exchanged by convection between inside air and the greenhouse cover

- the overall sensible heat loss through the side walls, excluding the side ventilation openings

From Table 2 which summarizes the calculated and measured values of all terms signaled before in the energy balance Equation (18), we can show a good agreement between measured and estimated net radiation values, , with an error about 3% [17] .

Verification respect to the average vegetation transpiration flux:

We have plotted a vertical cross-section of the vegetation transpiration flux field across the middle of the greenhouse (Figure 5). This figure shows that the transpiration flux value ranges between 2 (W∙m⁻²) at the top of canopy and 5 (W∙m⁻²) near the soil surface. The average value reaches 3 (W∙m⁻²). These results confirm the absence of crops transpiration during night time. This value is very low like as the calculated value of the latent heat flux exchanges between greenhouse inside and outside air which was very weaker (6 ´ 10⁻³ W∙m⁻²), Table 2.

5. Detailed Analysis of Inside Climate and Air Movement Analysis

Details of air wind speed, temperature and specific humidity patterns in vertical or longitudinal profiles of the studied greenhouse were presented (Figures 6-9).

5.1. Airflow Circulation

For a wind speed perpendicular to the roof and the west-east sidewalls ventilation openings surface. Figure 6 illustrates a vertical cross-section along the flow direction in the middle of the greenhouse. Generally, wind speed is important in the top of canopy with a similar direction to outside wind from west to east. Contrary, on the level of crop cover, air velocity is much lower, together with the development of a reverse flow from the leeward end to windward end part of the greenhouse. Similar observations were reported by Majdoubi [8] .

The air velocity profiles across the middle of the greenhouse at heights of 1, 3 and 4 m were reported in Figure 7.

The results confirm that the airflow in the high of canopy, especially between 0 and 2 m, circulates in the

opposed direction to the prevailing outside wind “from east to west” and is even much lower (1/25 of outside wind velocity U_ext). This strong air decrease mainly was caused by the drag effect due to the high crop cover density and insect screens resistance to the airflow, which can reduce the greenhouse ventilation rate by 50% and 46% respectively [1] . On the other hand, the change of direction is due mainly to the air penetration by roofs ventilation openings at the gutter (5 m). Also, the same figure, show that the air velocity at 3 and 4 m high, between the top of the tomato crop and the roof cover, within the experimental greenhouse in the leeward end is higher than that of the windward end, it is ranges between 0.16 U_ext and 0.4 U_ext respectively, result which can confirms that greenhouse air velocity is governed mainly by roof ventilation openings at the gutter. It appears clearly that the roof openings are alternatively acted as air entrances and exits. Which confirm the results observed byMajdoubi [8] .

5.2. Temperature and Humidity Pattern

From Figure 8 which represent the vertical profile of simulated temperature in the center of the greenhouse respect to the height. we can show from the previous figure that the temperature is high at the level of the soil surface then it decreases with the height to 4 m approximately, Within the crop canopy air temperature rises from 291 K (19˚C) at its top to 293 K (20˚C) near the soil surface. In continuation one observes a slight decrease in temperature which reaches its minimum on the level of the plastic roof cover 290 K (17˚C). Contrarily, in the daytime when the temperature reaches its maximum on the level of plastic cove, which explains the presence of thermal inversions during night time when plastic cover significantly cooler the greenhouse inside air, the same results have been recorded by Montero [18] for a three-span sloped-roof greenhouse in the south of Spain with a single layer of polyethylene cladding material. Also, Majdoubi’s [19] experimental data show that the canary greenhouse air temperature by night is lower than outside, this difference is approximately equal to 1.5˚C. This temperature inversion is due to the radiative cooling of the greenhouse promoted the confinement of the shelter and the high transmission of the plastic cover to the infra-red radiations. At the top of the canopy, there is no significant difference between the greenhouse air temperature and outside.

The air Specific humidity distribution vertical profile was reported in Figure 9. The figure shows clearly that the greenhouse air specific humidity above the canopy is similar to the air outside. But it is lower than that at the level of the canopy, especially near soil surface where the maximum was observed. General, specific humidity distribution was homogenous and there is not significant difference between inside and outside air specific humidity.

6. Conclusions

A three dimension model was developed to predict the nocturnal airflow and thermal and hydrous field distribution in a large scale canary greenhouse (1.125 ha) planted with a tomato crop and protected by insect screen (20/10) by means of CFD software. The model numerically solves the mass, momentum and energy conservation equations after customization of the software to also simulate atmosphere radiation and plastic cover temperature and water vapour and thermal exchanges between the crop cover and the inside air. The mechanisms induced by the presence of insect screens and the drag effects of the crop cover were also satisfactorily simulated by using the porous medium approach.

The model was first validated by checking numerical results with respect to the measured or calculated ones, as air temperature and air specific humidity values, the net radiation measured by using the calculated ventilation rate, and the transpiration flux calculated value always by using the air exchange rate, stresses the relatively good accuracy of the model and its reliability for simulation.

Results show that air velocity is important in the top of crop cover with the same direction respect to external prevailing wind, but it is higher in the leeward end than in the windward end. That means that air wind speed is governed by penetrating of cold air by roof ventilation openings at the gutter. Also, there is a strong decrease of air speed in the height of canopy with an opposite direction respect to external airflow circulation. It also shows that inside air temperature and humidity fields were homogeneous along the greenhouse area and there is no significant difference between inside and outside air temperature and humidity. So, there is the necessity to use an internal curtain or/and close the roofs or sidewall ventilation openings in order to limit the energy drops during night time.

Acknowledgements

The authors are highly thankful to the AUF “Agence Universitaire de la Francophonie” and the project PROTARS III “Programme Thématiquesd’Appui à la Recherche Scientifiques” CNRST Morocco, for the financial support during this research work at experimental and simulation period.