Numerical Simulation of Multi-Span Greenhouse Structures

Greenhouses had to be designed to sustain permanent maintenance and crop loads as well as the site-specific climatic conditions, with wind being the most damaging. However, both the structure and foundation are regularly empirically calculated, which could lead to structural inadequacies or cost ineffectiveness. Thus, in this paper, the structural assessment of a multi-tunnel greenhouse was carried out. Firstly, wind loads were assessed through computational fluid dynamics (CFD). Then, the buckling failure mode when either the European Standard (EN) or the CFD wind loads were contemplated was assessed by a finite element method (FEM). Conversely to the EN 13031-1, CFD wind loads generated a suction in the 0–55° region of the first tunnel and a 60% reduction of the external pressure coefficients in the third tunnel was not detected. Moreover, the first-order buckling eigenvalues were reduced (32–57%), which resulted in the need for a different calculation method (i.e., elastoplastic analysis), and global buckling modes similar to local buckling shape were detected. Finally, the foundation was studied by the FEM and a matrix method based on the Wrinkler model. The stresses and deformations arising from the proposed matrix method were conservative compared to those obtained by the FEM.


Introduction
Although greenhouses extend all over Europe (Figure 1), the number of these structures are especially predominant in some regions due to the favorable local meteorological circumstances (e.g., Mediterranean climate characterized by warm and humid winters, hot and dry summers and clear days with fairly high radiation even during fall and winter) in which crop production in greenhouses is a rising agricultural sector.
Among European nations (Figure 1), Spain is the country with the highest greenhouse area and accounts for the 35% of the total reported statistics on protected cultivation for the latest available dataset (i.e., 2016) provided by Eurostat [1]. Nonetheless, Italy, France, the Netherlands, Poland and Greece also display greenhouse area values greater than the European average (4729 ha). Among European nations (Figure 1), Spain is the country with the highest greenhouse area and accounts for the 35% of the total reported statistics on protected cultivation for the latest available dataset (i.e., 2016) provided by Eurostat [1]. Nonetheless, Italy, France, the Netherlands, Poland and Greece also display greenhouse area values greater than the European average (4729 ha).
Due to the relative importance in Europe, Figure 2 illustrates the evolution of the crop production area under greenhouse in Spain between 2002 and 2019. Although not homogenously, protected horticulture increased up to 28.6% in the last 18 years.     Among European nations (Figure 1), Spain is the country with the highest greenhouse area and accounts for the 35% of the total reported statistics on protected cultivation for the latest available dataset (i.e., 2016) provided by Eurostat [1]. Nonetheless, Italy, France, the Netherlands, Poland and Greece also display greenhouse area values greater than the European average (4729 ha).
Due to the relative importance in Europe, Figure 2 illustrates the evolution of the crop production area under greenhouse in Spain between 2002 and 2019. Although not homogenously, protected horticulture increased up to 28.6% in the last 18 years. Last available estimates show that 70,744 ha of Spanish crop production occurred in greenhouses [2]. Nonetheless, a great disparity in figures exists among the different Autonomous Communities (Figure 3a). Thus, it is worth mentioning the relevance of Andalusia with 53,902 ha, which represents 76.2% of the country total. Moreover, the Almeria province is the largest contributor with 32,958 ha, i.e., 61.1% of the total Andalusian territory (i.e., 46.56% of the national area) (Figure 3b). The growth of greenhouse area in Almeria, commonly known as the "sea of plastic", has been unprecedented from 1970 onwards, when only 3000 ha was reported [20]. 55017   Greenhouses are structures used for the cultivation and/or protection of plants and crops under controlled conditions (i.e., light, temperature, humidity and air composition) to create a favorable growing environment for a wide range of products that could increase the produced yield and also extend the production out of season.
In order to present a better adaptation to the site-specific climatic conditions, different typologies of greenhouses have been designed. An exhaustive analysis on the types of greenhouses could be found in von Elsner et al. [21,22].
In Spain, greenhouses are predominantly lightweight and low-cost structures [23]. Historically, the most popular Spanish greenhouse was the flat parral-type, which originated from Almeria and is characterized by an almost non-existent slope. However, due to the increase in the level of automation, the multi-span or multi-tunnel greenhouse (as could be observed in Figure 4) is the most widely greenhouse structure employed to date [24]. Although greenhouses are subjects of interest in the technical and scientific community, the vast Greenhouses are structures used for the cultivation and/or protection of plants and crops under controlled conditions (i.e., light, temperature, humidity and air composition) to create a favorable growing environment for a wide range of products that could increase the produced yield and also extend the production out of season.
In order to present a better adaptation to the site-specific climatic conditions, different typologies of greenhouses have been designed. An exhaustive analysis on the types of greenhouses could be found in von Elsner et al. [21,22].
In Spain, greenhouses are predominantly lightweight and low-cost structures [23]. Historically, the most popular Spanish greenhouse was the flat parral-type, which originated from Almeria and is characterized by an almost non-existent slope. However, due to the increase in the level of automation, the multi-span or multi-tunnel greenhouse (as could be observed in Figure 4) is the most widely greenhouse structure employed to date [24]. Greenhouses are structures used for the cultivation and/or protection of plants and crops under controlled conditions (i.e., light, temperature, humidity and air composition) to create a favorable growing environment for a wide range of products that could increase the produced yield and also extend the production out of season.
In order to present a better adaptation to the site-specific climatic conditions, different typologies of greenhouses have been designed. An exhaustive analysis on the types of greenhouses could be found in von Elsner et al. [21,22].
In Spain, greenhouses are predominantly lightweight and low-cost structures [23]. Historically, the most popular Spanish greenhouse was the flat parral-type, which originated from Almeria and is characterized by an almost non-existent slope. However, due to the increase in the level of automation, the multi-span or multi-tunnel greenhouse (as could be observed in Figure 4) is the most widely greenhouse structure employed to date [24]. Multi-span or multi-tunnel greenhouse on an irregularly-shaped plot [24].
Although greenhouses are subjects of interest in the technical and scientific community, the vast majority of papers in this regard are focused on the management of the controlled environment (e.g., heating systems [25][26][27][28] and natural/forced ventilation [29][30][31][32]). On the contrary, to date, the number of research works related to the structural calculation of greenhouses is more limited. There are some  Although greenhouses are subjects of interest in the technical and scientific community, the vast majority of papers in this regard are focused on the management of the controlled environment (e.g., heating systems [25][26][27][28] and natural/forced ventilation [29][30][31][32]). On the contrary, to date, the number of research works related to the structural calculation of greenhouses is more limited. There are some studies based on the structural analysis of greenhouses [22,24,[33][34][35][36][37][38][39][40][41][42][43][44][45][46][47][48][49], whereas, to the best of the authors' knowledge, only one paper focused on the analysis of foundations in greenhouses [50].
Furthermore, the aforementioned researchers have shown their concerns about the widespread greenhouse structure construction practice based on pre-engineered solutions without involving a specific technical project [22,38,42,45,46,51,52]. Similarly, Ha et al. [39] also reported that greenhouses foundations are regularly empirically calculated regardless of the specific soil mechanics. This approach could either lead to structural inadequacies or cost ineffectiveness.
In Europe, the EN 13031-1 [53] standard establishes the requirements for the design and construction of commercial production greenhouses, which are designed to sustain the characteristic permanent and concentrated vertical loads, as well as the structural actions imposed by the distinct crop production and the site-specific climatic conditions. Among the literature, wind loads are reported as one of the principal sources for damage and failure of lightweight greenhouses structures [34,39,44,[54][55][56]. In fact, as lightweight structures located in the Mediterranean Basin, both the dead and snow load are low, whereas wind load is the predominant design factor. Thus, in-depth knowledge of the wind effect on the structure and foundation of greenhouses is essential to ensure the overall structural safety. In this regard, wind action on structures of Eurocode 1-4 [57] should also be considered. However, as the recommendations of the external wind pressure coefficients on curved roofs are different to those conveyed by the EN 13031-1 [53] standard, an analysis is required to assess the implications of such differences on the greenhouse design.
Therefore, this research was firstly focus on the study of the discrepancies between both aforementioned standards and the simulation of the external wind pressure coefficients by means of Computational Fluid Dynamics (CFD). And then, the paper addressed the influence of the wind load consideration on the structural behavior, as well as the foundation through the development of a three-dimensional (3D) matrix calculation method based on the Winkler model to design the cylindrical footings.

Geometry
The greenhouse dimensions are greatly influenced by the crop type and the desired yield as well as the constraints imposed by the plot of land. For instance, Spanish multi-span greenhouses, specifically those in the south-east of the country, range between 0.06 m 2 and 3.25 ha, with an average area of 1.28 ha [58].
The vast majority of the multi-span greenhouses in the Mediterranean basin comprise two to 15 modular structures [38]. Nonetheless, Carreño-Ortega et al. [36] and Vázquez-Cabrera et al. [48] reported that greenhouses that comprise three tunnels are subjected to greater stresses than those that comprise a higher number of modules. The authors also noted that the stress decreased for increasing numbers of tunnels and stated that hall lengths above 50 m marginally affected the stress outcomes. Therefore, the structural assessment of multi-span greenhouses could be limited to the analysis of three tunnels.
Usually, each module has a 6-9.3 m width, 3-4.5 m eave height and a 4.35-6.2 m ridge height [24,41]. It is worth mentioning that the greenhouse structure also presents a unit of repetition through the length of the structure. The central portal frames of the greenhouses commonly consist of two different constructions (i.e., primary and secondary frame), which are separated 2.5 m, alternating between each other until the final length of the greenhouse is reached ( Figure 5). The primary portal frame, which serves to sustain the crop, comprises trussed arches supported by central and lateral columns, whereas the secondary portal frame consists of simple arches supported by the lateral columns and is connected to the rest of the structure through the gutters and purlins.
Despite a greenhouse structure requires a 3D definition to include the influence of the external loads as well as the effect of the structural elements responsible for the connections (i.e., purlins and gutters), Briassoulis et al. [33] demonstrated that the study of the main planes (two-dimensional (2D) analysis) was a valid approach in the assessment of greenhouses. In fact, such simplification has been commonly employed in both industrial and agricultural constructions [59][60][61][62][63]. In the same way, Ren et al. [47] also performed a similar simplification by means of 3D models in order to decrease the computational cost.
Thus, in this research work, a multi-tunnel greenhouse located in Almeria, due to the relative importance of this type of structure as well as this region regarding greenhouse density area in Europe, was employed as the basis of the calculations. The structure (namely, type B10 according to EN 13031-1 [53]) was a multi-span metal-frame covered with a single-sheet of light-translucent plastic material (such as low density polyethylene), comprised of three tunnels and two differentiated portal frame types ( Figure 6). Despite a greenhouse structure requires a 3D definition to include the influence of the external loads as well as the effect of the structural elements responsible for the connections (i.e., purlins and gutters), Briassoulis et al. [33] demonstrated that the study of the main planes (two-dimensional (2D) analysis) was a valid approach in the assessment of greenhouses. In fact, such simplification has been commonly employed in both industrial and agricultural constructions [59][60][61][62][63]. In the same way, Ren et al. [47] also performed a similar simplification by means of 3D models in order to decrease the computational cost.
Thus, in this research work, a multi-tunnel greenhouse located in Almeria, due to the relative importance of this type of structure as well as this region regarding greenhouse density area in Europe, was employed as the basis of the calculations. The structure (namely, type B10 according to EN 13031-1 [53]) was a multi-span metal-frame covered with a single-sheet of light-translucent plastic material (such as low density polyethylene), comprised of three tunnels and two differentiated portal frame types ( Figure 6).
In order to analyze an average commercial multi-tunnel greenhouse, the selected dimensions were within the commonly encountered in the zone (i.e., 8 m width, 4.5 m eave height and 6 m ridge height). Since the modules were contained in two consecutive planes of symmetry, the 2D calculations were made possible without an excessive computational cost. The primary portal frame consisted of a truss consisting of vertical bars separated 2 m on the bottom chord, which was connected to the central purlin of each module through a 45 • stay. Moreover, the lateral columns of both primary and secondary portal frames were connected by purlins located at 3/8 and 6/8 from the base of the eave height. No limitations in the depth of the greenhouse were contemplated, as the external loads applied to the structure (permanent, concentrated vertical, crop, snow and wind actions according to EN 13031-1 [53], see Section 3.1) are constant for all central portal frames (i.e., except for the front and back-end frames) regardless of the total length. importance of this type of structure as well as this region regarding greenhouse density area in Europe, was employed as the basis of the calculations. The structure (namely, type B10 according to EN 13031-1 [53]) was a multi-span metal-frame covered with a single-sheet of light-translucent plastic material (such as low density polyethylene), comprised of three tunnels and two differentiated portal frame types ( Figure 6).   Figure 7 displays the different types of connections established between the structural elements of the greenhouse as considered in the Ansys model [64]. In this regard, a rigid connection was present in the joint between the columns and the arches as well as the columns with the terrain, whereas most of the remainder nodes were articulated ( Figure 8). Such articulated connections were employed to connect a bar to a plate or to other bars by means of bolts or clamp fasteners.

Boundary Conditions
Agriculture 2020, 10, x FOR PEER REVIEW 6 of 32 In order to analyze an average commercial multi-tunnel greenhouse, the selected dimensions were within the commonly encountered in the zone (i.e., 8 m width, 4.5 m eave height and 6 m ridge height). Since the modules were contained in two consecutive planes of symmetry, the 2D calculations were made possible without an excessive computational cost. The primary portal frame consisted of a truss consisting of vertical bars separated 2 m on the bottom chord, which was connected to the central purlin of each module through a 45° stay. Moreover, the lateral columns of both primary and secondary portal frames were connected by purlins located at 3/8 and 6/8 from the base of the eave height. No limitations in the depth of the greenhouse were contemplated, as the external loads applied to the structure (permanent, concentrated vertical, crop, snow and wind actions according to EN 13031-1 [53], see Section 3.1) are constant for all central portal frames (i.e., except for the front and back-end frames) regardless of the total length. Figure 7 displays the different types of connections established between the structural elements of the greenhouse as considered in the Ansys model [64]. In this regard, a rigid connection was present in the joint between the columns and the arches as well as the columns with the terrain, whereas most of the remainder nodes were articulated ( Figure 8). Such articulated connections were employed to connect a bar to a plate or to other bars by means of bolts or clamp fasteners.

Structural Analysis
Ansys [64] software was used to apply the finite element method required for the numerical simulation of the greenhouse structure. The model was based on steel beam elements (Beam 4) with an elasticity module of 210 GPa and a Poisson coefficient of 0.3. Beam4 is an uni-axial element with compression, tension, torsion and bending capabilities, which present six degrees of freedom at each node. In this regard, the type of beam element is similar to those employed in Roux and Motro [65] and Roux et al. [66,67].
A sensitivity analysis showed that the use of eight elements to define a bar resulted in a correct estimation of both the stresses and the buckling modes as well as adequately discretizing the arch as a polygonal. In this regard, the selection of eight elements is in line with the five-element limit established in Simões et al. [68] to correctly estimate the critical loads. Due to the relatively low computational cost, the number of elements was also maintained for all the bars in the structure (columns, purlins, stays, etc.). Figures in Section 3.2 show the geometry of the greenhouse before and after buckling, which in their non-deformed state serve as an evidence of the mesh employed according to the number of elements considered. Regarding the cross-section of the structural elements, the columns consisted of square hollow Section 80/3 (i.e., 80 mm in external dimensions and 3 mm in thickness), the arches were designed with a 60 mm diameter and 2 mm thickness tubular section and the bottom chord and the vertical members of the truss and the stay were made up of 40 mm diameter and 2 mm thickness circular section.
The research carried out by Roux et al. [66] set the basis for the buckling study of greenhouse arches as it was included as a regulatory annex for the stability check of greenhouses arches in EN 13031-1 [53]. The evaluation of either the first-order buckling eigenvalue (λ cr ), the second-order elastic critical load factor (α cr ) or the second-order elastoplastic critical load factor (α u ) was proposed as a reference in the estimation of the collapse load of greenhouses. In this regard, information concerning the λ cr /α u ratio could be found across the literature. For the greenhouses evaluated by Roux et al. [66], λ cr / α u ratios close to 1.7 were found. Similarly, for a single-span round arched greenhouse, Ren et al. [47] observed a lower λ cr /α cr ratio, around 1.5.
In this paper, to assess the influence of the external actions on the greenhouse, the buckling analysis was carried out through the study of λ cr , which was calculated by means of the Subspace Method though a subroutine created in the Ansys Parametric Design Language (APDL) language. The study of the λ cr was preferred since the determination of α u is more complex and multi-span greenhouses display a quite stable λ cr /α u ratios. Moreover, as the λ cr was determined for the different load combinations on the same structure, the λ cr was used as a strength index and allowed the investigation of the most unfavorable conditions. Furthermore, as λ cr is a parameter commonly employed by the different regulatory standards on steel structures, a general extrapolation and comparison approach is attainable between different authors, regardless of the standard employed in their study or the calculation method (i.e., elastic or plastic).

According to the Current Standards
In this research work, the European standards (EN 13031-1 [53], Eurocode 1-4 [57]) and Spanish standards (UNE 76209 IN [69] and UNE 76210 IN [70]) were employed to analyze the different external loads contemplated:

1.
Permanent loads: Those produced by the self-weight of the structural metal frame and the film cladding system (0.092 kg·m −2 for a 10 −4 m film thickness).

2.
Concentrated vertical action: In order to account for the weight of the worker in the maintenance labor, it is considered as a normal force of 0.35 kN. The force is applied to the middle of the secondary portal frame, which is the most unfavorable situation.

3.
Crop load: Usually, the greenhouse structure supports the weight of the plants and the crop products. The standard establishes this load as a uniformly distributed vertical force applied to the bottom chord of the trussed roof frame. A typical value for tomatoes and cucumbers of 0.15 kN·m −2 was considered in this research work.

4.
Snow load: According to the Spanish standard UNE 76210 IN [70], the corresponding snow load for Almeria, which is located at an altitude of less than 200 m, should be contemplated as an accidental load (A k ) with a value of 0.19584 kN·m −2 , as determined by Equation (1): where µ i is the snow load shape coefficient (µ i = 0.8 for a curve multi-span greenhouse as defined in the EN standard), c t is the thermal coefficient (c t = 1 due to the lack of heating system as defined in the EN standard) and c n is the correction coefficient for a specified return period (c n = 0.612 for a 10-year lifespan as defined in the EN standard).

5.
Wind loads: According to EN 13031-1 [53], Eurocode 1-4 [57] and UNE 76209 IN [69], greenhouses in which windows could be open and closed must be designed as a structure without openings. In fact, the lack of interior wind loads is the most realistic option for greenhouses that are professionally constructed, as the air leakage, which causes loss of heat and carbon dioxide, is completely controlled.
This action is dependent on the location of the structure that, according to the Eurocode 1-4 [57], is considered terrain category II, i.e., areas with low vegetation such as grass and isolated obstacles (trees or buildings) with separations of at least 20 obstacle heights. The basic probable wind velocity (v b ) reached 23.563 m·s −1 and was defined as a function of wind direction and time of year at 10 m above ground of terrain category II was determined by Equation (2): where c dir is the directional factor (c dir = 1 as defined in the EN standard), c season is the seasonal factor (c season = 1 as defined in the EN standard), c prob is the probability factor (c prob = 0.873 for a 10-years return period, which corresponds to a type B10 greenhouse, as defined in the EN standard, and a K factor of 0.33 according to the Spanish national guidelines [69]) and v b,0 is the fundamental value of the basic wind velocity, which is established as 27 m/s for Almería in UNE 76209 IN [69]. Moreover, the mean wind velocity (vm(z)) of 20.74 m·s −1 was determined by Equation (3) at a reference height (z) of 5.25 m (i.e., the mean of the eave and ridge height values and greater than the where c r(z) is the roughness factor (c r(z) = 0.88 as defined in the EN standard) and c 0(z) is the orography factor (c 0(z) = 1 as defined in the EN standard). Then, the peak velocity pressure (q p(z) was calculated to account for both the mean and short-term velocity fluctuations according to Equation (4): where I v(z) is the turbulence intensity (I v(z) = 0.215 as defined in the EN standard) and ρ is the air density (ρ = 1.25 kg·m −3 ). Thus, at the reference height of 5.25 m, the wind load reached a peak velocity pressure of 0.679 kN·m −2 .
Finally, in order to calculate the wind pressure in the different zones on the surface of the greenhouse, the peak velocity pressure was multiplied by the external pressure coefficient corresponding to the different zones. Table 1 shows the different values of the external pressure coefficients (cp e ) conforming to EN 13031-1 [53], which depend on the wind direction and the area of the greenhouse subjected to the wind action (i.e., number of tunnel and roof angle) and the related wind load values. Furthermore, Figures 9 and 10 display the wind loads for the different zones of the greenhouse structure.

Wind Load Action Simulated Through CFD
The findings reported by Kim et al. [44], who observed a good correlation between the results arising from a CFD numerical simulation and a wind tunnel experiment, are used as basis for the CFD simulation carried out in this paper. The authors in [43] analyzed different greenhouse geometries, mesh and computational domain sizes and turbulence models.
Regarding the domain size, accurate and efficient computation was achieved for values of 15 times the height of the greenhouse in the downstream zone and three times the height of the greenhouse in the windward zone [43,44]. The reported domain size was consistent with the values used in other numerical models [71,72].
In the literature [29,[73][74][75][76], the k-ε turbulence model is commonly employed in CFD simulations of greenhouses with or without insect-proof screens. Nonetheless, it is worth mentioning that these simulations are focused on a crop level, i.e., an interior zone of the greenhouse defined by a porous jump that is the sink or source of humidity as a transport medium. In turn, when the CFD capabilities are used to define the wind loads to be considered in the design of greenhouses, the recommendations regarding turbulence models put forward by Kim et al. [43,44] are more suitable. The authors indicated that the k-ω shear stress transport (k-ω SST) turbulence model was the best representation for the experimental results within the Reynolds-averaged Navier-Stokes (RANS) models. The same conclusion was also reported by Xing et al. [77] when the k-ω SST was compared to the large eddy simulation (LES) models.
In this paper, the use of a 2D numerical simulation with CFD is proposed as an alternative to the calculation of the 0 • -wind load determined by the EN 13031-1 [53].
Firstly, the numerical simulations carried out in the CFD model proposed in this research work had to be validated. Then, Ansys Fluent software [78] was employed to compare the new CFD-modeled results and the wind tunnel experimental results reported by Kim et al. [44].
Therefore, a previously CFD-tested greenhouse structure (i.e., the so-called even-span by Kim et al. [44], hereinafter referred as EV-D22, defined by a 7 m width, 2 m eave height, 3.41 m ridge height and 22 • roof slope) was initially implemented to the model. Moreover, in the velocity inlet, the same turbulence, energy dissipation and wind logarithmic profile were employed to carry out the validation. Similarly, a k-ω SST turbulence model, a first cell height of 1.5·10 −4 m and a y+ < 1 were also adopted. A pressure-based solver using the simple-algorithm and identical viscosity and air density values were assumed as well. Figure 11a illustrates the mesh size employed in the CFD model applied to the gabled-roof greenhouse defined in Kim et al. [44].  Then, in order to simulate the transversal wind load applied to the specific greenhouse studied, the validated CFD model was adapted to the structural dimensions proposed in this research work, which resulted in the mesh size distribution shown in Figure 11b, as well as the wind profile and turbulence conditions ( Figure 12) for the terrain category II established by the Eurocode 1-4 [57]. Then, in order to simulate the transversal wind load applied to the specific greenhouse studied, the validated CFD model was adapted to the structural dimensions proposed in this research work, which resulted in the mesh size distribution shown in Figure 11b, as well as the wind profile and turbulence conditions ( Figure 12) for the terrain category II established by the Eurocode 1-4 [57].
Then, in order to simulate the transversal wind load applied to the specific greenhouse studied, the validated CFD model was adapted to the structural dimensions proposed in this research work, which resulted in the mesh size distribution shown in Figure 11b, as well as the wind profile and turbulence conditions (Figure 12) for the terrain category II established by the Eurocode 1-4 [57]. The new CFD model domain is based upon the same principles as the one proposed by Kim et al. [43,44]. Thus, the same relative distance between the velocity inlet and the greenhouse, the pressure outlet, the upper zone and y+ were applied. Nonetheless, the implementation of the specific geometry and dimensions of the multi-span greenhouse made necessary the redefinition of the boundary conditions regarding to the turbulence velocity and dissipation in the velocity inlet.

Foundation
The foundations of greenhouse structures are generally composed of cylindrical concrete footings of approximately 0.5 m in diameter and a depth between 1 or 2 times the diameter. Thus, once the diameter is fixed, the objective is to assess the optimum depth in order to decrease the overall costs.
It is worth mentioning that greenhouse foundations display characteristics between those of shallow and deep foundations, since part of the stresses received from the column are transmitted to the ground both through the lateral walls and the base of the footing. Usually, the distribution of the The new CFD model domain is based upon the same principles as the one proposed by Kim et al. [43,44]. Thus, the same relative distance between the velocity inlet and the greenhouse, the pressure outlet, the upper zone and y+ were applied. Nonetheless, the implementation of the specific geometry and dimensions of the multi-span greenhouse made necessary the redefinition of the boundary conditions regarding to the turbulence velocity and dissipation in the velocity inlet.

Foundation
The foundations of greenhouse structures are generally composed of cylindrical concrete footings of approximately 0.5 m in diameter and a depth between 1 or 2 times the diameter. Thus, once the diameter is fixed, the objective is to assess the optimum depth in order to decrease the overall costs.
It is worth mentioning that greenhouse foundations display characteristics between those of shallow and deep foundations, since part of the stresses received from the column are transmitted to the ground both through the lateral walls and the base of the footing. Usually, the distribution of the bending moment between the walls and the base of the footing is hardly detectable by analytical methods. Thus, a matrix method is a more suitable procedure to assess the stress distribution.
A new matrix model (MM) of 3D articulated bars programmed with Matlab [79] and based on the Winkler model was developed to study the foundation of greenhouse structures. Some research works focused on the Wrinkler model applied to caisson foundations can be found in the literature [80][81][82]. These works, which include lateral, rotational springs and dashpots models, have been successfully employed in the calculation of bridge foundations under dynamic conditions and layered soils. Nevertheless, they rely on a complexity of calculation that escapes from the technology frequently used for the calculation of greenhouse structures. Thus, in this paper, the proposed methodology is easier and compatible with the typical commercially available software based on a matrix analysis with 3D elements. Furthermore, the use of this simple methodology favors a seamless implementation in optimization tools such as genetic algorithms [83] while avoiding the convergence problems that other non-linear calculation procedures would generate.
The MM was defined by typical 3D truss elements that, in local coordinates, only present an axil force (either a compression or tensile force) and a degree of freedom (i.e., displacement on the directrix of the element). However, when global coordinates are considered, the displacement on the bar directrix decomposes into displacements along the global axes. Although 3D truss elements are articulated, a rigid solid is created when several of these elements are joined together, which allows the interaction with the Wrinkler's springs.
Moreover, the proposed MM model was compared to a solid-elastic contact model through the finite element method (FEM) programmed with Ansys APDL [64]. The assessment between models was performed on a specified cylindrical concrete footing of 0.5 m diameter and 0.8 m depth.
In the subsequent sections, the methodologies employed for both the matrix model and the solid-elastic model are detailed. Once the geometry of the network of bars was characterized, the relationship between the solid and the ground was defined through articulated springs based on the Winkler model. In this approach, the soil is represented by a number of elastic springs whose rigidity depends on the subgrade reaction modulus and their area hinges on the range of influence around the perimeter of the foundation.
To avoid a false redistribution of stresses caused by springs subject to a tensile load, which would avoid the separation between the footing and the ground, the calculation was programmed iteratively by modifying the modulus of elasticity of each aforementioned spring to an almost zero value. A Matlab [79] subroutine was created to perform a matrix calculation of the footings and to identify the elements representing a terrain subjected to tensile stresses. Then, an almost zero crosssection area was assigned to those elements, which effectively deleted the previously considered spring, and a new matrix calculation was carried out. Once again, the recalculated elements subjected to tensile stresses were singled-out, their area-value was changed to virtually zero and a subsequent matrix calculation was implemented. This loop was repeated until the number of compressed springs remains stable after two consecutive iterations (i.e., the end of the calculation is reached when no variations in the number of springs subjected to tensile loads is detected after two iteration cycles).
From the results arising from the springs, it was possible to ascertain the pressure generated by the concrete footing on the ground, as well as the specific stress redistribution generated by the lateral walls and base of the element, which ultimately would allow for the calculation of the greenhouse In order to achieve an overall behavior as a rigid solid body, each subsector is composed of articulated bars arranged to form a cross on each plane to ensure a total rigidity as well as the transmission of forces. Whereas Figure 13a only illustrates the external bars in the subsectors, the real number of bars conforming two 45 • horizontal subsectors is displayed in the close-up detail (Figure 13b). Note that the intermediate bars cross all the existing rectangles contained in all the horizontal planes of the structure. In this manner, the force transmission from each node to all directions was guaranteed. The network of articulated bars was defined with an area of 1 m 2 and a modulus of elasticity of 1000 MPa to ensure the non-deformable and rigid-block behavior of the assembly.
Once the geometry of the network of bars was characterized, the relationship between the solid and the ground was defined through articulated springs based on the Winkler model. In this approach, the soil is represented by a number of elastic springs whose rigidity depends on the subgrade reaction modulus and their area hinges on the range of influence around the perimeter of the foundation.
To avoid a false redistribution of stresses caused by springs subject to a tensile load, which would avoid the separation between the footing and the ground, the calculation was programmed iteratively by modifying the modulus of elasticity of each aforementioned spring to an almost zero value. A Matlab [79] subroutine was created to perform a matrix calculation of the footings and to identify the elements representing a terrain subjected to tensile stresses. Then, an almost zero cross-section area was assigned to those elements, which effectively deleted the previously considered spring, and a new matrix calculation was carried out. Once again, the recalculated elements subjected to tensile stresses were singled-out, their area-value was changed to virtually zero and a subsequent matrix calculation was implemented. This loop was repeated until the number of compressed springs remains stable after two consecutive iterations (i.e., the end of the calculation is reached when no variations in the number of springs subjected to tensile loads is detected after two iteration cycles).
From the results arising from the springs, it was possible to ascertain the pressure generated by the concrete footing on the ground, as well as the specific stress redistribution generated by the lateral walls and base of the element, which ultimately would allow for the calculation of the greenhouse foundation.
The model was studied to reduce its computational cost regarding the number of perimeter elements as well as in height. The optimal solution was reached for six sectors in height and four sectors at the horizontal plane. It should be noted that Figure 13 is a simplified representation of the model and does not depict the optimum solution.
In order to verify the accuracy of the model, the maximum values of the reactions (i.e., vertical, horizontal and moment loads determined through the EN standard specifications, as well as the CFD calculations, see Section 3), but with the opposed sign, were applied as external loads. Firstly, each type of stress was applied separately. Then, the effect of a simultaneous application was studied and used to compare the Wrinkler model's results against those obtained through FEM.
In Figure 13a, the numbered spots (1-6) are the start and end-points of the lines drawn in the retrieval of the results in this model, which were employed in the comparison.

Solid Model with FEM
The FEM was programmed in Ansys [64] and developed with 3D solid elements, an elastic soil and a surface-to-surface contact model between the footing and the ground.
As shown in Figure 14a, the model includes the half-volume of a cylindrical concrete footing and a half cylinder volume of the surrounding ground. Since both the geometry and load are symmetrical in the vertical plane, it was possible to reduce the domain in half by following this approach. For both the footing and the ground, the mesh was generated through eight-node three-dimensional solid elements (solid45). Furthermore, to simulate the connection between the walls and base of the footing with the ground, a surface-to-surface contact model was established by means of conta173 and targe170 elements, which support the separation between both contact surfaces and which is in line with the method followed by Vidal et al. [84,85].
Similar to the previous section, in Figure 14, the numbered spots (1 and 3, 2 and 4 and 5 and 6) are the start and end-points of the lines (1, 2 and 3, respectively) used for the retrieval of results in this model that would serve as the validation basis of the matrix model.
Both the foundation and the ground were considered as linear elastic materials. On the one hand, the foundation was defined through a modulus of elasticity of 3·10 7 kPa and a Poisson coefficient of 0.2. On the other hand, to support the coherence between the modulus of elasticity of the soil and the subgrade reaction modulus employed in the Winkler model, a simulation was carried since both parameters are related. Thus, the modulus of elasticity was determined as the load required to generate a vertical deformation of 0.00127 m (Figure 15). A force was applied to a 0.5 m load plate (i.e., the same diameter as the footing) until a deformation of 0.05 inches was registered in the simulated soil with a Poisson coefficient of 0.3. Finally, the soil was characterized by a modulus of elasticity of 10,500 kPa and a subgrade reaction modulus of 38,400 kN·m −3 , which denotes a sandy soil type. symmetrical in the vertical plane, it was possible to reduce the domain in half by following this approach. For both the footing and the ground, the mesh was generated through eight-node threedimensional solid elements (solid45). Furthermore, to simulate the connection between the walls and base of the footing with the ground, a surface-to-surface contact model was established by means of conta173 and targe170 elements, which support the separation between both contact surfaces and which is in line with the method followed by Vidal et al. [84,85]. Similar to the previous section, in Figure 14, the numbered spots (1 and 3, 2 and 4 and 5 and 6) are the start and end-points of the lines (1, 2 and 3, respectively) used for the retrieval of results in this model that would serve as the validation basis of the matrix model.
Both the foundation and the ground were considered as linear elastic materials. On the one hand, the foundation was defined through a modulus of elasticity of 3·10 7 kPa and a Poisson coefficient of 0.2. On the other hand, to support the coherence between the modulus of elasticity of the soil and the subgrade reaction modulus employed in the Winkler model, a simulation was carried since both parameters are related. Thus, the modulus of elasticity was determined as the load required to generate a vertical deformation of 0.00127 m (Figure 15). A force was applied to a 0.5 m load plate (i.e., the same diameter as the footing) until a deformation of 0.05 inches was registered in the simulated soil with a Poisson coefficient of 0.3. Finally, the soil was characterized by a modulus of elasticity of 10,500 kPa and a subgrade reaction modulus of 38,400 kN·m −3 , which denotes a sandy soil type.

CFD Simulation for the 0° Wind Loads
For the so-called for EV-D22 greenhouse, Figure 16 shows the different external pressure coefficient values expressed as function of the ratio between the distance from windward side of greenhouse cross-section (S) and the total length of greenhouse cross-section (Smax). In addition, the figure includes both the wind tunnel experiment results reported by Kim et al. [44] and the values resulting from the CFD model simulation proposed in this research work.

CFD Simulation for the 0 • Wind Loads
For the so-called for EV-D22 greenhouse, Figure 16 shows the different external pressure coefficient values expressed as function of the ratio between the distance from windward side of greenhouse cross-section (S) and the total length of greenhouse cross-section (Smax). In addition, the figure includes both the wind tunnel experiment results reported by Kim et al. [44] and the values resulting from the CFD model simulation proposed in this research work.
As could be observed in Figure 16, the experimental and numerical results exhibit a good correlation. Therefore, it was ascertained that the CFD simulation allowed a correct prediction of the areas with external pressure coefficients positive and negative.
In line with the European standard [53], positive cp e values were observed in the windward façade for both the experimental and numerical approaches. Nonetheless, the CFD simulation reported a mean cp e of 0.53, whereas the values arising from the experimental test [44] and the standard [53] were greater, 0.72 and 0.6, respectively. In this regard, it is worth mentioning that the velocity and turbulence profiles employed by Kim et al. [44] and EN 13031-1 [53] were different.
For the so-called for EV-D22 greenhouse, Figure 16 shows the different external pressure coefficient values expressed as function of the ratio between the distance from windward side of greenhouse cross-section (S) and the total length of greenhouse cross-section (Smax). In addition, the figure includes both the wind tunnel experiment results reported by Kim et al. [44] and the values resulting from the CFD model simulation proposed in this research work. As could be observed in Figure 16, the experimental and numerical results exhibit a good correlation. Therefore, it was ascertained that the CFD simulation allowed a correct prediction of the areas with external pressure coefficients positive and negative.
In line with the European standard [53], positive cpe values were observed in the windward façade for both the experimental and numerical approaches. Nonetheless, the CFD simulation reported a mean cpe of 0.53, whereas the values arising from the experimental test [44] and the standard [53] were greater, 0.72 and 0.6, respectively. In this regard, it is worth mentioning that the velocity and turbulence profiles employed by Kim et al. [44] and EN 13031-1 [53] were different.
For the windward side of the gabled-roof, the CFD results correctly predicted the change from wind pressure to suction due to the edge of the eave. Although the sign change began with a large suction figure, the cpe value was stabilized around the middle of the run. Once again, a good For the windward side of the gabled-roof, the CFD results correctly predicted the change from wind pressure to suction due to the edge of the eave. Although the sign change began with a large suction figure, the cp e value was stabilized around the middle of the run. Once again, a good correlation was found for the set of values. The CFD simulated mean cp e reached −0.22 and the wind tunnel test yielded a cp e of −0.3 [44]; a cp e of −0.2 has been established by EN 13031-1 [53].
Similarly, the CFD simulation also reported the existence of suction on the leeward side of the gabled-roof with a mean cp e value of −0.5. In this case, the EN standard [53] indicated the same −0.5 cp e value, However, the greater discrepancy was found with the wind tunnel test results as the cp e value reported by Kim et al. [44] was −0.61. Once correctly verified the proposed model (y+ < 1, SST k-ω turbulence model, domain size and geometry) on the EV-D22 greenhouse, a new simulation was carried out taking into account the new geometry (i.e., multi-span greenhouse defined in Section 2) and the different velocity and turbulence profiles in the inlet according to the terrain category II in Eurocode 1-4 [57] (Section 2.3). Figure 17 illustrates the wind velocity distribution on the streamlines. A greater velocity was observed on the zone of connection between the windward column and the first arch. As expected, the wind, which was tangential to that point of contact, created a greater suction effect on this initial zone, especially on the higher part of the arch. Moreover, regions of low velocity wind and suction were apparent on the connection zones between arches. In this regard, the greater suction appeared on the assembly region of the third arch and the leeward column. Similarly, the CFD simulation also reported the existence of suction on the leeward side of the gabled-roof with a mean cpe value of −0.5. In this case, the EN standard [53] indicated the same −0.5 cpe value, However, the greater discrepancy was found with the wind tunnel test results as the cpe value reported by Kim et al. [44] was −0.61.
Finally, for the leeward façade, the CFD simulation found a mean cpe value of −0.47, which is in good agreement with the cpe of −0.3 established in EN 13031-1 (2001) and the cpe of −0.54 indicated in the research carried out by Kim et al. (2017b).
Once correctly verified the proposed model (y+ < 1, SST k-ω turbulence model, domain size and geometry) on the EV-D22 greenhouse, a new simulation was carried out taking into account the new geometry (i.e., multi-span greenhouse defined in Section 2) and the different velocity and turbulence profiles in the inlet according to the terrain category II in Eurocode 1-4 [57] (Section 2.3). Figure 17 illustrates the wind velocity distribution on the streamlines. A greater velocity was observed on the zone of connection between the windward column and the first arch. As expected, the wind, which was tangential to that point of contact, created a greater suction effect on this initial zone, especially on the higher part of the arch. Moreover, regions of low velocity wind and suction were apparent on the connection zones between arches. In this regard, the greater suction appeared on the assembly region of the third arch and the leeward column. The external pressure coefficients of both the roof and the lateral columns are represented in Figure 18. For the roof (Figure 18a), the cpe are expressed as a function of the ratio between the distance from windward side of greenhouse cross-section (W) and the width of the tunnel module (Warch). Similarly, the cpe values corresponding to the lateral columns (Figure 18b,c) of the greenhouse  The external pressure coefficients of both the roof and the lateral columns are represented in Figure 18. For the roof (Figure 18a), the cp e are expressed as a function of the ratio between the distance from windward side of greenhouse cross-section (W) and the width of the tunnel module (W arch ). Similarly, the cp e values corresponding to the lateral columns (Figure 18b,c) of the greenhouse are expressed as the ratio between the elevation on the column (H) and its total height (H max ). To favor the comparisons, the CFD calculated cpe values were discretized into the regions defined by EN 13031-1 [53] as could be observed in Table 3. Correspondingly, the 0°-direction wind loads were also defined for each region (Table 3) and represented graphically in Figure 19, as the average value of those obtained within each defined area.  To favor the comparisons, the CFD calculated cp e values were discretized into the regions defined by EN 13031-1 [53] as could be observed in Table 3. Correspondingly, the 0 • -direction wind loads were also defined for each region (Table 3) and represented graphically in Figure 19, as the average value of those obtained within each defined area.  The comparison between the cpe values from the EN standard [53] and the proposed CFD simulation (Tables 1 and 3, respectively) revealed similar situations for the lateral columns, regardless of the calculation method. For the windward column, the CFD model yielded a cpe of 0.4, while EN 13031-1 [53] shows a 0.6 cpe value. Both the EN 13031-1 [53] and the CFD results showed a −0.5 cpe value for the leeward column.
However, important differences were observed between the cpe values reported by the EN standard [53] and the proposed CFD simulation, specifically for the first and third tunnels. In Marveas [49], the author also discussed the differences found between the EN standard [53] and the CFD simulation carried out by Kwon et al. [54] for single span greenhouses.
The most significant discrepancy was found at the first point of contact between the wind and the roof. Contrarily to the cpe value of 0.3 indicated by EN 13031-1 [53], the CFD results pointed to a significant suction behavior on the first arch with a cpe of −1.2. This finding is in line with the values stated in Eurocode 1-4 [57] for single-span vaulted roof structures, which also show a cpe of −1.2 for the geometrical conditions of the greenhouse studied in this paper (f/d = 0.188 and h/d = 0.56). However, caution should be taken in this justification, as single-span and multi-span results cannot be extrapolated due to the effect of the wind recirculation on the roof. Nonetheless, the differences caused by the recirculation effect are hardly noticeable in the wind entrance area. In any case, it should also be noted that Eurocode 1-4 [57] established a different angle for the first roof region, i.e., 55° instead the 45° stated in EN 13031-1 [53], which was considered negligible.
In this regard, Blackmore and Tsokri [86] have previously mentioned the discrepancies between EN 13031-1 [53] and Eurocode 1-4 [57]. The wind tunnel experiments carried out by the authors, who conducted their research on single-span curved roofed buildings (i.e., with dimensions not comparable to those frequently found in greenhouses), also indicated the existence of suction zones at the contact point of the wind with the roof. Furthermore, the authors detected a dual pattern of behavior for this initial zone, i.e., exhibiting both pressure and suction, which is in line with the information provided by Eurocode 1-4 [57] for duopitch roofs. Thus, it could also be appropriate to consider these two alternatives for multi-tunnel greenhouses.
For the superior part of the arch (55-115°), both the standard [53] and the simulation pointed to negative coefficients of −1 and −1.3, respectively. Similarly, suction was also detected in the remainder region of the first arch. A cpe value of −0.9 was reported by the simulation compared to the cpe of −0.4 The comparison between the cp e values from the EN standard [53] and the proposed CFD simulation (Tables 1 and 3, respectively) revealed similar situations for the lateral columns, regardless of the calculation method. For the windward column, the CFD model yielded a cp e of 0.4, while EN 13031-1 [53] shows a 0.6 cp e value. Both the EN 13031-1 [53] and the CFD results showed a −0.5 cp e value for the leeward column.
However, important differences were observed between the cp e values reported by the EN standard [53] and the proposed CFD simulation, specifically for the first and third tunnels. In Marveas [49], the author also discussed the differences found between the EN standard [53] and the CFD simulation carried out by Kwon et al. [54] for single span greenhouses.
The most significant discrepancy was found at the first point of contact between the wind and the roof. Contrarily to the cp e value of 0.3 indicated by EN 13031-1 [53], the CFD results pointed to a significant suction behavior on the first arch with a cp e of −1.2. This finding is in line with the values stated in Eurocode 1-4 [57] for single-span vaulted roof structures, which also show a cp e of −1.2 for the geometrical conditions of the greenhouse studied in this paper (f/d = 0.188 and h/d = 0.56). However, caution should be taken in this justification, as single-span and multi-span results cannot be extrapolated due to the effect of the wind recirculation on the roof. Nonetheless, the differences caused by the recirculation effect are hardly noticeable in the wind entrance area. In any case, it should also be noted that Eurocode 1-4 [57] established a different angle for the first roof region, i.e., 55 • instead the 45 • stated in EN 13031-1 [53], which was considered negligible.
In this regard, Blackmore and Tsokri [86] have previously mentioned the discrepancies between EN 13031-1 [53] and Eurocode 1-4 [57]. The wind tunnel experiments carried out by the authors, who conducted their research on single-span curved roofed buildings (i.e., with dimensions not comparable to those frequently found in greenhouses), also indicated the existence of suction zones at the contact point of the wind with the roof. Furthermore, the authors detected a dual pattern of behavior for this initial zone, i.e., exhibiting both pressure and suction, which is in line with the information provided by Eurocode 1-4 [57] for duopitch roofs. Thus, it could also be appropriate to consider these two alternatives for multi-tunnel greenhouses.
For the superior part of the arch (55-115 • ), both the standard [53] and the simulation pointed to negative coefficients of −1 and −1.3, respectively. Similarly, suction was also detected in the remainder region of the first arch. A cp e value of −0.9 was reported by the simulation compared to the cp e of −0.4 stated in the standard [53]. Except for the 80-100 • region of the second tunnel, a similar trend was observed throughout the greenhouse roof with greater cp e values arising from the simulation. This pattern was particularly apparent for the third tunnel as the EN 13031−1 [53] establishes these values as 60% reduction of those occurring in the second tunnel. Based on the CFD results, such an approach seems questionable as increases in the cp e values were detected compared to those of the second tunnel. In this regard, the simulation carried by Kim et al. [43] also reported cp e values closer to those of the central tunnel.
As expected, the second and third tunnels exhibited suction effects in the entirety of the arches. In addition, regardless of the calculation method, maximum suction was found on the crown region (80-110 • ) of the arch. Despite the numerical differences, both the standard [53] and the simulation resulted in similar cp e developments through the second and third tunnels.

Effects of the External Loads on the Greenhouse Structure
The influence of the external loads applied on the greenhouse were assessed through the λ cr value. Table 4 shows the λ cr values for the different combinations of actions considered. Load combinations numbered 1 to 20 correspond to those established by EN 13031-1 [53] and UNE 76210 IN [70] (Table 2). In addition, six additional load combinations were considered to also account for the effect of the 0 • -direction wind loads determined by CFD. Those are identified by a number (i.e., 1, 2, 7, 8, 11 and 12), which correspond to the number of the original combination, as per Table 2, and the CFD suffix. Moreover, Table 4 displays a concise description on the buckling mode shape resulting from the buckling mode, which could also be observed in Figures 20-22, for each specific load combination.
Agriculture 2020, 10, x FOR PEER REVIEW 20 of 32 43.86% were detected for load combinations 11 and 12 (i.e., combinations of the permanent load and the 0° wind load), respectively. For load combinations 1 and 1CFD, a global failure mode was observed with the brunt of the buckling pattern located within the windward arch ( Figure 20). Nonetheless, both load combinations could be set apart by the resulting shape. Whereas the wind load conforming to EN 13031-1 [53] resulted in a deformation that could be described as an up-down displacement throughout the arch (Figure 20a), the CFD wind load produced an opposed shape, which is a down-up displacement throughout the arch (Figure 20b). Nevertheless, both load combinations resulted in a significant reduction in the overall resistance of the structure as a consequence of the failure on the lower chord of the windward truss.     Conversely, the CFD wind (i.e., load combination 7CFD) generated a first arch buckling that caused the descent and subsequently ascent of the arch deformed geometry. A similar buckling behavior was again noticed when load combinations 8 and 8CFD, respectively, were considered. Thus, all the CFD results pointed to unsafe failure analysis when suction loads on the windward arch were not contemplated. characterized by the displacement of the top of all columns. Conversely, the CFD wind (i.e., load combination 7CFD) generated a first arch buckling that caused the descent and subsequently ascent of the arch deformed geometry. A similar buckling behavior was again noticed when load combinations 8 and 8CFD, respectively, were considered. Thus, all the CFD results pointed to unsafe failure analysis when suction loads on the windward arch were not contemplated.  Firstly, to assess the influence of the CFD wind loads in the buckling modes of the structure, the discussion was solely focused on the combinations including the 0 • -direction wind.
In general (Table 4), the use of the wind loads from CFD instead of the values established in the standard [53] resulted in significant λ cr decreases. The greater declines were observed for load combinations 1 and 2 (i.e., combinations of the permanent load, the 0 • wind load and the crop load). A 56.90% reduction was registered for load combination 1 and a 56.42% decrease was noticed for load combination 2. In the case of load combinations 7 and 8, which also grouped the permanent load, the 0 • wind load and the crop load, a lesser effect was found on the λ cr with reductions of 31.59% and 36.82%, respectively. The differences between the two sets of load combinations could be traced back to the 0.6 coefficient applied to the latter set of combinations. Finally, declining figures of 44.06% and 43.86% were detected for load combinations 11 and 12 (i.e., combinations of the permanent load and the 0 • wind load), respectively.
For load combinations 1 and 1CFD, a global failure mode was observed with the brunt of the buckling pattern located within the windward arch ( Figure 20). Nonetheless, both load combinations could be set apart by the resulting shape. Whereas the wind load conforming to EN 13031-1 [53] resulted in a deformation that could be described as an up-down displacement throughout the arch (Figure 20a), the CFD wind load produced an opposed shape, which is a down-up displacement throughout the arch (Figure 20b). Nevertheless, both load combinations resulted in a significant reduction in the overall resistance of the structure as a consequence of the failure on the lower chord of the windward truss.
As expected, load combinations 2 and 2CFD exhibited a similar buckling pattern ( Figure 21) to the aforementioned load combinations 1 and 1CFD (Figure 20), respectively.
It is worth mentioning that when the wind loads established by EN 13031-1 (2001) are responsible for the collapse of the lower chord on the windward arch, the displacements of the top of all columns, as well as the remainder nodes in all arches, were greater than those occurring in the CFD load combinations. Furthermore, the significant strength reduction of this lower chord generated a global buckling mode, but certainly closer to a failure that could be depicted as local (i.e., without displacement of the top of the columns) within the lower chord of the first arch. Figure 22 displays the buckling failure modes for load combinations 7 and 7CFD. The wind loads established in EN 13031-1 [53] (i.e., load combination 7) resulted in global buckling failures characterized by the displacement of the top of all columns. Conversely, the CFD wind (i.e., load combination 7CFD) generated a first arch buckling that caused the descent and subsequently ascent of the arch deformed geometry. A similar buckling behavior was again noticed when load combinations 8 and 8CFD, respectively, were considered. Thus, all the CFD results pointed to unsafe failure analysis when suction loads on the windward arch were not contemplated.
Both for the standard and the CFD wind load values, load combinations 11 and 11CFD, as well as 12 and 12CFD, exhibited the same buckling mode as in the aforementioned CFD load combinations. Similar to the pattern represented in Figure 22b, important displacements were registered in the first tunnel with the initial half of the arch moving downward as the second part moves upward. In addition, it is worth mentioning that, among all considered combinations, load combinations 11CFD and 12CFD displayed the lowest λ cr . The registered figures of 2.59 and 2.56, respectively, barely left any margin to the reduction between the λ cr and the α u . It is worth mentioning that load combinations 1CFD and 2CFD obtained the next lowest λ cr values.
Finally, for those combinations in which the lateral wind load was not of influence, the failure mode was described as general (Table 4) since generalized displacements were observed on the top of all columns, which corresponds to the same behavior as the one presented in Figure 22a. Among these combinations, the inclusion of the snow load (i.e., load combinations 19 and 20) resulted in the most unfavorable λ cr due to the closeness to the 3.6 limit value [53].
As established in EN 13031-1 [53], the stability study of the greenhouse arches could be carried out through fist-order analysis when λ cr ≥ 3.6. For the conditions studied in this paper, all the load combinations that incorporate the action load values corresponding to EN 13031-1 [53] complied the requirement. However, the consideration of the CFD wind load values resulted in lower values of λ cr than those imposed by the standard [53], which require a second-order analysis. This duality further reinforces the need for consideration of different wind loads in the behavior of multi-span greenhouse.
The analysis of the failure modes showed that buckling due to the compression of the arch and the lower chord of the truss on the first tunnel occurred easily for the CFD wind loads. Thus, for the conditions set in the present paper, the lack of consideration of the CFD wind suction loads under the 0-55 • zone of the first tunnel would trigger unsafe results.

Concrete Footings
The maximum force values at the base of the columns were extracted from the envelope of reactions arising from the load combinations and reached a maximum moment load of 5.53 kN·m, a maximum shear load of 2.61 kN and a maximum axial (compression) load of 7.58 kN, which was considered jointly with the self-weight of the footing. Note that the maximum values did not occur at the same column. Therefore, in a first approach, the effect of each force was considered separately on both the FEM and MM. Nonetheless, the influence of a simultaneous application of all maximum reactions was also studied in both models to assess the maximum stress exhibited by the ground. Figure 23 illustrates the stress results on the base of the footing (i.e., line 3 defined between points 5 and 6, as stated in Figure 14) determined by both models in the plane of symmetry. According to the MM results, the joint action of the axial force and the foundation self-weight are responsible for a 57.75 kN·m −2 constant value pressure under the footing. Contrarily, stress peaks up to 300 kN·m −2 are noticeable on the FEM. The stress concentration occurs near the edges in both lateral walls of the footing, which is common behavior in elastic models, known since the contact theory was developed, and previously reported [87,88]. Nonetheless, in the comparison between the normal stresses reported by each model, FEM values were lower due to the influence of the friction in the vertical walls on the transmission of the vertical stresses. Therefore, the stress on the base of the footing reported by the MM resulted in slightly greater values compared to both the real situation and the FEM results, which in turn could be considered an asset for the model as it lies on the safety side in terms of axial calculations.  Figure 24 shows the pressures on the ground arising from the application of a 2.61 kN shear load. The maximum pressure on the left wall (i.e., line 1 defined between points 1 and 3 as stated in Figure 14) was observed on the top and reached 32.01 kN·m −2 and 23.89 kN·m −2 for the MM and FEM, respectively. The difference between models (33.99%) points again to a conservative approach when employing the MM. On the opposed wall (i.e., line 2 defined between points 2 and 4 as stated in Figure 14), a sudden local stress increase could be observed at the base of the footing. Despite that both models detected the stress concentration effect, the maximum value registered by the FEM was greater (30.54 kN·m −2 ) than the one arising from the MM (15.36 kN·m −2 ), which in this case represent an underestimation (−49.70%) for the matrix method. Moreover, both the MM and the FEM distinguished the presence of another local stress peak at the bottom of the left-hand wall of the footing, in this case, signified by a zero-value pressure to the end of the line graph.
Finally, it is worth mentioning that the MM resulted in a lack of stress on the base of the footing (i.e., line 3 defined between points 5 and 6 as stated in Figure 14) since the self-weight of the foundation was disregarded and solely the shear load was applied. Nevertheless, an 18.51 kN·m −2 was registered by the FEM on the left-hand edge of the footing.  Figure 24 shows the pressures on the ground arising from the application of a 2.61 kN shear load. The maximum pressure on the left wall (i.e., line 1 defined between points 1 and 3 as stated in Figure 14) was observed on the top and reached 32.01 kN·m −2 and 23.89 kN·m −2 for the MM and FEM, respectively. The difference between models (33.99%) points again to a conservative approach when employing the MM. On the opposed wall (i.e., line 2 defined between points 2 and 4 as stated in Figure 14), a sudden local stress increase could be observed at the base of the footing. Despite that both models detected the stress concentration effect, the maximum value registered by the FEM was greater (30.54 kN·m −2 ) than the one arising from the MM (15.36 kN·m −2 ), which in this case represent an underestimation (−49.70%) for the matrix method. Moreover, both the MM and the FEM distinguished the presence of another local stress peak at the bottom of the left-hand wall of the footing, in this case, signified by a zero-value pressure to the end of the line graph.
Finally, it is worth mentioning that the MM resulted in a lack of stress on the base of the footing (i.e., line 3 defined between points 5 and 6 as stated in Figure 14) since the self-weight of the foundation was disregarded and solely the shear load was applied. Nevertheless, an 18.51 kN·m −2 was registered by the FEM on the left-hand edge of the footing. Figure 25 shows the pressures resulting from the isolated application of a 5.53 kN·m determined by both models in the plane of symmetry. A similar pattern to the one exposed in the shear load section was found. On the superior-end of the left wall (i.e., line 1 defined between points 1 and 3 as stated in Figure 14), the MM registered a pressure of 125.01 kN·m −2 , which lies again on the side of safety (31.34%) as the FEM reported a 95.18 kN·m −2 value. At the bottom of the right wall (i.e., line 2 defined between points 2 and 4 as stated in Figure 14), a peak stress was detected. Whereas a 248.34 kN·m −2 was determined through the FEM, the MM pointed to a lower value (−49.27%) and reached 125.98 kN·m −2 . Agriculture 2020, 10, x FOR PEER REVIEW 24 of 32  Figure 25 shows the pressures resulting from the isolated application of a 5.53 kN·m determined by both models in the plane of symmetry. A similar pattern to the one exposed in the shear load section was found. On the superior-end of the left wall (i.e., line 1 defined between points 1 and 3 as stated in Figure 14), the MM registered a pressure of 125.01 kN·m −2 , which lies again on the side of safety (31.34%) as the FEM reported a 95.18 kN·m −2 value. At the bottom of the right wall (i.e., line 2 defined between points 2 and 4 as stated in Figure 14), a peak stress was detected. Whereas a 248.34 kN·m −2 was determined through the FEM, the MM pointed to a lower value (−49.27%) and reached 125.98 kN·m −2 .
The stress concentration effect was also noticed in the base of the foundation (i.e., line 3 defined between points 5 and 6 as stated in Figure 14). A peak of 267.79 kN·m −2 was determined by FEM on the left-hand edge of the footing. Once again, it was revealed that an overestimation of the pressure results occurred for the FEM due to the stress concentration effect, which did not occur in the MM. Figure 26 illustrates the pressures generated against the ground when all stresses are considered simultaneously. On the superior-end of the left wall (i.e., line 1 defined between points 1 and 3 as stated in Figure 14), a pressure of 137 kN·m −2 was determined in the MM, whereas a value of 109.39 kN·m −2 observed through FEM, which constitutes a 25.24% deviation of the MM compared to FEM. Contrarily, on the inferior-end of the right wall (i.e., line 2 defined between points 2 and 4 as stated in Figure 14), the FEM results were greater (139.21 kN·m −2 ) than those found in the MM (120 kN·m −2 ). It should be noted that the highest FEM value occurred due to the effect of the stress concentration. Nonetheless, the difference between the models was lower (−13.80%). Finally, the base of the footing (i.e., line 3 defined between points 5 and 6 as stated in Figure 14) showed a greater stress concentration, 681.41 kN·m −2 according to the FEM compared to the 136 kN·m −2 determined in the MM, which represents a −80.04% variation.  The stress concentration effect was also noticed in the base of the foundation (i.e., line 3 defined between points 5 and 6 as stated in Figure 14). A peak of 267.79 kN·m −2 was determined by FEM on the left-hand edge of the footing. Once again, it was revealed that an overestimation of the pressure results occurred for the FEM due to the stress concentration effect, which did not occur in the MM. Figure 26 illustrates the pressures generated against the ground when all stresses are considered simultaneously. On the superior-end of the left wall (i.e., line 1 defined between points 1 and 3 as stated in Figure 14), a pressure of 137 kN·m −2 was determined in the MM, whereas a value of 109.39 kN·m −2 observed through FEM, which constitutes a 25.24% deviation of the MM compared to FEM. Contrarily, on the inferior-end of the right wall (i.e., line 2 defined between points 2 and 4 as stated in Figure 14), the FEM results were greater (139.21 kN·m −2 ) than those found in the MM (120 kN·m −2 ). It should be noted that the highest FEM value occurred due to the effect of the stress concentration. Nonetheless, the difference between the models was lower (−13.80%). Finally, the base of the footing (i.e., line 3 defined between points 5 and 6 as stated in Figure 14) showed a greater stress concentration, 681.41 kN·m −2 according to the FEM compared to the 136 kN·m −2 determined in the MM, which represents a −80.04% variation.    (Figure 27b).
According to the MM, the vertical displacement reached a value of −3.56·10 −3 m on both the superior and inferior left-hand points (i.e., points 1 and 3 in Figure 27b, respectively). However, the FEM results showed deformations of −0.282·10 −3 m and −0.246·10 −3 m for points 1 and 3, respectively. A similar trend was observed on both the superior (i.e., point 2 in Figure 27b) and inferior (i.e., point 4 in Figure 27b) right-hand points. The displacement values reported by FEM were 4.39·10 −5 m and 7.71·10 −6 m, respectively, whereas a significantly greater figure (6.23·10 −3 m) was determined through the MM.
The aforementioned results suggest that the matrix model based on the Wrinkel model is a suitable method as well as safe option in the design of greenhouse foundations. Nonetheless, it is worth mentioning that the contact model employed was based on a standard frictional model in Ansys, which considers a zero normal pressure if separation occurs. The stiffness contact was adjusted to reach both exact and convergent results. Moreover, the MM constitutes an easier calculation alternative against more complex approaches such the solid-elastic contact model through the finite element method. the MM.
The aforementioned results suggest that the matrix model based on the Wrinkel model is a suitable method as well as safe option in the design of greenhouse foundations. Nonetheless, it is worth mentioning that the contact model employed was based on a standard frictional model in Ansys, which considers a zero normal pressure if separation occurs. The stiffness contact was adjusted to reach both exact and convergent results. Moreover, the MM constitutes an easier calculation alternative against more complex approaches such the solid-elastic contact model through the finite element method.

Conclusions
In this paper, the structural assessment of a multi-span greenhouse was carried out from multiple numerical approaches:

Conclusions
In this paper, the structural assessment of a multi-span greenhouse was carried out from multiple numerical approaches: 1.
The wind loads affecting a multi-span greenhouse, such as the ones defined in this research paper, were numerically simulated by means of CFD. In this regard, a greenhouse composed of three tunnels was proven to be of interest from a scientific point of view, since such a geometry allowed to capture differences among the arches on the wind loads, as well as the effects of the loads on the buckling mode of the structure.

2.
For the columns on the greenhouse structure, the consideration of the 0 • -direction wind load led to similar results to those of the corresponding wind action stated in EN 13031-1 [53]. Nonetheless, significant discrepancies were detected for cp e values on the arches of the first and third tunnel.

3.
Conversely to the EN 13031-1 [53], the CFD wind loads resulted in suction effects on the initial windward zone of the first tunnel of a multi-span greenhouse. In this regard, the consideration of both pressure and suction effects for this region (0-55 • ) seems appropriate to analyze multi-tunnel greenhouses with this type of loads, which have been also supported by the Eurocode 1-4 [57] regulations and the results found in the literature.

4.
Regarding the third tunnel, the EN standard [53] establishes the cp e as the 60% value of those in the second tunnel. However, the CFD simulation did not registered the mentioned behavior but exhibited cp e values closer to those in the second arch. 5.
The CFD wind load values were also employed in the structural assessment of the structure. Due to the increase in the compression axial applied to the arches and the lower chord of the windward tunnel, decreases of the λ cr ranging from 31.59% to 56.90% were registered for those combinations of actions in which the 0 • -CFD wind loads were involved. 6.
In this regard, the application of the CFD wind loads would affect the stresses and strain of the structure, but also to the calculation method required to verify the stability of the arches. Since the λ cr for those load combinations did not exceed the 3.6 requirement established in EN 13031-1 [53], a second-order analysis should be performed instead of the first-order analysis required for the load combinations that only include the actions established in EN 13031-1 [53]. 7.
The reported declines in the λ cr due to the consideration of the CFD wind loads were responsible for global buckling modes, which were certainly closer to a typical local buckling shape. A comparison between the influence of the wind loads according to EN 13031-1 [53] and the CFD simulation was carried out for several load combinations in order to evaluate the differences in the resulting failure modes. For those combinations in which the lateral wind load was not of influence, the failure mode was described as general buckling and characterized by the generalized displacement of the top of the columns. 8.
Finally, a simple calculation approach (i.e., similar to the typical matrix structural analysis calculation procedure followed by any commercial structural software) was proposed for the determination of a non-superficial greenhouse foundation. Cylindrical concrete footings were simulated based on a large rigid block supported on the sides and bottom based on the Winkler model.
The results, both stresses and deformations, arising from the proposed method were conservative compared to those obtained by a finite element calculation method with contact in an elastic medium.
Regarding the stress determination, the values obtained through MM are more in line with the analytical results arising from the calculation of foundations, which is especially noticeable for the axial load. The presence of peaks on the base of the footing pointed to a stress concentration and skewed the FEM results due to the consideration of the elastic model instead of an elastoplastic model. The deformation results obtained through the MM are greater compared to those arising from the FEM; therefore, the MM is more conservative.