CFD Simulation of the Increase of the Roof Ventilation Area in a Traditional Colombian Greenhouse: Effect on Air Flow Patterns and Thermal Behavior

CFD Simulation of the Increase of the Roof Ventilation Area in a Traditional Colombian Greenhouse: Effect on Air Flow Patterns and Thermal Behavior

Edwin Andres Villagran MunarCarlos Ricardo Bojacá Aldana

Department of Basic Sciences and Modelling, Faculty of Natural Sciences and Engineering, Jorge Tadeo Lozano University, Bogotá 110311, Colombia

Corresponding Author Email: 
edwina.villagranm@utadeo.edu.co
Page: 
881-892
|
DOI: 
https://doi.org/10.18280/ijht.370326
Received: 
9 July 2019
|
Revised: 
24 August 2019
|
Accepted: 
5 September 2019
|
Available online: 
28 September 2019
| Citation

© 2019 IIETA. This article is published by IIETA and is licensed under the CC BY 4.0 license (http://creativecommons.org/licenses/by/4.0/).

OPEN ACCESS

Abstract: 

The protected horticultural production in tropical countries such as Colombia takes place mainly under passively ventilated greenhouses. The major drawback of the traditional Colombian greenhouse is its low ventilation rate leading to inadequate microclimates. The purpose of this study was to evaluate the airflow behavior, ventilation rates and thermal distribution effects of a commercial-size greenhouse, due to the gradual increase of the roof ventilation areas. A numerical simulation model was developed under a computational fluid dynamics (CFD) approach following an experimental validation of the simulated scenarios. The results showed that increasing the roof ventilation areas reduced the internal average temperatures by up to 16 %, improved the temperature homogeneity and increased the ventilation rates by up to 191 % with respect to the reference scenario. The findings of this research may serve as a basis to develop alternative greenhouse designs that consider increasing the ventilation areas but within the low-tech context that constrains innovation in the Colombian horticultural sector.

Keywords: 

finite volumes, passive ventilation, protected horticulture, thermal difference, ventilation rate

1. Introduction

Colombia is the largest producer of cut flowers in the South American continent and ranks second worldwide in export volumes after the Netherlands [1]. Ornamental production is developed mainly in the cool highlands of the Colombian Andes, with 78 % of cultivated area (ca. 7800 hectares) established in the Bogota plateau [2]. Since its establishment, more than 60 years ago, the cut flower sector promoted the use of the traditional Colombian greenhouse (CTG), which is a low-tech low-cost greenhouse model, built in wood or metal and covered with plastic. Still today, despite its limitations and low ventilation rates, the CTG greenhouse is widely established in approximately 70 % of the production areas [3]. This greenhouse model has also been used for tomato cropping across several production areas of the country [4]. The main issue with the CTG greenhouse is its restricted airflow leading to high temperatures or humidity levels during critical periods along the day. This inadequate microclimate results in strong incidence of diseases that can cause relevant production losses such as downy mildew and grey mold, caused by Peronospora sparsa and Botrytis cinereal, respectively.

This type of greenhouse is characterized as a passive greenhouse, meaning that the climate control method is the natural ventilation, which operates through openings located on the side and roof areas of the structure [5]. The aims of this ventilation method are to regulate the internal thermal and hygrometric behavior, and to maintain the greenhouse CO2 level close to that present in the exterior. These climate parameters affect the optimal development of crops if they are not within certain range [6].

The driving forces of the natural ventilation comprises a static component related to the wind effect, created by the speed of the outside air, and a dynamic component generated by the so-called buoyancy effect. The latter effect is the result of changes in the internal air density as a consequence of increasing temperatures and turbulent flows due to pressure changes in the cross-sectional and longitudinal greenhouse areas [7]. This phenomenon has been extensively studied since the second half of the 20th century [8]. During this period, various methodologies have been developed to analyze natural ventilation including modeling and simulation techniques as well as real and reduced scale experimentation, as depicted by Bournet and Boulard [9].

Natural ventilation is highly dependent on external meteorological variables such as wind speed, intensity and direction, greenhouse design factors such as roof geometry, greenhouse area and volume, and size and location of ventilation areas. According to the above, natural ventilation varies over daytime generating sometimes an inadequate microclimate condition inside the greenhouse [10, 11]. The CTG greenhouse is generally located in regions with low wind speeds and its ventilation areas are less than the recommended figure of being at least 25 % of the covered ground area [12]. These two characteristics of the Colombian greenhouse production system compromise the performance of the CTG greenhouse, resulting in low ventilation rates that rarely exceed the minimum value of 0.04 m3s-1m2 established for passive greenhouses [13]. These low ventilation rates lead to the occurrence of inadequate microclimatic conditions characterized by heat and humidity accumulation and, a significant reduction in the internal CO2 levels. This situation occurs mainly during high radiation moments, that in tropical regions, happens for most of the daylight hours (9:00-16:00).

One of the most extensively modeling and simulation technique applied on aerodynamic and microclimate studies in greenhouses has been computational fluid dynamics (CFD), which is known as a robust, mature and advanced design technique in many engineering fields. CFD modeling and simulation offer fast numerical solutions to real or prospective scenarios considering structural modifications in the case of buildings. Being a computer-aided simulation tool, CFD has also become a tool to design and optimize greenhouse ventilation systems [14-17]. The use of CFD has contributed to the optimization of greenhouse ventilation systems in countries such as Spain [15, 18], Italy [19], China [20, 21] and Mexico [22].

An adequate comprehension of the natural ventilation phenomenon is essential to optimize or redesign the ventilation systems of greenhouses. Therefore, the objective of this work was to evaluate the gradual increase of the roof ventilation areas and its effect on the thermal behavior and airflow patterns of a CTG greenhouse using a 2D CFD model. The numerical model considered the main meteorological variables of the Bogota plateau and a virtual geometric representation of the studied greenhouse.

2. Materials and Methods

2.1 Description of the greenhouse and climatic conditions

The experiment was developed on a commercial size CTG devoted to cut flowers production, with wooden structure and polyethylene cover located in the Bogota plateau (5º03’31.81’’ N, 73º55’49.22 ’W, 2571 m.a.s.l.). The covered area was 5353 m2, equivalent to a greenhouse composed by 12 spans, 6.8 m wide (81.6 m) and 65.6 m long each, with its longitudinal axis located in the NW-SE direction. The greenhouse was equipped with side vents on its four side walls, with an effective opening area of 1.7 m forming a side ventilation area of 500 m2 (9.3 % of the covered area) and a top ventilation of 0.19 m opening in 11 of the spans, for a roof ventilation area of 137.1 m2 (2.56 % of the covered area).

As can be seen in Figure 1a, the region presents climatological characteristics where the mean multiannual average temperature over a 15-year period is 13.19 °C, with an average of maximums and minimums of 21.90 and 6.21 °C and high humidity conditions generally higher than 78 %. The annual precipitation reaches a value of 830.7 mm, exhibiting a bimodal behavior with maximum values grouped in the months of March-May and September-November. The analyzed wind speeds were grouped in 5 frequency intervals, with 84.5 % of total data located in the first 3 intervals which present speeds lower than 1 ms-1 (Figure 1b). The predominant wind directions were found in the S-W quadrant (180-270°).

a) Multiannual climograph

b) Wind speed distribution frequencies

Figure 1. Meteorological characteristics of the study area.

2.2 Model and numerical procedure

Numerical simulations solved the equations that describe the flow behavior of a fluid between the exterior and interior environment of a greenhouse, using the second-order finite volume discretization (FVM) method. The equations of moment, mass, energy and concentration can be represented for a two-dimensional flow in steady state with the following equation:

$\frac{\partial \phi}{\partial t}+\frac{\partial}{\partial x_{j}}\left(u_{j} \phi\right)=\frac{\partial}{\partial x_{j}}\left(\Gamma_{\phi} \frac{\partial \phi}{\partial x_{j}}\right)+S_{\phi}$       (1)

where, ϕ represents the variable of interest, in this case the vertical and horizontal components of speed uj (ms-1) and temperature T (°C). Sϕ y Γϕ represent the source term and the diffusion coefficient of ϕ [23].

The numerical model considers the turbulence of the flow through the standard model k - ε [24]. This semi-empirical model is based on the transport equations that solve turbulent kinetic energy k and the dispation of this energy per unit volume ε, this model has been the most used and widely validated in studies of air flow in greenhouses showing to be efficient with the use of computational resources and providing realistic solutions. [25]. The transport equations for k and ε can be modeled as:

$\frac{\partial}{\partial x}(\rho k)=\frac{\partial}{\partial x_{j}}\left[\left(\mu+\frac{\partial k}{x_{j}}\right) \frac{\partial k}{\partial x_{j}}\right]+G_{k}+G_{b}-\rho \epsilon-Y_{M}$        (2)

$\frac{\partial}{\partial t}(\rho \varepsilon)=\frac{\partial}{\partial x_{i}}\left[\left(\mu+\frac{\mu_{t}}{\sigma}\right) \frac{\partial \epsilon}{\partial x_{i}}\right]+\rho C_{1} S_{\epsilon}-\rho C_{2} \frac{\epsilon^{2}}{k+\sqrt{v \epsilon}}+C_{1 \epsilon} \frac{\epsilon}{k} C_{3 \epsilon} G_{b} k$       (3)

where, $\mu$ is viscosity and $\mu_{t}$ is turbulent viscosity in (kg m-1s-1), Gb is turbulent kinetic energy generation due to buoyancy, Gk is turbulent kinetic energy generation due to velocity gradients, $\sigma_{k} y \sigma_{\epsilon}$  are Prandtl's turbulent numbers for $k$ y $\varepsilon, v$  is the coefficient of kinematic viscosity, YM  is the fluctuating expansion in turbulence due to the global dissipation rate and $C_{1 \epsilon}, C_{2 \epsilon}, C_{\mu}, \sigma_{k} y \sigma_{\epsilon}$  are constant with empirically determined values and established by default in the simulation software [26].

To take into account the buoyancy forces the energy equations and the buoyancy model were activated by the Boussinesq approximation applied to the entire computational domain, this approximation provides faster convergence, compared to numerical simulations where the density variable is considered in all equations, this model takes gravity into account and is added to the momentum equation as a source term that models changes in air density due to increases in temperature value. Roy, et al., and Camacho, et al., [27, 28] are described by the following equation:

$\rho=\rho_{0}\left(1-\beta\left(T-T_{0}\right)\right)$         (4)

where, $\beta$  is the coefficient of thermal expansion of the air (°C), $T$ the temperature at the instant analyzed $\left(^{\circ} \mathrm{C}\right), \rho_{0} y T_{0}$ are the reference values of the density (Kg m-3) and the temperature (°C) respectively, this approximation can be applied under the following conditions, the temperature differentials in the flow field of the computational domain are less than 30 °C, the density differentials are required only by the buoyancy term of the moment equation and there is a linear relationship between temperature and density, and all other extensive fluid properties are constant [29, 30].

The semi-implicit solution method for the pressure-velocity equation (SIMPLE) was applied to solve the flow field of the simulated fluid, with second order discretization schemes, this method of solution has a working scheme which starts with some initial values to solve the moment equations and obtain the components of the velocity variable, successively calculates coefficient and source terms through the pressure correction equation, calculates the velocity corrections and updates the pressure and velocity components to obtain the solution. The model convergence criteria were established with values for the residuals of 10-6 for the energy equation and 10-4 for the continuity, momentum and turbulence equations, values that have been shown to be accurate in other similar studies [31].

2.3 Computational domain and mesh size

The two-dimensional (2D) computational domain and the generation of the mesh was done through the use of commercial pre-processing software ANSYS-ICEM (v.17.2). The computational domain included in detail the geometrical characteristics of the middle cross section of the greenhouse and a size large enough to allow adequate airflow development and a coupling between the exterior and interior environment of the greenhouse. This domain was built following the recommendations included in the wind effect study of Tominaga et al. [32] and successfully applied in numerical studies such as that of Peren et al. [33]. The selected domain dimensions were 381 m for the length (x-axis) and 60 m for the height (y-axis). Subsequently the computational domain was divided into an unstructured mesh with finer resolutions imposed near the floor, ceiling and walls of the greenhouse, where thermal gradients can be more pronounced. The mesh was configured with square elements composed of a total of 551,935, discrete volumes in space. This size was defined once a mesh refining analysis was carried out where a total of 7 squares were analyzed with a quantity of elements that varied between 102,871 and 1,254,329, this to obtain an adequate precision and a total independence of the results to the size of the mesh.

The quality of the mesh was evaluated by means of the relative determinant 2 X 2, which showed a behavior superior to 93.72 % with cells between 0.9 and 1 which indicates a perfectly regular mesh element, additionally the asymmetry factor was evaluated, where 90.8 % of the active volumes exhibited a value within the range 0- 0.25, range considered of excellent quality [34, 35].

2.4 Boundary conditions

The numerical model considers the atmospheric characteristics of the experimental site, such as gravity force, air viscosity and atmospheric pressure. The symmetry limit condition was imposed on the upper limit of the computational domain in order to describe the gradients of all variables studied and the normal velocity to a symmetrical plane. In the solid regions of the computational domain such as the soil and the inverandero structure, the contour condition of the non-slip wall type was established, as well as the properties of polyethylene and agricultural soil, such as density (ρ), thermal conductivity (k) and specific heat (Cp), were established, according to the study developed by Villagran et al. [36]. The effect of solar radiation was not established to be solved by some resolution model for the equation of radiative transfer, although the effect of it was included in the numerical model indirectly by establishing the surface temperatures in the greenhouse walls and on the ground surface inside and outside the greenhouse, following the procedure described by Flores [37], where a heat flow is established on the surface of the soil depending on the incident solar radiation, the transmission of the covering material and the textural class of the soil, which for this case is clay loam. 

The air inlet condition for this study varies according to NE-SW (S1-S7) or SW-NE (S8-S14) wind directions, directions that impact perpendicularly to the side wall and roof ventilation openings, were established as a logarithmic wind profile coupled to the CFD numerical model by means of a user defined function (UDF) [31, 38]. Calculated by the following equation of law for the logarithmic wind profile.

$V_{2}=V_{1} \frac{\left(\frac{h_{2}}{z_{0}}\right)}{\left(\frac{h_{1}}{z_{0}}\right)}$       (5)

where, $V_{1}$ is the wind speed measured in a height $h_{1} ; V_{2}$ is the wind speed calculated for a height $h_{2}$ y $z_{0}$ is the length of the roughness coefficient expressed in meters, value that this in function of the type of surface of the surroundings of the greenhouse, for this case I consider a value of 0.1 that is used for agricultural lands [39].

The air outlet limit depending on the simulation scenario and the corresponding wind direction was set under the pressure output contour condition. The greenhouse was modeled without the presence of plants to obtain results independent of the type, state and disposition of the crop. Outdoor wind conditions were simulated parallel to the cross-section of the greenhouse (x-axis) and perpendicular to the ventilation openings. These two-dimensional modelling approaches have proved to be a precise and agile solution in the study of natural ventilation of protected environments [40].

2.5 Experimental data collection

The CFD model was validated between April 1 and May 5, 2017 with a recording frequency of ten minutes, for which a data capture procedure was established inside the real greenhouse and the external environment. The hourly climatic variables in the external environment of the greenhouse necessary to feed the numerical model such as air temperature, wind speed and direction were recorded by means of a climate system Vantage Pro2 (Davis Instruments, Hayward, CA, EE. UU.). Inside the greenhouse at a height of 1.75 m above ground level (y-axis) and above the cross section (x-axis), a total of 17 sets of two copper thermocouples used to measure dry bulb temperatures were evenly distributed every 5 m and stored in a data logger (Cox-Tracer Junior, Escort DLS, Edison, NJ, EE. UU.).

2.6 Numerical simulations and validation of the numerical model

The data obtained during the measurement period were analyzed, selecting those of a specific time of the 14:00 hour day, which presented a large number of wind direction records in NE-SW direction, this simulation was developed with the average values obtained of temperature, radiation and wind speed which were 21.7 °C, 412 Wm-2 and 0.83 ms-1. Once the simulation was developed, the temperature data generated in the cross section of the greenhouse at a height of 1.75 m above ground level were extracted. These were compared with the data captured experimentally, through a statistical analysis that included the calculation of the following criteria: absolute mean error and error of the mean quadratic root (MAE, RMSE, °C), these are measures commonly used to verify the concordance between measured and simulated values and their interpretation of the value obtained says that the closer to 0 , the greater the degree of fit of the model there is. These measures of goodness-of-fit are calculated by means of the following relationships.

$M A E=\frac{1}{m} \sum_{i=1}^{m}|D m i-D s i|$        (6)

$R M S E=\sqrt{\frac{\left.\sum_{i=1}^{m}|D m i-D s i|\right)}{m}}$        (7)

where, Dmi and Dsi, are the value of the observed data and the value of the simulated data respectively and m the number of data compared.

Once the CFD simulation model was validated, it was considered an agile and accurate tool to investigate the effects of increased ventilation surface on thermal distribution and air flow patterns in the greenhouse, ventilation rates were quantified through the widely used mass flow rate method. [41]. In accordance with the objective of this research, 7 simulation scenarios were proposed, including reference scenario S1. For the remaining S2-S7 scenarios the gradual increase of the roof ventilation surface obtained through the generation of a folding opening of 1.2 meters in the deck areas of the different greenhouse spans was included, as can be seen in Table 1. In addition, in order to evaluate the effect of wind direction, the same S1-S7 scenarios with inverted wind direction were evaluated and in their order they were named S8-S14 (Table 1). The initial conditions to determine the temperature values in the computational domain and the heat source and which were established as constants for each simulated case, were taken from the historical maximum mean values of the study site where a value of 21.9 °C was obtained for the temperature (Figure 1a) and the solar radiation 605 Wm-2, other variables such as wind speed and direction are observed in Table 1.

Table 1. Simulation scenarios evaluated

3. Results and Discussion

3.1 Validation of the model

Experimental and simulated data during the validation of the CFD model can be seen in Figure 2a. The trend on the greenhouse cross axis (x-axis) shows a similar behavior for both data sets. The thermal distribution shown in Figure 2b allows us to observe how the greenhouse is suffering a warming in the direction of the air flow, warming generated due to the wind effect of natural ventilation and the limited roof ventilation surface of the CTG greenhouse. The values of MAE and RMSE obtained were 0.84 and 0.91 °C, values that are equivalent to 3.8 and 4.1 % of the simulated mean temperature and which are values very similar to those obtained in similar studies such as the one developed by Villagran et al., [42]. Therefore, it can be concluded that the predictions of the CFD model are satisfactory and this model can be used to carry out the study for optimizing the air flow and temperature patterns inside the CTG greenhouse.

Figure 2. Measured and simulated temperature patterns

3.2 Airflow patterns and ventilation rates

3.2.1 Wind direction NE-SW

In Figure 3, the behavior of the flow pattern can be observed for scenarios S1, S4 and S7 and for velocities V1, V3 and V5. The average speeds obtained inside the greenhouse for S1V1, S2V3 and S3V5 were 0.13, 0.56 and 2.16 ms-1, which represents a reduction of 35, 44 and 40% with respect to the external wind speed. The characteristic flow pattern for this scenario shows a behavior in which air flows from the windward side window to the leeward side ventilation with a loss of speed from the second span, this may be due to frictional force of some inertial components of the flow pattern. There is also a clear differentiation in S1V1 and S1V3 where there are some movements in the direction of some of the roof vents, this is caused by the chimney effect generated by the difference in temperature between the fresh air coming in from outside and the warm air inside the CTG and which is very relevant under outside wind conditions below 1 ms-1, additionally I can find that although there are these types of movements are deficient and with low wind speeds, which can be generated by the low roof ventilation surface that has CTG. For S4V1, S4V3 and S4V5 mean CTG velocities of 0.15, 0.73 and 2.58 ms-1 were obtained, which represents an increase in velocity compared to the critical scenario S1 of 25, 30.3 and 19.4 % for V1, V3 and V5 respectively.

In the case of S1V1 it can be observed how the air flow pattern presents two movement cells, one from the windward side ventilation to the CTG span 7, where the flow is accelerated by the folding ventilation created in the roof and another flow from the leeward side window directed towards hall 9 where it accelerates horizontally towards the roof sale, this type of behavior is similar to the one found in previous studies such as the one developed by Molina-Aiz et al., [43]. Another common feature in this scenario is the increased airflow movement that occurs in the roof area and the folding vents generated. On the other hand, for scenario S7 a similar behavior to that obtained in S4 is observed, with small improvements in the flows at the ceiling ventilation level and increases in the average air velocities in CTG where values of 0.18, 0.81 and 2.61 ms-1 were obtained for V1, V3 and V7 respectively, these being higher than those obtained in S1 and S4.

The ventilation rates as a function of the covered surface (Φ, m3m-2s-1) for each evaluated configuration were calculated and are shown in Figure 4. For the reference scenario S1 it can be observed that the values of Φ ranged between 0.0085 and 0.067 m3m-2s-1 for V1 and V5, as previously demonstrated in multiple studies, the dependence of the ventilation rate on wind speed is linear, under these results it can be clearly identified that the ventilation rates for the dominant speeds in the study region V1, V2 and V3, are below the recommended minimum value of Φ = 0.04 m3m-2s-1. It is also perceptible that as the roof ventilation area is increased, an increase in ventilation rates is generated for each scenario compared to S1, obtaining ventilation rates with a higher mean value of 52, 78, 97, 118, 150 and 191 % for S2 to S7 respectively, which allows to verify that the ventilation rate is also highly dependent on the ventilation configuration and the arrangement of ventilators used in the greenhouse, as previously demonstrated in several studies and more recently in the one developed by Espinoza et al. [11].

Figure 3. Simulated airflow patterns for S1, S4 and S7

Figure 4. Ventilation rate Φ (m3 m-2 s-1) calculated for each simulated case

3.2.2 Wind direction SW-NE

In general terms the results of air flow and ventilation rates showed to be independent of the direction of the external air flow, this is because the proposed scenarios sought to alternate the directions in which the folding vents were generated in each building, this was done in order to seek the best distribution of the flow pattern and independence to the direction of the wind that is not achieved when these types of windows are installed in a single direction either windward or leeward as can be found in the study conducted by Baeza et al. [44].

Figure 5 shows the flow patterns for S8, S11 and S14 under outside wind speeds V1, V3 and V5. For scenario S8 it is observed that the flow patterns present a displacement behavior that enters CTG through the windward side ventilation and a small flow through the fixed roof ventilation of the second span in contact with the air flow, the mean velocities obtained were 0.14, 0.61 and 2.29 ms-1 being 7.6, 8.9 and 6.0 % with respect to S1. For scenarios S11 and S14 the mean velocity values obtained inside CTG showed increases compared to S1 that ranged between 2.3 and 6.1 % for S11 and 1.2 and 2.3 % for S14.

Figure 5. Simulated airflow patterns for S8, S11 and S14

Figure 6. Ventilation rate Φ (m3 m-2 s-1) calculated for each simulated case

With respect to the calculation of the ventilation rates Φ these can be observed for all scenarios S8-S14 in Figure 6. For the reference scenario S8 the values of Φ ranged between 0.009 and 0.071 m3m-2s-1 for speeds V1 and V5 respectively, this value is 6.8 and 5.5 % higher than S1. In general, this tendency to obtain ventilation rates slightly higher than those of scenarios S1-S7 is observed.  For scenarios S9 to S14 the mean values of Φ were 41, 60, 105, 127, 164 and 223 % with respect to those obtained in S8.

3.3 Temperature behavior and homogeneity

The evaluation of the thermal behaviour of the greenhouse was carried out by calculating the average temperature and the thermal gradient (ΔT) generated between the interior environment of the greenhouse and the exterior environment of the computational domain. Additionally, the homogeneity of the temperature was evaluated by calculating the standard deviations of ΔT (SD-ΔT). This parameter indicates that when the value approaches zero, the most homogeneous is the behavior of the evaluated variable.

3.3.1 Wind direction NE-SW

The qualitative behavior of the temperature variable inside the greenhouse for scenarios S1, S2 and S3, for the respective simulated outdoor wind speeds can be observed in Figure 7. In general terms it was found that there is a relationship between the behavior of the temperature variable, with the outdoor wind speed and the roof ventilation surface, as the values of these two variables increase a gradual reduction of the

temperature value inside the greenhouse is observed. The most critical case was found for S1V1 where it can be observed that up to 5 zones differentiated in thermal value are generated, with a cold zone near the windward side ventilation wall and a warm zone of greater magnitude in the 4 adjacent spans of the leeward side wall, behavior that is repeated although in smaller magnitude for S1V3 and S1V5.

As the roof ventilation surface in CTG increases, it is observed that the temperature distribution in the interior is undergoing qualitative and quantitative changes, finding that for S3V1 and S7V1 the colder zones are in the zones adjacent to the leeward and windward wall while the warmer zones are in the central region of the greenhouse just below spans 8, 9 and 10. On the other hand for S3V3 and S7V3 it is observed that although a larger ventilation surface gradually reduces the value of the temperature, its behavior is similar, finding that the coldest zones are located in the first 3 naves adjacent to the windward wall and the warmest zone in the last 4 naves near the leeward wall.

Figure 7. Simulated temperature distribution contours for S1-S7

The average temperatures calculated from the simulation data at a height of 1.75 m above ground level, showed a decreasing behavior depending on the external wind speed, for a large percentage of the scenarios evaluated, as can be seen in Table 2, except for scenarios S4V2, S5V2, S6V2 and S7V2, where it is observed that the mean temperature value is higher than the values found for V1, this can be strongly influenced by the thermal effect of the natural ventilation that takes on greater relevance at low wind speeds. Thus, the scenarios that presented the greatest magnitude of the temperature variable were S1V1, S2V1 and S3V1 with an average value of 27.51, 26.80, 24.94 °C respectively and S4V2, S5V2, S6V2 and S7V2 with values of 24.20, 23.87, 23.23 and 23.18 °C. On the contrary, the lowest temperature values were found in the scenarios evaluated with V5 velocities, obtaining average temperatures close to 22.2 °C in S5, S6 and S7, although it should be remembered that these high velocities are only present in a percentage lower than 6 % in the region of study.

Comparing scenarios S1 and S7 it was found that the temperature reduction percentage was 16.84, 8.69, 7.04, 3.82 and 2.21 % for V1, V2, V3, V4, V5, V6 and V7, so it is valid to say that for the predominant speeds in the study region V1, V2 and V3 is where most benefit is obtained from increasing the surface area of the roof vent, thus becoming an alternative climate optimization for the CTG type greenhouse, although it is also clear that the relationship of temperature decrease versus roof ventilation surface increase is not linear, since when analyzing the data it was found that the highest percentages of reduction occurred between S2, S3 and S4, the above may be influenced by the ventilation configuration used in CTG, which in this case is combined between the sides and the roof openings, according to the study of the CTG Kittas et al. [45], the greatest benefit of natural ventilation is obtained when the lateral and roof ventilation surfaces are equal, that scenario in our case occurs in the intermediate between S3 and S4.

Table 2. Thermal parameters obtained for simulated scenarios S1-S7

Scenarios

Tmean (°C)

ΔT (°C)

SD-ΔT (°C)

S1V1

27,51

5,61

3,04

S1V2

25,39

3,49

1,98

S1V3

24,63

2,73

1,61

S1V4

23,36

1,46

0,94

S1V5

22,74

0,84

0,55

S2V1

26,80

4,90

3,35

S2V2

25,04

3,14

1,94

S2V3

24,32

2,42

1,61

S2V4

23,27

1,37

0,99

S2V5

22,71

0,81

0,57

S3V1

24,94

3,04

1,89

S3V2

24,28

2,38

1,66

S3V3

24,19

2,29

1,54

S3V4

23,21

1,31

0,99

S3V5

22,69

0,79

0,58

S4V1

23,63

1,73

1,20

S4V2

24,20

2,30

1,47

S4V3

23,63

1,73

1,17

S4V4

22,89

0,99

0,74

S4V5

22,49

0,59

0,44

S5V1

23,38

1,48

1,02

S5V2

23,87

1,97

1,26

S5V3

23,32

1,42

0,95

S5V4

22,71

0,81

0,61

S5V5

22,40

0,50

0,39

S6V1

23,09

1,19

0,81

S6V2

23,23

1,33

0,74

S6V3

23,08

1,18

0,78

S6V4

22,59

0,69

0,51

S6V5

22,31

0,41

0,31

S7V1

22,88

0,98

0,67

S7V2

23,18

1,28

0,71

S7V3

22,90

1,00

0,64

S7V4

22,46

0,56

0,41

S7V5

22,23

0,33

0,24

In Figure 8, the behavior of ΔT can be observed in the cross section of the greenhouse for all the scenarios evaluated. The results show that for the CTG reference ventilation configuration S1 and for the most frequent speeds of the study region V1 to V3 which are equivalent to 83.5 % of the recorded wind speeds (Figure 1b), there is a marked climatic heterogeneity with average values of ΔT of 5.61 ± 3.04, 3.49 ± 1.98 and 2.73 ± 1.61 °C. This type of gradients affects the processes of transpiration and photosynthesis of the plants being able to generate a significant reduction in the final production of the crops as much in quality as in quantity [46, 47]. On the contrary, the most favourable scenarios were those simulated under speed V5 where the values of ΔT and its standard deviation S1 to S7 are below 1 °C, although it should be remembered that this speed is presented in a low percentage in the study area.

Figure 8. Thermal differential behavior for S1-S7

3.3.2 Wind direction SW-NE

The heat distribution fields for scenarios S8, S11 and S14 under external wind speed simulations V1, V3 and V5 are shown in Figure 9. How the temperature distribution depends on the flow field and ventilation rates [48], it can be observed that as these parameters did not present significant changes compared to scenarios S1 to S7, this translates into similar thermal distributions, where the colder zones agree with the fresh and less dense air intake zones of the environment and outside and the warmer zones with the zones of lesser air movement or air flow exit zones, which is consistent with similar studies such as the one carried out by Senhaji et al., [49].

The mean temperature values in CTG for the evaluated scenarios ranged between 27.66 and 22.27 °C for S8V1 and S14V5 respectively showing trends similar to those discussed in the previous numeral, the data obtained can be found in Table 3.

Figure 9. Simulated temperature distribution contours for S8-S14

Table 3. Thermal parameters obtained for simulated scenarios S8-S14

Scenarios

Tmean (°C)

ΔT (°C)

SD-ΔT (°C)

S8V1

27,66

5,76

3,09

S8V2

25,40

3,50

1,98

S8V3

24,60

2,70

1,59

S8V4

23,34

1,44

0,93

S8V5

22,73

0,83

0,54

S9V1

26,35

4,45

3,11

S9V2

24,93

3,03

1,88

S9V3

24,23

2,33

1,55

S9V4

23,20

1,30

0,95

S9V5

22,68

0,78

0,56

S10V1

24,48

2,58

1,62

S10V2

24,82

2,92

1,71

S10V3

24,10

2,20

1,44

S10V4

23,14

1,24

0,91

S10V5

22,65

0,75

0,54

S11V1

23,68

1,78

1,24

S11V2

24,09

2,19

1,27

S11V3

23,56

1,66

1,09

S11V4

22,86

0,96

0,71

S12V1

23,41

1,51

1,05

S12V2

23,76

1,86

1,08

S12V3

23,31

1,41

0,92

S12V4

22,72

0,82

0,60

S12V5

22,40

0,50

0,37

S13V1

22,95

1,05

0,84

S13V2

23,49

1,59

0,91

S13V3

23,10

1,20

0,78

S13V4

22,59

0,69

0,51

S13V5

22,32

0,42

0,30

S14V1

22,81

0,91

0,71

S14V2

23,32

1,42

0,80

S14V3

22,98

1,08

0,69

S14V4

22,52

0,62

0,45

S14V5

22,27

0,37

0,26

Figure 10. Thermal differential behavior for S8-S14

The distribution of the temperature variation for scenarios S8-S14 can be seen in Figure 10. The scenarios that exhibited a greater ΔT and greater heterogeneity were S8V1, S8V2 and S9V1 showing values of 5.73 ± 3.09, 3.50 ± 1.98 and 4.45 ± 3.11 °C respectively, while those of lower value were the high velocity V5 scenarios, with values of 0.37 ± 0.26, 0.42 ± 0.30 and 0.50 ± 0.37 °C for S14, S13 and S12 respectively, which ends up confirming the independence of the thermal behavior of CTG to the wind direction studied. It should also be pointed out that this greenhouse exceeds the recommended transverse length for passive greenhouses which, according to design criteria, must be less than 50 m in length [44], when analysing the graphs obtained from ΔT for S1-S7 and S8-S14, it can be deduced that if the greenhouse had been constructed respecting this design parameter, it would present a higher homogeneity in a large majority of the scenarios evaluated.

4. Conclusions

The optimization of air flow patterns in naturally ventilated greenhouses is a relevant factor of study in the search for improving the microclimatic conditions generated for plants, due to the direct relationship between microclimate and crop growth and development. In the present research, a numerical study was developed from a CFD model successfully validated with experimental data, which allowed to analyze the ventilation phenomenon in a traditional Colombian greenhouse. The current roof ventilation areas of the CTG greenhouse are inadequate and with the predominant wind speeds of the study region generates a heterogeneous and inadequate thermal behavior with average temperatures that can reach 27.51 ± 3.04 °C. Additionally, it was found that as the roof ventilation areas increase, the ventilation rates also increase within the range of 51 to 191 % with respect to the reference situation. This generates a positive effect on the thermal behavior of the CTG greenhouse obtaining average temperatures up to 16 % lower and a homogeneous profile with standard deviations of ± 1 °C. Likewise, the independence of the CTG greenhouse behavior with respect to the wind direction was verified.

Nowadays there are proposals for alternative designs of greenhouses adapted to the predominant climatic conditions of the Colombian cut flower production areas, such as those proposed by Villagran et al. [36, 42]. These proposals due to economic and logistical factors are likely to be welcomed in the medium or long term. The results obtained in this study showed that increasing the roof ventilation in the CTG greenhouse is a feasible technical alternative in the short term, which will improve the thermal behaviour of this type of greenhouse widely used in the Colombian ornamental sector. Future works that can complement the findings of this study should couple the increment of the roof ventilation areas with increasing the greenhouse height, consider the presence of crop in the simulations to get more realistic results and analyze the effect on other variables such as the CO2 exchange with the exterior.

Nomenclature

CTG

CFD

Cp

Colombian greenhouse, type chapel

Computational fluid dynamics

specific heat (J. kg-1. K-1)

$C_{1 \epsilon}, C_{2 \epsilon}, C_{\mu}$

coefficients in turbulence model

Dmi

observed temperature data (°C)

Dsi

Simulated temperature data (°C)

FMV

finite volume method

g

k

gravitational acceleration, (m.s-2)

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

k

turbulence kinetic energy (m2.s-2)

MAE

absolute mean error (°C)

RMSE

error of the mean quadratic root

S

scenari

SIMPLE

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

SD-ΔT

the standard deviations of ΔT

Sϕ

source term

Svc/Ssc

 Roof Ventilation Surface/Covered Floor Surface

T

Air temperature (°C)

Tmean

Average temperature (°C)

uj

components of speed (ms-1)

UDF

user defined function

V

Air speed (m.s-1)

WS

Wind Speed (m.s-1)

$Z_{0}$

the length of the roughness coefficient (m)

Greek symbols

Γϕ

the diffusion coefficient 

ΔT

thermal gradient (°C)

$\beta$

thermal expansion coefficient (K-1)

$\boldsymbol{\varepsilon}$

turbulent kinetic energy dissipation rate (m2.s-3)

µ

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

$\mu_{t}$

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

$\rho_{0}$

density (Kg.m-3)

$\sigma_{k} y \sigma_{\epsilon}$

Prandtl's turbulent numbers

$\Phi$

the ventilation rates as a function of the covered surface (m3m-2s-1)

  References

[1] Mena-Vásconez, P., Boelens, R., Vos, J. (2016). Food or flowers? Contested transformations of community food security and water use priorities under new legal and market regimes in Ecuador's highlands. Journal of Rural Studies, 44: 227-238. https://doi:10.1016/j.jrurstud.2016.02.011

[2] García-Cáceres, R.G., Felknor, S., Córdoba, J.E., Caballero, J.P.,  Barrero, L.H. (2012). Hand anthropometry of the Colombian floriculture workers of the Bogota plateau. International Journal of Industrial Ergonomics, 42(2): 183-198. https://doi:10.1016/j.ergon.2011.12.002

[3] Villagran Munar, E.A., Bojacá Aldana, C.R., Rojas Bahamon, N.A. (2018). Determinación del comportamiento térmico de un invernadero espacial colombiano mediante dinámica de fluidos computacional. Revista U.D.C.A. Actualidad & Divulgacion Cientifica, 21: 415-426. https://doi.org/10.31910/rudca.v21.n2.2018.1070

[4] Díaz, D.C., Bojacá, C.R., Schrevens, E. (2018). Modelando la idoneidad del invernadero de plástico tradicional para la producción de tomate en las regiones colombianas. Acta Horticulturae. 1205: 857–864. https://doi.org/10.17660 / actahortic.2018.1205.109 

[5] López, A., Valera, D.L., Molina-Aiz, F. (2011). Sonic anemometry to measure natural ventilation in greenhouses. Sensors, 11(10): 9820-9838. https://doi.org/10.3390/s111009820

[6] Bartzanas, T., Boulard, T., Kittas, C. (2004). Effect of vent arrangement on windward ventilation of a tunnel greenhouse. Biosystems Engineering, 88(4): 479-490. https://doi.org/10.1016/j.biosystemseng.2003.10.006

[7] Haxaire, R., Boulard, T., Mermier, M. (2000). Greenhouse natural ventilation by wind forces. In: Acta Horticulturae, 534: 31-40. https://doi.org/10.17660/ActaHortic.2000.534.2

[8] Critten, D.L., Bailey, B.J. (2002). A review of greenhouse engineering developments during the 1990s. Agricultural and Forest Meteorology, 112(1): 1-22. https://doi.org/10.1016/S0168-1923(02)00057-6

[9] Bournet, P., 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. https://doi.org/10.1016/j.compag.2010.08.007

[10] McCartney, L., Orsat, V., Lefsrud, M.G. (2018). An experimental study of the cooling performance and airflow patterns in a model Natural Ventilation Augmented Cooling (NVAC) greenhouse. Biosystems Engineering, 174: 173-189. https://doi.org/10.1016/j.biosystemseng.2018.07.005

[11] Espinoza, K., López, A., Valera, D.L., Molina-Aiz, F.D., Torres, J.A., Pena, A. (2017). Effects of ventilator configuration on the flow pattern of a naturally-ventilated three-span mediterranean greenhouse. Biosystems Engineering, 164: 13-30. https://doi.org/10.1016/j.biosystemseng.2017.10.001

[12] Baeza, E., Montero, J.I., Pérez-Parra, J., Bailey, B.J., Hernández, J.C., Gázquez, J.C. (2012). Avances en el estudio de la ventilación natural. http://www.publicacionescajamar.es/pdf/series-tematicas/centros-experimentales-las-palmerillas/avances-en-el-estudio-de-la-ventilacion.pdf, accessed on Feb. 27, 2018.

[13] Climate Change 2013 – The Physical Science Basis. Cambridge University Press, pp. 1-30. https://doi.org/10.1017/CBO9781107415324.004

[14] Mesmoudi, K., Meguellati, K., Bournet, P.E. (2017). Thermal analysis of greenhouses installed under semi arid climate. International Journal of Heat and Technology, 35: 474-486. https://doi.org/10.18280/ijht.350304

[15] Baeza, E.J. (2007). Optimización del diseño de los sistemas de ventilación en invernadero tipo parral. Doctoral Dissertation, Universidad de Almería. Spain.

[16] Villagrán, E.A., Gil, R., Acuña, J.F., Bojacá, C.R. (2012). Optimization of ventilation and its effect on the microclimate of a colombian multispan greenhouse. Agronomía Colombiana, 30(2): 282-288.

[17] Flores-Velázquez, J., Villarreal-Guerrero, F. (2015). Diseño de un sistema de ventilación forzada para un invernadero cenital usando CFD. Revista Mexicana de Ciencias Agrícolas, 6(2): 303-316. http://www.redalyc.org/comocitar.oa?id=263138086007, accessed on Oct. 7, 2017.

[18] Molina-Aiz, F.D., Valera, D.L., Peña, A.A., Gil, J.A., López, A. (2009). A study of natural ventilation in an Almería-type greenhouse with insect screens by means of tri-sonic anemometry. Biosystems Engineering, 104(2): 224-242. https://doi.org/10.1016/j.biosystemseng.2009.06.013

[19] Benni, S., Tassinari, P., Bonora, F., Barbaresi, A., Torreggiani, D. (2016). Efficacy of greenhouse natural ventilation: Environmental monitoring and CFD simulations of a study case. Energy and Buildings, 125: 276-286. https://doi.org/10.1016/J.ENBUILD.2016.05.014

[20] He, X., Wang, J., Guo, S., Zhang, J., Wei, B., Sun, J., Shu, S. (2018). Ventilation optimization of solar greenhouse with removable back walls based on CFD. Computers and Electronics in Agriculture, 149: 16-25. https://doi.org/10.1016/j.compag.2017.10.001

[21] 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(3): 847-854. https://doi.org/10.18280/ijht.360310

[22] Rico-Garcia, E., Soto-Zarazua, G.M., Alatorre-Jacome, O., De la Torre-Gea, G.A., Gomez-Melendez, D.J. (2011). Aerodynamic study of greenhouses using computational fluid dynamics. International Journal of. Physical Sciences, 6(28): 6541-6547. https://doi.org/10.5897/IJPS11.852

[23] Boulard, T., Fatnassi, H., Majdoubi, H., Bouirden, L. (2008). Airflow and microclimate patterns in a one-hectare canary type greenhouse: An experimental and CFD assisted study. Acta Horticulturae, 801(2): 837-845. http://doi.org/10.1016/j.agrformet.2009.01.002

[24] Launder, B.E., Spalding, D.B. (1972). Lectures in Mathematical Models of Turbulence. Acad Press, p. 176.

[25] Piscia, D., Muñoz, P., Panadès, C., Montero, J.I. (2015). A method of coupling CFD and energy balance simulations to study humidity control in unheated greenhouses. Computers and Electronics in Agriculture, 115: 129-141. https://doi.org/10.1016/j.compag.2015.05.005

[26] ANSYS-Fluent User’s Guide. (2013). ANSYS Fluent Theory Guide. Release, 182. https://doi:10.1016/0140-3664(87)90311-2

[27] Roy, J.C., Boulard, T., Kittas, C., Wang, S. (2002). Convective and ventilation transfers in greenhouses, Part 1: The greenhouse considered as a perfectly stirred tank. Biosystems Engineering, 83: 1-20. https://doi.org/10.1006/bioe.2002.0114

[28] Camacho, J.I., 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, 11(1): 32-46. https://doi.org/10.5424/sjar/2013111-411-11

[29] 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. https://doi.org/10.1016/j.biortech.2006.11.025

[30] Kim, K., Yoon, J.Y., Kwon, H.J., Han, J.H., Son, J.E., Nam, S.W., Lee, I.B. (2008). 3-D CFD analysis of relative humidity distribution in greenhouse with a fog cooling system and refrigerative dehumidifiers. Biosystems Engineering, 100(2): 245-255.. https://doi.org/10.1016/j.biosystemseng.2008.03.006

[31] Saberian, A., Sajadiye, S.M. (2019). The effect of dynamic solar heat load on the greenhouse microclimate using CFD simulation. Renewable Energy, 138: 722-737. https://doi.org/10.1016/j.renene.2019.01.108

[32] 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. https://doi.org/10.1016/j.jweia.2008.02.058

[33] Perén, J,I., van Hooff, T., Leite, B.C.C., Blocken, B. (2015). CFD simulation of wind-driven upward cross ventilation and its enhancement in long buildings: Impact of single-span versus double-span leeward sawtooth roof and opening ratio. Building and Environment, 96: 142-156. https://doi.org/10.1016/j.buildenv.2015.11.021

[34] Baeza, E.J., Sapounas, A., Stanghellini, C., Bonachela, S., Hernández, J., Montero, J.I., López, J.C., Granados, M.R., Muñoz, P., Lorenzo, P., Fernández-Del Olmo, P. (2017). Numerical simulation of the effect of different mulches on the heat storage capacity of a Mediterranean greenhouse soil. In: Acta Horticulturae, 1170: 119-127. https://doi.org/10.17660/ActaHortic.2017.1170.13

[35] Ansys. (2016). ANSYS ICEM CFD User Manual. https://doi.org/10.1016/j.joen.2015.02.033

[36] 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. https://doi.org/10.4081/jae.2019.856

[37] Flores-Velazquez, J. (2010). Analisis de la ventilación en los principales modelos de invernaderos en Mejico mediante dinámica de fluidos computacional (CFD). PhD Thesis. Almeria Spain, University Almeria.

[38] Archer, C.L., Jacobson, M.Z. (2003). Spatial and temporal distributions of U.S. winds and wind power at 80 m derived from measurements. Journal of geophysical research: Atmospheres, 108(D9). https://doi.org/10.1029/2002JD002076

[39] Bañuelos-Ruedas, F., Angeles-Camacho, C., Rios-Marcuello, S. (2010). Analysis and validation of the methodology used in the extrapolation of wind speed data at different heights. Renewable and Sustainable Energy Reviews, 14(8): 2383-2391. https://doi.org/10.1016/j.rser.2010.05.001

[40] Baxevanou, C., Fidaros, D., Bartzanas, T., Kittas, C. (2010). Numerical simulation of solar radiation, air flow and temperature distribution in a naturally ventilated tunnel greenhouse. Agricultural Engineering International: CIGR Journal, 12: 48-67. 

[41] Lee, S.Y., Lee, I.B., Kim, R.W. (2018). Evaluation of wind-driven natural ventilation of single-span greenhouses built on reclaimed coastal land. Biosystems Engineering, 171: 120-142. https://doi.org/10.1016/j.biosystemseng.2018.04.015

[42] Villagrán, E.A.,  Bojacá, C.R. (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: 133-143. https://doi.org/10.14295/oh.v25i2.2013

[43] Molina-Aiz, F.D., Valera, D.L.,  Peña, A.A., Gil, J.A. (2005). Optimisation of Almería-type greenhouse ventilation performance with computational fluid dynamics. In: Acta Horticulturae, 691: 433-440. https://doi.org/10.17660/ActaHortic.2005.691.52

[44] Baeza, E., Montero, J.I., Pérez-Parra, J., Bailey, B.J., Hernández, J.C., Gázquez, J.C. (2014). Avances en el estudio de la ventilación natural. Documentos Técnicos, (7).

[45] Kittas, C., Boulard, T., Papadakis, G. (1997). Natural ventilation of a greenhouse with ridge and side openings: sensitivity to temperature and wind effects. Transactions of the ASAE, 40: 415-425. https://doi.org/10.13031/2013.21268

[46] Majdoubi, H., Boulard, T., Fatnassi, H., Bouirden, L. (2009). Airflow and microclimate patterns in a one-hectare Canary type greenhouse: An experimental and CFD assisted study. Agricultural and Forest Meteorology, 149: 1050-1062. https://doi.org/10.1016/j.agrformet.2009.01.002

[47] Bojacá, C.R., Gil, R., Cooman, A. (2009). Use of geostatistical and crop growth modelling to assess the variability of greenhouse tomato yield caused by spatial temperature variations. Computers and Electronics in Agriculture, 65: 219-227. https://doi.org/10.1016/j.compag.2008.10.001

[48] Roy, J.C., Vidal, C., Fargues, J., Boulard, T. (2008). CFD based determination of temperature and humidity at leaf surface. Computers and Electronics in Agriculture, 61: 201-212. https://doi.org/10.1016/j.compag.2007.11.007

[49] Senhaji, A., Mouqallid, M., Majdoubi, H. (2019). CFD Assisted study of multi-chapels greenhouse vents openings effect on inside airflow circulation and microclimate patterns. Open Journal Fluid Dynamics, 09: 119-139. https://doi.org/10.4236/ojfd.2019.92009