Computational fluid dynamics analysis of heat transfer in a greenhouse solar dryer “chapel‐type” coupled to an air solar heating system

The distribution of turbulent kinetic energy (TKE), temperature, and velocity of humid air inside the greenhouse solar dryer (GHSD) was numerically investigated using 3D CFD ANSYS FLUENT code. The effect of solar radiation was coupled with the energy equations using the discrete ordinate model. Numerical simulations were based on two geometric models: Real model and model with reduced height, the solution was in good agreement with experimental data of temperature. The results of the real model showed that the TKE is ranged between 1.27 m2/s2 and 6 m2/s2 with an average of about 1.6 m2/s2 for the entire greenhouse dryer (GHD) volume. The greatest TKE magnitude is in the paths of the diffusers, which caused a temperature drop of about 2 K in the areas near the walls. Consequently, almost homogeneous temperature distribution was obtained in the entire volume of the GHD, although the average temperature was 315 K, and a gradient with respect to ambient temperature was of 14 K, that is, suitable for drying. Also, the average air velocity at 1 m height was 0.71 m/s, which is a value near the lowest limit (0.6 m/s) of forced convection drying. The improvement in the GHD by 36.5% volume reduction allowed an increase in the average TKE of 3.8 m2/s2 (2.4 times more than the previous one) located in the middle of the greenhouse; the average temperature reached 316.5 K with a gradient of 15.5 K, which represents an increase of 1.5 K (11%) compared to the real geometric model. The air velocity at 1 m height increased to 0.9 m/s in the improved geometric model (a growth of 35.7% compared with the previous geometry). More than 95% of the improved GHD volume has a uniform temperature, which is very suitable for a good quality drying process with higher speed.


| INTRODUCTION
The unremitting increase in the price of fossil fuels and the negative environmental effects of their use has promoted the use of renewable energies in several sectors of activity. 1 Among renewable energy sources, solar energy is the most promising alternative to fossil fuels because it is abundant on the planet, can be easily transformed into other forms of energy, and can be used for diverse purposes. In the industry, solar energy has been introduced as an essential resource in the drying for agricultural crop preservation, dairy processing, textile production, wood processing, wastewater treatment, etc. 2,3 The food manufacturing is one of the industries in which solar energy is increasingly used. In this industry, solar radiation can be directly used to dry food (direct mode), or it can be converted by solar collectors to thermal energy in solar thermal collectors to heat the air or water for drying processes (indirect mode). There is a third mode of use of solar radiation in drying called mixed-mode, which consists of the combination of the two previous modes. In any of the three modes, the drying process can be conducted by natural or forced convection. 4 Several technologies are used for the solar drying process. [4][5][6] The GHSD is more advantageous compared to other technologies; it has a simple construction with low operational costs; it has a high convective heat transfer (CHT) coefficient; it gives a very high quality of the dried product and reduces the drying time; 7,8 it can be used in mixed-mode with natural or forced convection; and under forced convection can reduce moisture content in lower time as well as increase drying rate compared with natural convection mode. 9,10 The main disadvantages of GHSD are that drying efficiency depends on climatic conditions and seasons; the need to use additional heat sources; the fact that facilities occupy a large area; and the possibility of odor emission released during drying of some products. 7 However, the design of the GHD must be based on sound scientific principles that facilitate good quality control of products during the drying process. 11 This design is specially based on the greenhouse shape in order to be more efficient in any season of the year. 12,13 The design of GHD must also include the characterization of air flow and the temperature distribution, since a bad distribution of temperature inside a GHD not only causes nonuniform production and quality, but also causes problems with pests and diseases. 14 Knowledge of temperature distribution should help to improve its homogeneity by modifying the design of the greenhouse structure and climate air conditioning systems and by locally controlling these systems. This requires the characterization and modeling of the processes, particularly of convective heat mass transfers, involved in its elaboration. 15 In the GHD design, the use of computational tools is very important, since they can allow the numerical testing of several shapes of the greenhouse geometries using different operating conditions, without the need to build the prototype of each geometric model. 16 The CFD is one of the most suitable tools for the design and analysis of GHSDs. 17,18 Its advantages over conventional experimental studies are the low cost and development time, the availability to study situations where experiments are not possible, and the ease of performing a large range of parametric studies for optimization. 19 CFD can help to manufacturers to improve greenhouse designs. It can produce quantitative predictions of air circulation and temperature distribution inside the GHDs under different climatic conditions, considering the external incident radiation distributed in three wavelength bands and the spectral optical properties of the involved materials. 14 It can also predict the mass transfer behavior in the food products during the drying process. [20][21][22][23] CFD enables easier studies of scalar and vector variables, which determines greenhouse microclimate with respect to its structural specifications and used equipment. 14,15,[24][25][26] There are several studies conducted on the greenhouse solar in which programming software such as MATLAB, C/C ++, Fortran, and TRNSYS has been used to simulate the behavior of air flow inside GHDs, 16,27-33 but very few studies have been done using the CFD software. Among the studies in which CFD was used, several of them were based on the twodimensional (2D) GHD geometric model, 14,17,18,24,[34][35][36] which cannot provide detailed information on the behavior of moist air in the entire GHD volume. Few studies have been carried out based on GHD's 3D geometric model. 26,[37][38][39][40][41] The use of the 3D CFD model allows obtaining a real image of the entire volume of the greenhouse, including at canopy level. Regarding the results obtained in the previous works, many studies have focused on the effect of the air flow variation on the temperature distribution inside the greenhouse. 17,34,38,39 The distribution of air velocity inside the greenhouse was not thoroughly investigated. 18,24,26,[35][36][37] The effect of the direct and diffuse solar irradiation on the air temperature and velocity distribution inside the GHSD, which can be easily simulated using CFD discrete ordinate (DO) model, 14 has not been taken into account in several previous studies. Another important aspect that is not wisely investigated is the air turbulent flow behavior in the GHD and its effect on the temperature distribution. As it knows, the CHT is intimately related to the TKE. [40][41][42][43][44] Increasing TKE in a fluid results in better CHT. 45 In the case of a GHD, the increase in the TKE of the air could lead to a homogeneity of air temperature inside the greenhouse. These behaviors can be easily analyzed using CFD simulation. The TKE is the mean parameter based on which the present study will be performed. The study consists of modeling and simulating a chapel-type GHSD using CFD code FLUENT. Two greenhouse geometric models are used in this study. The first geometric model is a representation of a prototype that serves as an experimental setup for the validation of numerical results. The second geometric model is the reduced form (63.5% of the volume) of the first one, which is used to study the effect of greenhouse volume reduction on the distribution of TKE, temperature, and velocity of humid air inside the GHD and also the relationship between TKE and temperature behavior. Both models are 3D and are simulated based on the same operating conditions. The study is focused on the behavior of air properties in the greenhouse and does not take into account the crop to be dried. The numerical simulation is performed in steady state combining energy, radiation DO, and viscous realizable k-e turbulent models in ANSYS FLUENT software. The thermophysical properties of moist air are calculated using EES software. 46

| Site and greenhouse description
The experimental setup is located in Temixco/Morelos (18° 51′N latitude, 99°14′O longitude and 1280 m altitude), in the southwest of Mexico. The experimental setup is composed of a chapel-type greenhouse ( Figure 1A,B) and a solar thermal system ( Figure 1D). As shown in Figure 1C, the solar thermal system (point 1) allows heating of the airflow before sending it through two diffusers to the greenhouse. The greenhouse has two inlets (point 2) and two outlets (point 3), both placed on the same path. At each outlet, a ventilator is installed to draw the air out of the greenhouse in order to maintain a constant air flow during the drying process. The dimensional characteristics of the greenhouse are described in Table 1. The solar system is composed of three flat plate solar collectors for air heating ( Figure 1D), which dimensions are length of 2.10 m, width of 1.20 m, and thickness of 0.093 m, and a fan installed in the outlet channel of the manifold (point 4) to involve air a constant mass flow rate of air inside the solar collectors.

| Measurement systems
The experimental data used to validate the numerical simulations were performed in the same GHSD. Measurements were focused on (i) the external solar irradiation; (ii) the air temperature and relative humidity inside and outside of the greenhouse; and (iii) the air velocity at the inlet and outlet of the greenhouse were performed in Temixco/Morelos in November and December 2017, in the absence of rain.

| Solar radiation measurements
The external solar irradiance (W/m 2 ) was measured with two pyranometers (Kipp & Zonen, CMP6, with an uncertainty of ± 4%). The hourly variations of the average direct and diffuse solar irradiance were measured during the 2 months of tests are presented in Figure 2.
The drying time considered for the present GHSD is 4 hours; due to this, an average solar irradiance between 11:00 and 15:00 hours was considered for this study. The corresponding average values of direct solar irradiance and diffuse solar irradiance were 761.05 W/m 2 and 62.21 W/m 2 , respectively.

| Inside air temperature measurements
The temperature of the air inside, at the inlet and outlet section of the greenhouse, was measured with type T thermocouples (uncertainty of ± 1 K). Sixteen measurement positions were chosen to place the temperature sensors as it is shown in Figure 3. Eight temperature sensors were placed at 1 m of the greenhouse height, while eight were at a height of 2 m. The coordinates of each measurement position are described in Table 2. Table 3 presents the average air temperature measurement from 11:00 to 15:00, for each position.
Based on the average temperatures (Table 3) and the air pressure (87588 Pa) inside the GHD, the thermophysical properties of moist air were calculated using the computational tool EES. 46 It was observed that the thermophysical properties of the moist air do not vary significantly with temperature for the range considered since the temperature range is from 309.2 K to 323.5 K. Therefore, the average values considered are reported in Table 4.

| External air temperature measurements
The external air temperatures on the GHD walls were measured with RTDs type PT1000 sensor (uncertainty of ± 0.2 K). Experimental data are shown in Table 5.

| Air velocity measurements
The measurement of inlet and outlet air velocities was recorded using an anemometer (Extech Instruments, SDL 350, with an uncertainty of ± 0.2 m/s). The average values obtained were 0.72 m/s and 0.59 m/s for West inlet (near to position 1) and East inlet (near to position 9), respectively, and 5.1 m/s and 4.4 m/s for West outlet (near to position 4) and East outlet (near to position 12), respectively ( Figure 3).

| Mathematical model
The air flow inside the greenhouse is assumed to be 3D, steady state, incompressible, and turbulent. 24 Flow and transport phenomena for air flow and the heat transfer are described by the Navier-Stokes equations for the momentum, mass, and energy are given as follow: Mass conservation equation:

Momentum conservation equation:
Energy conservation equation: The effect of turbulence on the air flow was implemented via the k-ε realizable model, which is a semi-empirical model based on transport equations for the TKE and turbulent dissipation rate (TDR). It is suitable for the analysis of turbulent flows in large spaces. 47 Turbulent kinetic energy:

Turbulent dissipation rate
In the equations (4) and (5), x j is the 3D coordinate systems, P k represents the generation of TKE due to the average velocity gradients; P b is the generation of TKE due to buoyancy; Y M is the fluctuating dilatation in compressible turbulence, and S k and S are root definitions. The coefficients are included in this model and can be defined as shown below 48 (1)

Where in the turbulent viscosity, we have
where Ω ij is the tensor of the center of rotation speed in a rotating zone with angular velocity k . They have constants A 0 and A s that describe by 48 where In order to simulate the effect of solar incident irradiance on the greenhouse cover, the DO model is used. In this model, it was assumed that irradiance energy was "convected" simultaneously in all directions through the medium at its own speed. The DO model allows the solution of irradiance at semitransparent walls. It can be used to nongray radiation using a gray-band model. So, it is adequate for use with participating media with a spectral absorption coefficient that varies in a stepwise fashion across spectral bands. The DO model solves the radiative transfer equation. The energy equations coupled with the radiation are presented below: where where ∆V is the control volume, S h i is the source term, and T ij is the coefficient due to convection and diffusion. S h i and T ij are evaluated based on the approximate results of the diffusion and convection, as well as the nonradiative terms of the source.

| Boundary conditions
Mixed boundary conditions including solar irradiance, natural convection of external air, and conduction through the semitransparent polyethylene were considered in the upper wall of the greenhouse. For the north, south, east, and west walls, constant temperatures were considered. The lower wall is considered adiabatic. All walls are stationary and noslip shear conditions. The diffuser boundary conditions are mass flow rate at the inlet and pressure at the outlet section. The values corresponding to each boundary condition including complementary parameters of irradiance, convection heat transfer, and turbulent flow are described in Table 6.

Sensor position (x, y, z) [m]
Sensor number

| Geometric model
The geometric models of the GHD under study were drawn in 3D using the ANSYS Workbench DesignModeler. The real geometric model ( Figure 4A) was drawn according to the dimensional characteristics of experimental chapeltype greenhouse, while the reduced height geometrical model ( Figure 4B) which represents 63.5% (61 m 3 ) of the real geometry volume (96 m 3 ) has the same dimensional characteristics as the first one, except the height that was reduced by 1 m. The use of the 3D geometric model allows obtaining a real image of the entire volume of the greenhouse.

| Mesh generation and sensibility analysis
The geometric models were meshed in ANSYS Workbench Meshing tools, using a structured hexahedral meshing mode ( Figure 5). The mesh independence of the results was verified by analyzing seven different mesh sizes including 1, 2, 3, 4, 5, 6, and 7 million of elements. Total heat transfer rate was chosen as the main parameter based on which this analysis was performed, since it reflects the overall behavior of moist air inside the GHSD. The results are shown in Figure 6, where the total heat transfer rate profile is no longer influenced by the mesh size when the number of mesh elements is larger than 6 million. Therefore, the mesh having 6 million of elements is chosen for this study in order to reduce the computational time. The characteristics of the selected mesh size are described in Table 7. It can also be observed that the quality control factors (orthogonality) of the selected mesh are within a favorable range. [48][49][50][51] The model with reduced height was meshed based on the same features as those of the real geometric mesh, in order to have the number of mesh elements greater than 6 million ( Figure 5B). This technique avoids a second mesh sensitivity analysis that would be an additional computational time-consuming. The quality control factors of reduced height geometric mesh are more suitable as shown in Table 7.

| Numerical simulation procedure
The numerical simulation is performed using the CFD software ANSYS FLUENT 19.0. Energy, radiation DO, and viscous realizable k-e turbulent models were used. The simulation is carried out in steady state, and the thermophysical properties of moist air are independent on temperature during simulation; they have been calculated based on the average air temperatures and pressure inside the greenhouse using the EES software. 46 The SIMPLE algorithm is used for pressure-velocity coupling, and the first order upwind scheme is used to discretize the convective terms in the momentum, energy, and turbulence equations. Regarding the initialization of the simulation, it is based on the average value of pressure, velocity (x, y, z), TKE, TDR, and temperature in both streams: 8.7 × 10 5 Pa, (0, 0, -0.58) m/s, 0.00164 m 2 / s 2 , 0.0003 m 2 /s 3 , and 323.37 K, respectively. The residual values are (1 × 10 -6 , 1 × 10 -5 , 1 × 10 -5 ) for x, y, z velocity, 1 × 10 −8 for energy, 1 × 10 −6 for epsilon (e), and 1 × 10 −6 DO-intensity. The computational time was about 36 hours with 8500 spatial iterations using 64-bit processors (Intel (R) Xeon (R) CPU E5-2620 v4 @ 2.10 GHz, cores: 16, installed ram: 32.0 GB).

| Validation of CFD results
To ensure the reliability of the results of this study, a comparison is made between the numerical results and the experimental data. Air temperature is the parameter based on which this comparison is carried out, since it was the parameter measured in detail in the GHD, and it is also the most important parameter in the drying process. The moist air temperatures from the nodes of the numerical model, which coincide with the different sensor positions, were plotted together with the experimental data in Figure 7, for four different paths previously described. Figure 7A,B shows the temperature results of positions at 1 m height and 1 m from the west wall, and 5 m from the west wall, respectively, while Figure 7C,D shows the temperature results of positions at 2 m height and 1 m from the west wall, and 5 m from the west wall, respectively. In general, the numerical results agree well with the experimental data. The moist air temperatures from the numerical model in most positions are within the range of the standard deviation tolerance (from 1 K to 2 K). The discrepancies in the numerical results beyond the standard deviation  Figure 7A,B, respectively. This error is due to measurement sensor uncertainty (± 1K).
Therefore, the numerical model of the GHD can predict adequately the behavior of experimental setup under study.
The mean deviation (% error), the reduced chi-square value (χ 2 ), and the root mean square error (RMSE) between the experimental data and the numerical results were evaluated using equations (1), (2), and (3), respectively.  The result of the error is 0.3315%, leading a minimal deviation between the experimental data and the results produced by the model. Considering a significance level (α) of 0.05 and the calculated χ 2 value of 1.1097, the difference in the experimental data compared with the numerical is not significant. Finally, the RMSE value is 1.020, indicating a good adjustment between the model and the experimental data.

| RESULTS AND DISCUSSION
The results of both models are presented graphically according to the different planes described in Figure 8. The orientation of the selected drawings is based on the different positions of the sensors as described in Figure 3. Figure 8A shows three lateral planes (PA, PB, and PC) at 1 m, 3 m, and 5 m with respect to the west wall, which are obtained through the cuts AA, BB, and CC, respectively. The Figure 8B also shows three vertical planes (PD, PE, and PF) at 0.5 m, 2.5 m, and 4.5 m front of the south wall, obtained by means of the cuts DD, EE, and FF, respectively.

| Air turbulence inside GHD
The forced CHT in a GHD is intimately related to the TKE. 44 Increasing TKE in a fluid could result in better CHT 45 and thus can improve the homogeneity of the air temperature inside the greenhouse. Figure 9 shows the contours of the TKE of moist air inside the greenhouse at different positions. Figure 9A-C shows the contour of the TKE distribution along the GHD at 1 m (PA), 3 m (PB), and 5 m (PC) from the west wall, respectively, whereas Figure 9D-F shows the contour in the transverse plane at 0.5 m (PD), 2.5 m (PE), and 4.5 m (PF) from the north wall, respectively. It is observed that the central area of the GHD is where the TKE is greater (6 m 2 /s 2 ), it is located in the middle of the path of the diffuser 1 ( Figure 9C,E) due to the increased air velocity inlet (0.72 m/s) compared to the diffuser 2 (0.59 m/s). The zone with greater turbulence is wide with an average of 3 m 2 /s 2 ( Figure 9C), which means that the distribution of the airflow is adequate to favor the homogeneity of temperature within the GHD. On the other hand, the lower turbulence area is located near the front and back walls as shown in Figure 9D,F. This is due to the decrease in the intensity of the air flow near the walls. The average magnitude of TKE in the entire volume is 1.6 m 2 /s 2 .

| Air temperature distribution in GHD
The distribution of the temperature of the humid air inside the GHD is described in the same way as the TKE distribution previously analyzed. It is shown in Figure 10A,C how the temperature of the hot air from the thermal solar collector decreases rapidly after entering the GHD. Under the effect of turbulence, associated with solar irradiance through the upper semitransparent wall, the temperature of moist air increases while becoming homogeneous in the large volume of the GHD. In the central zone ( Figure 10B,E), the moist air temperature is completely homogeneous. This is due to the increase in TKE, as previously described. The large central zone has a homogeneous temperature of about 317 K ( Figure 10B). The area near the south wall ( Figure 10A,C,D) has an elevated temperature due to the hot air flow coming from the solar air heaters. The temperature drops 3 K near the north wall ( Figure 10A,C) and 2 K in areas near the east and west walls ( Figure 10F). The average temperature inside the GHD is about 315 K; this temperature is suitable for food drying process since it is 14 K higher than the ambient temperature. 49

| Air temperature and velocity profile at 1 m GHD height
The distribution of humid air temperature and velocity at 1 m height inside GHD was also analyzed in this work. The choice of 1 m is justified because it is the height of the tables on which the products must be stored to be dried. The results are represented graphically in Figure 11. The profile of the average moist air temperature at 1 m height along the GHD is shown in Figure 11A. The temperature decreases from 320 K to 315 K between 0 and 1 m due to the high air temperature from solar thermal collectors. The rest of the zone has an almost constant temperature of 315 K, before undergoing a small drop of 1 K near the north wall. The gradient with respect to the minimum value of the drying process and environment is 11 K and 14 K, respectively. This means that the greenhouse under study is suitable for a better drying process. 52 The profile of the average humid air velocity at 1 m height along the greenhouse is shown in Figure 11B. The humid air velocity inside the GHD is ranged between 0.5 m/s and 3.75 m/s. The distribution is good between 0 m and 4.7 m, but the average magnitude is about 0.71 m/s, a value close to the lowest limit of air velocity (0.6 m/s) for the forced convection drying process. 53 This can lengthen the drying time, since not only the increase in temperature contribute to a good drying, the air must also have a slightly elevated velocity (from 0.6 m/s to 3 m/s) to guarantee a reduced drying time. 53 These results lead to propose an improvement in the structure of the GHSD under study, in order to increase the air velocity and temperature.

| Results of the model with reduced height
The results of the model with reduced height are presented and discussed in this subsection. Figure 12 shows the contours of the TKE in different planes previously described; it is observed that the area of major magnitude ( Figure 12A-C) is wider compared to the previous geometry. The greatest magnitude (11 m 2 / s 2 ) is in the middle of the GHD ( Figure 12E) compared to the previous geometric model where the greatest magnitude was in the trajectory of the diffusers. The average magnitude of TKE is 3.8 m 2 /s 2 , almost 2.4 times higher than that of the previous geometry due to a volume reduction of 36.5%. This behavior may favor a better homogeneity of the moist air temperature in the large part of the GHD volume, consequently, a faster drying. The contours of the moist air temperature are shown in Figure 13. Apart from the GHD entrance ( Figure 13A,C,D), the entire volume has a homogenous temperature of about 316.5 K, which represents an improvement of 1.5 K (11%) compared to the average temperature in the real model. This is due to the increase in the kinetic energy turbulent caused by the reduction of the GHD volume. The temperature drops 1 K near the north, east, and west walls ( Figure 13A,C,F). The model with reduced height reached 15.5 K higher than the ambient temperature, which is suitable for improving the drying time. 7,8 In Figure 14, the average moist air temperature and velocity profiles at 1 m height are compared between the two geometric models. Temperature behavior is shown in Figure 14A, where the temperature in the model with reduced height decreases from 322.5 K to 316.5 K between 0 and 1 m and remains almost constant (at 316.5 K) in rest of the zone. The gradient with respect to the minimum value of the drying process and the environment is 12.5 K and 15.5 K, respectively. Thus, an increase of 1.5 K is reached in the model with reduced height compared to the real model.
Similarly, in Figure 14B, velocity profiles in both geometric models are compared. The average air velocity in the model with reduced height reaches 0.9 m/s, which represents an increase of 0.2 m/s (28.2%) compared to the real geometric model average velocity (0.71 m/s). This is due to the reduction in GHD volume. Therefore, the reduction in GHD volume should contribute to reduce the drying time. 53

| CONCLUSIONS
The distribution of TKE, temperature, and velocity of humid air inside a GHD solar dryer was numerically investigated in this paper. The numerical simulations of two 3D GHD geometric models were carried out in ANSYS Fluent 19.0 using energy, radiation DO, and viscous realizable k-e turbulent models.
The numerical results were validated with the experimental data of temperature at the same GHD showing a good agreement. The different results obtained lead to the following conclusions: • The numerical model of the GHD can predict adequately the actual behavior under study. • The average magnitude of TKE in the entire volume of GHD real height is 1.6 m 2 /s 2 . The TKE is greater (6 m 2 /s 2 ) in the middle of the path of the diffuser 1. This behavior caused that the temperature distribution be not homogeneous in the whole volume of the GHD. • The temperature drop near the north wall was 3 K and 2 K in areas near the east and west walls. The average temperature inside the GHD was about 315 K, and the gradient with respect to the ambient temperature was 14 K, which is suitable for the food drying process, but the air velocity inside the GHD is low (0.71 m/s). This does not allow a very fast drying process; hence, it is necessary to improve it. • The reduction of the GHD volume of 35 m 3 (36.5%) leads an improvement in the average TKE of 3.8 m 2 /s 2 (2.4 times more than the previous one). • More than 95% of the GHD reduced height has a homogenous temperature of about 316.5 K. The gradient of 15.5 K with respect to ambient temperature was obtained, which represents an improvement of 1.5 K (11%) compared to the average temperature in the real model. • The air velocity average at 1 m height in the model with reduced height reached 0.9 m/s, which represents an increase of 0.2 m/s (28.2%) compared to the air velocity average of the real geometric model (0.71 m/s). • Thus, the reduction of the GHD volume under the same operating conditions can cause a very important change in the drying conditions for obtaining a minor drying time and higher drying velocity.
Based on the results of this study, it can be concluded that a CFD is an important tool for the design and optimization of GHD applied to the dehydration of agroindustrial products.

ACKNOWLEDGEMENTS
This study was partially supported by the CONACyT, Cátedras-CONACyT 352 project, and the postdoctoral UNAM-fellowship of Ituna-Yudonago Jean-Fulbert. The authors would like to thank Daniel Hernández Tamayo for his participation in experimental work.