Equivalent Shell Model of Elastic Gridshells Including the Effect of the Geometric Curvature

: In this work, an equivalent continuum of a barrel gridshell is introduced. Constitutive identiﬁcation procedures based on periodic homogenization are provided in the literature for this purpose, based on a ﬂat Representative Element Volume (REV), notwithstanding that the geometry of the structures concerned is curved. Therefore, the novelty of the present study is the selection of a curved REV to obtain the equivalent elastic constants. The numerical validation of the identiﬁcation procedure is made comparing gridshell response to that of the equivalent shell under homogeneous load conditions. Finally, in order to highlight the effect of the curved geometry on the constitutive law of the continuum, the response of the proposed model is also compared to that of a continuum obtained from a ﬂat REV.


Introduction
A gridshell has a doubly curved geometry as a shell; however, its material is laid out in a grid pattern and that is the reason why it is also called a lattice shell or reticulated shell.
This kind of structure may be analysed as systems of beam elements, from now on referred to as "beam lattice" models, through extensions of the ordinary theories for structural frames, or alternatively they may be analysed as continua.
Although the propagation of highly efficient algorithmic tools for the analysis of the beam lattices is constantly growing (FEM and DEM software [1,2] or Cad tools endowed with Isogeometric analysis [3]), the deeper understanding of their mechanics is increasingly difficult to achieve. As a result, the designers are not protected against errors originating from the inability to check the solution of automatic computations, especially at the initial phases.
Analysing the static behaviour of beam lattices in parallel with that of their equivalent continua may be a solution to this problem.
In the second half of the twentieth century, the attraction for this approach led to defining different methods for establishing beam lattice-continuum relations. One of them is the equivalent stiffness method, conceived by Wright in [4]; the procedure provides the properties of the material and the effective thickness of the continuum model as a function of the geometry and the mechanical properties of the beam lattice shell, paying attention only on the membrane behaviour. Another method, called split rigidity method and introduced by K. P. Buchert in [5], is based on two different rigidities, one for the membrane deformations and the second for the bending ones. A further method is the orthotropic equivalent continuum method introduced in [6], based on an orthotropic constitutive law to find the rigidities of the equivalent shell.
The same results found in [6] are obtained by the authors of the present work in [7,8], where the orthotropic equivalent continuum, used to study the consequences of different orientations of gridshell laths, is obtained by homogenization using a flat REV.
One of the first validations of the equivalent continuum defined by the equivalent stiffness method has been carried out in [9]. In particular, the study compares the accuracy of the continuum for describing the buckling phenomenon and it concludes that the equivalent shell model could be employed only as preliminary tool to estimate the critical buckling load of a lattice shell.
Defining an equivalent continuum together with its validation is the dominant paradigm nowadays, as also shown in other application fields, such as periodic brickwork [10,11] and, recently, in nano-structures [12,13].
It is evident that the continuum model provides a less accurate description of the behavior of a gridshell than a beam lattice model. In this sense, using the denominations coarse and fine, respectively, to indicate the two models is justified. However, even if the continuum (coarse) model is less accurate than the beam lattice (fine) model, the former can be used if the loss of accuracy is compensated for by a considerable increase in usability [14].
Therefore, the continuum model is not considered as an alternative to the beam lattice one, but rather as a complementary tool; it is highly preferred in the first phase of a design, when the effectiveness of the global geometry of the gridshell has to be evaluated, an overall view of the distribution of internal forces is needed, the topologies of the grid and its orientation are still to be decided, and approximate values of bars forces are looked for.
The aim of this work is to provide a proper definition of the equivalent stiffness describing the transition from the fine model to the coarse one, that is, from the beam lattice model to the continuum shell. The aforementioned methods are addressed to beam lattice shells approximating shell surfaces within a framework of relatively short linear pieces; therefore, they do not consider the influence of the curved geometry. However, the equivalent continuum defined in [4] has also been used in the initial design phase of the Mannheim Multihalle, a timber gridshell assembled from a lattice of curved continuous laths [15]. That was justified since the ratio of the mesh size over to the sgrid span is sufficiently small.
However, how small does that ratio have to be so that the equivalent model is reliable? The present authors look for an answer to this question introducing, in the context of linearized elasticity, a new equivalent continuum for a simple barrel gridshell, deriving the constitutive coefficients through an identification procedure taking into account the effects of the curved geometry. Generally speaking, the identification of curved structures runs into complexities that depend on the value of the Gaussian curvature. Addressing the problem of curvature for the first time, having worked on the same theme in the past by using a flat REV, what we are doing now is to take into account the curvature on an extremely simple case, characterized by null Gaussian curvature. The geometric choice is accompanied by a further simplification, which limits the analysis to the coupling between axial and bending stresses, focusing attention only on symmetrical homogeneous states. Obviously, when we used the flat REV these coupling effects were neglected a priori.
To recap, we identify only some constitutive coefficients, with the ambition, on the one hand, to verify the effects of unidirectional curvature and, on the other hand, to focus on the aforementioned coupling. This identification is the subject of Section 2.1. Then, the numerical validation of the identification procedure is made comparing beam lattice response to that of the equivalent shell under homogeneous load conditions in Section 3.1. In order to show the actual influence of the curved geometry in the continuum formulation, the obtained equivalent continuum is finally compared, in Section 3.2, with the one proposed in [6], which neglects the curvature of the reference surfaces. The conclusions follow.

Methodological Premise
To determine the elastic coefficients of the continuum model equivalent to the beam lattice shell, a continualization procedure of a periodic system is used, which follows the lines of the homogenization methods of heterogeneous continuous models.
This procedure starts from the selection of an elementary reference volume (REV) in the fine model (the beam lattice), determined by the periodicity of the system, and identifies in the coarse model (the continuum) a part that occupies the same space region of the REV.
At this point, a point-wise correspondence is defined between the internal forces in the two models, coarse and fine ones, though a localization operator and, then, equality is imposed between the respective strain energies, expressed in terms of complementary energy. The correspondence between fine and coarse inner forces allows the strain energy of the fine model be expressed in terms of the continuous inner forces; the elastic coefficients of the continuum model are finally obtained by successive differentiations.
The determination of the strain energy of the fine REV passes through the resolution of a series of hyperstatic elastic problems, to which external forces of the bulk type and force boundary conditions are imposed. These problems are solved by the force method, after the imposition of constraints sufficient to remove the rigid motion. The simplicity of the beam lattice model allows the exploitation of elastic solutions in the literature ( [16]), which provide the displacement values in analytical terms. For reasons of simplicity, therefore, the evaluation of the strain energy of the fine model is carried out by estimating the external work, using the principle of virtual works.

The Fine and the Coarse Models: The REV
The beam lattice REV of a barrel gridshell is shown in Figure 1 on the left. The selected gridshell has a quadrangular mesh, whose members are modelled as Bernoulli beams (see [16] for equations). In particular, members can be divided into two groups: the longitudinal ones, that are straight beams, and the transversal ones, that are circular arches. Arches and beams are connected through panthographic joints at the intersections, in order to guarantee the continuity of the members. Finally, the quadrangular grid is triangulated by truss elements to provide in-plane shear stiffness. The model is immersed in a three-dimensional space, described by a global coordinate system xyz. In addition, each member is endowed with the own local coordinate system e 1 e 2 e 3 ; in particular, e 1 is tangent to the centroid axis of the member and e 2 and e 3 are the principal axes of inertia of the cross-section. Global and local coordinate systems of any member of the gridshell module are shown in Figure 1.
For the coarse model, a Donnell shell is selected, whose equations are described in [17] and recalled later. Even the continuum model is endowed with an own local coordinate system e 1 e 2 e 3 ; the first two local axes belong to tangent plane to the surface while the third one is perpendicular to it ( Figure 1). Local axis e 1 is tangent to the curved generator of the cylindrical shell, while e 2 is parallel to the straight one.
The identification procedure starts from the choice of a Representative Elementary Volume (REV) in the beam lattice model and the selection of a region in the continuum model of the same area as that occupied by the REV. Here, the module of the beam lattice shown in Figure 1 on the left is chosen as REV. The geometry of the REV is defined as having the length of the longitudinal straight beams and of the arches equal to l. Moreover, the arches are circular with radius R, and center angle 2β = l R . The same figure shows the corresponding continuum region on the right. To recap, the REV is composed of six elements: two curved beams, two straight beams and two pin-jointed bars, i.e., the diagonals. In addition, all the beams have been considered with half the area of their cross-sections and half the second moment of area, for reasons of periodicity, since they are common to two contiguous modules.

The Localization Operator
The objective is defining the constitutive law which describes membrane and bending behaviour of the continuum model, so that it is capable of providing strain states equivalent to the beam lattice model ones, under equivalent stress states. This is done by equalizing the internal work of the fine model to the internal work of the coarse one, obtained from the equivalent stress states. Figure 2 shows the uniform membrane and bending stresses acting on the sections of normal vectors e 1 and e 2 , which we collect in the vector s T = {N 11 , N 22 , M 11 , M 22 } T . Likewise, membrane and bending strains along the same directions are collected in the vector e T = { 1 , 2 , χ 1 , χ 2 } T . Their relationship is supposed to be linear, that is: where D is the flexibility matrix, whose coefficients are the unknowns of the identification procedure. In very general terms, it is written as follows: At this point, the continuum stress states, depicted in Figure 2, are made equivalent to the stress states of the beam lattice model and shown in Figure 3. By means of the following equations: where the radial pressure p on the arches corresponds to the radial bulk force p/l in the coarse model. This equation, often referred to as localization operator in standard homogenization methods, induces four elastic problems on the REV, shown in Figure 3: these are traction problems, due to the force identification procedure used up to now. The four load conditions are the following: 1. Uniform pressure p applied to the arches in the plane e 2 , e 3 , in equilibrium with the tensile forces pR at the end sections of the arches ( Figure 3a); 2.
Tensile forces F applied at the end sections of the straight beams ( Figure 3b); 3.
Couples m 1 applied at the end sections of the arches ( Figure 3c); 4.
Couples m 2 applied at the end sections of the straight beams ( Figure 3d).
To any traction problem kinematic constraints are added to eliminate any rigid motions and for reasons of periodicity. It is hardly worth pointing out that the end sections of arches are located on planes of symmetry of the fine model, on which all antisymmetric quantities, such as rotations, axial displacements and shear forces cancel out, for obvious reasons. Another plane of symmetry for the REV passes through the two center points of the straight beam axes (see Figure 4). These geometric conditions impose specific restrictions on the displacements, which are defined in Section 3, where, as those problems statically indeterminate, the force method is used to solve them.

Work Equality and the Elastic Coefficients
The solution of the elastic problems defined in Section 2.1.3 is presented in some detail in Section 2.2 and allows the strain energy of the fine model be defined. As mentioned in the methodological premise, in this case it is easier to determine the strain energy through the work done by the external forces on the corresponding displacements. External forces are those defined by the localization operator Equation (3). The only bulk force is the radial pressure, which acts in the plane of the arches and spends work for the radial displacement w of the axis line, while the remaining forces and torques spend work on the displacements and the rotations of the four nodes A, B, C and D. Therefore, the external work takes the following expression: where u is the displacement tangent to the axis line of the arches' end sections, v the axial displacement of the straight beams' end sections, and φ and θ are the components of nodes' rotations around the axis e 2 of the arch and the axis e 2 of the beam, respectively. Symmetry reasons impose the equality in norm of the homologue displacement components of the four edge nodes and transform the Equation (4) into the following: which in turn, by replacing the Equation (3) in it, is modified as follows: At this point the solutions of the elastic problems, defined in Section 2.1.3 and which are being treated extensively in Section 2.2, can be used to simplify Equation (6). First of all, strict periodicity reasons impose u = 0, canceling the second addend of Equation (6).
In addition, anticipating fragments of the solution, we have: with obvious meaning of the symbols used. Thus, Equation (6) transforms into: which, after some algebra, turns out to be a quadratic form of the coarse stress state The external work of the fine model W e f is ready to be equated with the internal work of the coarse model, given by where e, s and D are defined in Section 2.1.3. Taking into account Equations (9) and (10) assumes the following form: with the elastic constants D 11 , D 12 , D 13 , D 22 , D 23 , D 33 , D 44 given by: It is hardly worth noting that the constitutive coefficients in Equations (12) (7) do not have the physical dimension of displacements and rotations, but of derivatives of generalized displacements with respect to the internal forces of the coarse model, N 11 , . . . , M 22 , and, therefore, they contain the constitutive coefficients of the fine model, as is going to be shown later. Moreover, this is the reason for some apparent contradictions, like the one in Equation (16).
In addition, the mechanical couplings that can be deduced from Equations (13), (15) and (16) are worth to be underlined. The coefficient D 12 , coupling N 11 and N 22 , takes into account the Poisson effect in the membrane behavior, while the coefficients D 13 and D 23 , coupling N 11 and M 11 , and N 22 and M 11 , represent the answer to one of the questions raised by this paper. Furthermore, if in the plane of the arch the coupling measured by D 13 is practically discounted, the coupling in the perpendicular plane, highlighted by D 23 , is much less obvious.
This last coefficient represents the fact that when M 11 , while bending the arch in its plane, it lengthens its axis line, at the same time it also causes the elongation of the diagonals which, in turn, compress the longitudinal beams. On the contrary, the bending moment M 22 produces effects that remain confined in its plane, because since it does not lengthen the longitudinal beams, it does not lengthen the diagonals either, which are the coupling element of the two perpendicular planes.
The detailed expression of the elastic constants is given in Section 2.4. A final consideration is work making. The elastic constants determined above make Donell's shell model equivalent to the barrel beam lattice under examination. Having used a force-based approach, the stress states of the two models correspond point-wise via the localization equation Equation (3), while the deformation of the coarse model provides only an average measure of the deformations of the fine model. In fact, it is well known that the deformed shape of the fine model is characterized by local phenomena, which are lost in the coarse model.
For the determination of the displacements of the fine model, starting from the solutions of elastic problems of the equivalent Donnell shell, Equation (7) is used.

The Solutions of the Fine Model
The analytical solution of the four elastic problems defined on the REV can be easily found by exploiting the symmetries shown in Figure 4. These are due to the symmetry of the REV, with respect to the vertical planes π a , π c together with the global cylindrical symmetry that makes the circumferential components of the nodes displacements null.
In particular the symmetries allow to solve the whole problem by: 1.
Solving the sub-problem consisting of one half arch, one half diagonal and one half beam converging in the node B (see the rectangular hatched area in Figure 4); 2.
Assuming that the displacement of the node B belong to the plane π b .

The Solution Method
The symmetry and periodicity conditions allow each of the four 3D problems be solved as a linear combination of plane 1D problems. The 3D structure is statically undetermined with one degree of redundancy. The Force Method is used, where the axial force in the diagonal is assumed as redundant reaction. Each problem is then solved as the superposition of five elastic 1D plane statically determined sub-problems, by making use of the following scheme, where the index i = 1, . . . , 4 represents the quantities referred to the ith problem, while the apices a, b, d are referred to (arch, beam, diagonal), respectively: • The axial force in the diagonal is assumed as redundant reaction X i (see Figure 5). • The redundant reaction X i is projected onto the vertical plane π a containing the arch and in the plane π b containing the beam, by obtaining the components X a i = X i cos(α) and X b i = X i sin(α), respectively, where α is the angle formed between the diagonal and the plane π a in the reference configuration. The compatibility condition between the axial displacement of the diagonal and the projection of the axial displacements of arch and beam on the x-y plane is used in order to find the value of X i .
The whole displacement field for each element is finally written in terms of the external load only. Problem (1) Problem (2) Problem (3) Problem (4) REDUNDANT X PRIMARY STRUCTURE

The Equations of the Three Plane Problems
As well known, the elastic problems for the arch, beam and diagonal are defined by the following systems of equations: Elastic problem for the arch: Elastic problem for the beam: Elastic problem for the diagonal: The equilibrium equations are solved as a superposition of some basic solutions (BS) (see Figure 6) of statically determined problems.
The used basic equilibrium solutions for the arch (BS ai ), for the beam (BS bi ) and for the diagonal (BS di ) are the following:     (1) In this case on the arch both external load and X a 1 are applied, while on the beam only X b 1 is applied, then the basic equilibrium solutions to be used for the arch are BS a1 and BS a2 , by putting H = −X a 1 , while BS b1 must be used for the beam, by putting H = −X b 1 and BS d1 for the diagonal, by putting H = X 1 . Then, by integrating the compatibility equations for arch, beam and diagonal, we obtain the following results for the displacement fields (see Figure 7 for deformed configuration):

Solution of the Problem
Equation (19) is then solved to find the redundant reaction X 1 . Figure 7. Deformed configuration of the problem (1).

Solution of the Problem (2)
In this case on the arch only X a 2 is applied, while on the beam both F and X b 2 are applied, then the basic equilibrium solutions to be used for the arch is BS a2 , by posing H = −X a 2 , while BS b1 must be used for the beam, by posing H = F − X b 2 and BS d1 for the diagonal, by posing H = X 2 . Then, by integrating the compatibility equations for arch, beam and diagonal, we obtain the following results for the displacement fields (see Figure 8 for deformed configuration): Equation (19) is then solved to find the redundant reaction X 2 .

Solution of the Problem (3)
In this case on the arch both m 1 and X a 3 are applied, while only X b 3 is applied on the beam, then the basic equilibrium solutions to be used for the arch are BS a2 and BS a3 , by posing H = −X a 3 , and m = m 1 , respectively, while BS b1 must be used for the beam, by posing H = X b 3 and BS d1 for the diagonal, by posing H = X 3 . Then, by integrating the compatibility equations for arch, beam and diagonal, we obtain the following results for the displacement fields (see Figure 9 for deformed configuration): Equation (19) is then solved to find the redundant reaction X 3 .

Solution of the Problem (4)
In this last case the redundant reaction X 4 vanishes (see Figure 10 for deformed configuration).

The Identified Constitutive Coefficients
Once the four fine problems are solved, it is possible to put solutions in Equations (12)-(18) by means of the following substitutions: Then, by substituting β = l 2R , it is possible to write all geometrical parameters as functions of β: Finally, the identified constitutive coefficients, as functions of the angle β can be represented as follows: where the dimensionless coefficients c 1 (β), c 2 (β), c 3 (β) are defined as follows: It easy to show that: Consequently, taking the limit of the identified coefficients for small angles β, we obtain the same coefficients introduced in [6]: Figure 11 shows the behaviour of the identified coefficients as functions of β, superimposed with the value obtained in [6]. It is worth noting that: • For β = 0 (flat REV) the values coincide with that of [6].

Homogeneous Load Case
In order to validate the identification procedure, the response of coarse model is numerically compared to that of the fine one.
The numerical analyses are carried out with the finite element software SAP2000 [18]. A cylindrical gridshell having a 1.90 m radius is chosen as fine model, divided in 36 modules (six longitudinal and six circumferential, see Figure 12). The characteristic length of each module is 1.00 m. All structural elements are modelled as timber laths having a rectangular cross-section, 50 mm wide and 30 mm deep; the chosen material belongs to the strength class D30 [19], with the mean modulus of elasticity parallel to the fibers equal to 11,000 MPa. The fine model has a global orientation as defined in Figure 1. It is analysed under two homogeneous load conditions: a uniform pressure load distributed on the arches (Figure 12a) and a distributed moment about the axis of the straight edge beams (Figure 12b).
Boundary conditions (b.c.) are assigned with reference to the displacement components defined in Section 2.4. In both cases of Figure 12a,b all the nodes of the gridshell depicted as a circle have b.c. := {u, φ, θ} = 0 ∪ {v, w} = 0. However, in order to avoid the rigid translation in y direction, the arch where nodes are depicted with squares have b.c.
Coarse model is analysed under equivalent load and boundary conditions; the results in terms of deformed shapes for both models are shown in Figure 13. In the first case (a), the coarse model behavior can be characterized just by 1 and 2 . These strain measures are constant all along whole shell; in particular, 1 = 1.02 × 10 −4 and 2 = −1.38 × 10 −5 . By applying Equation (7), the w component of the displacement vector is calculated equal to 6.4 × 10 −2 m. Thus, the percent error between the coarse model and the fine model is below 1%.
In the second case (b) instead, the strain measures χ 1 , 1 and 2 characterize the coarse model behavior, and their numerical values are calculated: χ 1 = 0.20 m −1 , 1 = −1.03 × 10 −3 , rotation vector, this is evaluated as 0.102 rad. However, the obtained value of φ corresponds to the rotation of the arches' end sections of one module.
Being χ 1 constant in the problem of Figure 12b it can be deducted that rotation φ is linear. Thus, φ = 0.612 rad at the end section of the whole arch of the fine model when composed of six modules. By the numerical analysis of the fine model, φ = 0.622 at the end section of the arches; therefore, the error of the coarse model is below to 2%.
The v component of the displacement vector can be verified in the same way for both load cases; here, the error of the coarse model is again negligible.

Influence of the Curvature: A Non-Homogeneous Example
A gridshell can be analysed using continuum models obtained by different identifications. In this section, the focus is on the comparison of two continua obtained by constitutive identifications based on different REVs. The continuum model proposed in the present is now compared to that defined in [14]; for the sake of simplicity, the first one is named continuum A (CA) and the second one continuum B (CB).
The responses of the different models are compared in terms of mechanical work (W), numerically evaluated as the work spent by the external forces over the displacements.
A first comparison (shown in Figure 14) is done basing on the results of five tests, characterized by the same radius (R) but a different number of modules (n); in particular, the radius is 20 m and the number of the modules ranges from 12 × 12 to 36 × 36. Note that in this case β is fixed, having l and R fixed values.
Each model, namely the fine and the coarse, is analyzed under the same uniform vertical load equal to 1.0 N/m 2 . The shell is supported along the boundary straight edges by pins on one side and by rollers on the other side.   Figure 14 depicts the ratio (W c /W f ) of the mechanical work of the continuum model (W c ) over the mechanical work of the fine model (W f ). As could be expected, when the number of modules increases the response of both continua tend to that of the corresponding fine model. Although the trends are similar, the response of CA is closer to that of the fine model than CB's one. Since both curves take values greater than 1, the equivalent continua reveal to be more flexible than the corresponding gridshells.
A second comparison is performed in order to study the effect of the geometric curvature κ = 1/R. Seven numerical tests are performed on purpose, where the gridshells are subject to the same loads and b.c. as in the previous test. This time the number of modules is fixed to (n = 24), while the curvature κ ranges from 0 to 0.1125 m −1 (Figure 15). When κ is zero, i.e., flat geometry, the response of the two continua clearly provides the same error compared to the fine model. However, interestingly the responses of the two continua rapidly diverge as the curvature increases, as Figure 15 shows.

Conclusions
In this work, we demonstrated that a proper constitutive identification of curved gridshells leads to the expected coupling of membrane and flexural behavior for an equivalent continuum model. This result is due to a novel approach taking into account the geometric curvature in the equivalent continuum model. Even the geometrically simplest case of a barrel gridshell shows how much the solution can diverge while adopting an equivalent continuum derived from a flat REV.
As highlighted in the introduction, the transition from the fine model to the continuum one is a useful tool in a phase of generation and selection of the global form of gridshells, when a precise definition of the grid is probably useless and cumbersome. Therefore, the constitutive identification adopted in this work can be a useful tool for the preliminary design stage. Moreover, this work can be considered a first step towards a methodologically correct approach aiming at defining equivalent continua of gridshells with double curvature.

Conflicts of Interest:
The authors declare no conflict of interest.