Emilio Turco
email: emilio.turco@uniss.it

Anil Misra
email: amisra@ku.edu

Marek Pawlikowski
email: m.pawlikowski@wip.pw.edu.pl

AND FRANC ¸OIS HILD Francesco Dell'isola
email: francesco.dellisola@uniroma1.it

ENHANCED PIOLA-HENCKY DISCRETE MODELS FOR PANTOGRAPHIC SHEETS WITH PIVOTS WITHOUT DEFORMATION ENERGY: NUMERICS AND EXPERIMENTS

Keywords: April 17, 2018 Mechanics of materials, pantographic sheets, discrete models, nonlinear analysis

, is performed and shows that a remarkable agreement between theoretical predictions and experimental evidence can be obtained. This circumstance is easily explained by observing that said numerical code is based on a discrete model directly inspired by the mechanical properties of pantographic sheets and that, therefore, the passages to a continuum model via homogenization [8] and then to the subsequent re-discretization, via the introduction of more or less suitable finite elements, are avoided. In our opinion a theory driven formulation of a directly discrete numerical model presents many advantages and it seems suitable for attacking future structural optimization problems.

Introduction

One of the frontiers of research in mechanics must be situated at the border sharing the models introduced for describing standard materials and those for exotic materials. There are many difficulties in recognizing where such a frontier is located, especially because the adjectives standard and exotic are very difficult to make precise. A critical reader may even start discussing preliminarily about the most fundamental concept of material and the difficulties involving in its definition. Questions arise: Can a complex mechanical system as the one depicted in Figs. [START_REF] Alibert | Truss modular beams with deformation energy depending on higher displacement gradients[END_REF], 2 and 3 be called a material? And if it actually can be called so, when such a description is appropriate? The aforementioned system is constituted at the micro-level by matter distributed in a Figure 1. Rendering of the pantographic sheet frontal view, the unit cell is highlighted in the red box (courtesy of M. Golaszewski). refined and complex microstructure where, for instance, micro-gaps divide different deformable micro-parts, which, in some cases, may undergo large localized relative displacements.

The critical reader may maintain that if there is a material at all in such a system one has to find it at the micro-level, that is at the level at which the characteristic length is a fraction of the dimensions of the pivot (i.e. the small cylinders that link orthogonal beams) depicted in Fig. 3. To this reader another one, even more critical, may object that, by magnifying a little bit the image (see Fig. 4) another microstructure may appear, which is made of partially melt and agglomerated grains of polyamide powder used as raw material of the 3D printing process to produce the considered specimens whose theoretical analysis is reported in Ref. [START_REF] Milton | On the possible effective elasticity tensors of 2-dimensional and 3-dimensional printed materials[END_REF].

These critical remarks have often been made to try to understand the ultimate nature of physical phenomena [START_REF] Berryman | The Stanford Encyclopedia of Philosophy[END_REF][START_REF] Feynman | The Feynman Lectures on Physics[END_REF][START_REF] Dell'isola | The complete works of Gabrio Piola: Volume I -Commented English Translation[END_REF][START_REF] Dell'isola | At the origins and in the vanguard of peridynamics, non-local and higher-gradient continuum mechanics: An underestimated and still topical contribution of Gabrio Piola[END_REF]). Therefore, the term homogeneous material is defined relative to a length scale L of the corresponding Representative Volume Element (RVE), a cubic volume whose sides are L such that, by moving the RVE within the specimen the overall (macro-)mechanical response of the material included in it does not change and can be described exclusively in terms of overall (macro-)kinematic descriptors that are assumed to be uniform for every RVE. Polyamide, when considering an RVE including a set of grains, is a material in Figure 2. Technical drawing with the geometry of the pantographic sheet and the detail of the designed pivot, dimensions are expressed in mm (courtesy of M. Golaszewski).

the sense of the previous definition. The pantographic sheet, when considering an RVE including a set of cells (which is the set of four beam elements and four pivots shown in Fig. 1) is another material, although it is made ultimately, at a lower length scale, of another material (i.e. polyamide).

1.1. Definition of metamaterials. A metamaterial is defined herein as a material that has been designed to meet a specific purpose, by combining more elementary materials (characterized by a smaller micro-length scale) and by shaping them with geometrical structures and mechanical interactions (i.e. microstructure) characterized by the same micro-length scale. The micro level corresponds to the scale at which the considered structure shows all its (geometrical and mechanical) inhomogeneity and complexity. At the macro level, it behaves as a homogeneous material.

It is worth noting that the more interesting cases (i.e. the macro-metamaterial shows a completely different behavior as its micro-constituents) are represented by microstructures in which extremely high contrasts of mechanical and geometrical properties occur [START_REF] Pideri | A second gradient material resulting from the homogenization of an heterogeneous linear elastic medium[END_REF][START_REF] Boutin | Generalized inner bending continua for linear fiber reinforced materials[END_REF][START_REF] Alibert | Second-gradient continua as homogenized limit of pantographic microstructured plates: a rigorous proof[END_REF][START_REF] Dell'isola | Large deformations of planar extensible beams and pantographic lattices: Heuristic homogenisation, experimental and numerical examples of equilibrium[END_REF]. In this sense the metamaterial studied herein is extreme in all aspects, namely, very thin structural elements interconnect stiffer elements, some elements have zero stiffness, others have a soft behavior, others are much stiffer. The geometry is, also, extremely heterogeneous with empty spaces used to partially detach the elastic response of interconnected elements (this is, for instance, the function of perfect pivots).

In order to design a metamaterial one has to choose its governing equations and then find the microstructure which, at the macro-level, is described by given equations. This formulation is often called the problem of synthesis of a given set of governing equations. 1.2. Continuum versus discrete models. The ideas of continuum mechanics are used, in general, to model the deformation phenomena occurring in any kind of material body. Therefore, very often, the mechanical systems designed to get metamaterials are described both at the micro-and macro-levels with continuous models. In particular, the problem of the synthesis of second gradient (meta)materials has been confronted and solved with homogenization methodologies involving continuous descriptions both at macro-and micro-levels, see [START_REF] Eugster | Isola. Exegesis of the Introduction and Sect. I from "Fundamentals of the Mechanics of Continua" by E. Hellinger[END_REF][START_REF] Eugster | Isola. Exegesis of Sect. II and III.A from "Fundamentals of the Mechanics of Continua[END_REF][START_REF] Eugster | Isola. Exegesis of Sect. III.B from "Fundamentals of the Mechanics of Continua[END_REF]. The need for synthesizing a second gradient (meta)material (i.e. to design a suitable microstructure showing a second gradient behavior at macro-level) was formulated in a purely theoretical context in which even the consistency of mathematical theories was seriously questioned. This scientific controversy is not new. Piola did formulate mathematical models for continuous materials by introducing deformation energies depending on higher order gradients of displacement [START_REF] Dell'isola | The complete works of Gabrio Piola: Volume I -Commented English Translation[END_REF] and Peridynamics [START_REF] Dell'isola | At the origins and in the vanguard of peridynamics, non-local and higher-gradient continuum mechanics: An underestimated and still topical contribution of Gabrio Piola[END_REF]). When asked about the soundness of his continuum models and about the true physical content of his theories, he resorted to the study of a homogenization problem. He proved the validity of his macroscopic models by i) basing them on the study of a micro-structured mechanical system characterized by simply interacting elementary constituents and a well specified geometry, ii) describing such micro-system with a Lagrangian discrete finite dimensional model, iii) finally homogenizing micro-models (via the so-called Piola Ansatz) to get the corresponding macro-response [START_REF] Dell'isola | At the origins and in the vanguard of peridynamics, non-local and higher-gradient continuum mechanics: An underestimated and still topical contribution of Gabrio Piola[END_REF]).

Under the push for the need of new metamaterials, the paradigm set up by Cauchy seems to need an evolution, as already demanded by Piola in 1825 even before Cauchy's full formulation of his ideas. More generally one could attempt to solve a class of problems formulated as follows [START_REF] Dell'isola | Variational principles are a powerful tool also for formulating field theories[END_REF][START_REF] Alibert | Second-gradient continua as homogenized limit of pantographic microstructured plates: a rigorous proof[END_REF]). Given: a macro-length scale, set of kinematic macro-descriptors, some strain energy functional and Rayleigh dissipation functional; find: a micro-length scale, a set of micromaterials and a microstructure such that the complex micro-system behaves at the macro-level as specified by the given functionals.

Coming back to the particular problem of synthesis of second gradient materials by using architectured microstructures, it was solved with a double scale homogenization process [START_REF] Alibert | Truss modular beams with deformation energy depending on higher displacement gradients[END_REF][START_REF] Seppecher | Linear elastic trusses leading to continua with exotic mechanical interactions[END_REF]. Slender continua are considered and, after homogenization, beam elements are obtained, and subsequently interconnected with ideal pivots, and a second homogenization step is performed. Ideal pivots represent preferred constraints in the chosen class of micro-lattices that appear to be essential to produce second gradient macro-continua. The homogenized macro-equations show some pathologies that reflect the "exotic" behavior of considered metamaterials [START_REF] Boutin | Linear pantographic sheets. Asymptotic micro-macro models identification[END_REF][START_REF] Eremeyev | Linear pantographic sheets: Existence and uniqueness of weak solutions[END_REF]. These mathematical pathologies are of interest by themselves [START_REF] Eremeyev | Linear pantographic sheets: Existence and uniqueness of weak solutions[END_REF] as the theoretical issues that they provoke not only represent an intellectual challenge but also means for disclosing potentially interesting new phenomena.

In order to make such disclosure evident, and to experimentally investigate the special role that perfect pivots played in synthesized microstructures, the following technological demand arose: "to find an innovative design and production process to fabricate, via additive manufacturing, pantographic sheets made for instance of polyamide, whose pivots actually twist without deformation and with negligible dissipation" [START_REF] Golaszewski | Metamaterials with relative displacements in their microstructure: technological challenges in 3D printing, experiments and numerical predictions[END_REF]. The specimens could be printed in a monolith and required no post-assembly but only an easily standardized run-in procedure.

A potential vicious circle immediately appears under first examination of the previously described conceptual procedure involving continuous systems: i) as underlined by Piola [START_REF] Dell'isola | The complete works of Gabrio Piola: Volume II -Commented English Translation[END_REF] a clear understanding of the mechanical features of considered microstructures is attained mainly by formulating finite dimensional Lagrangian models, possibly with dissipation accounted for by a Rayleigh potential; ii) a homogenization procedure can be used to get continuum models but the range of applicability is restricted by the finite size of the periodic cell constituting the microstructure; iii) nowadays, to get predictions one needs to use a numerical code and therefore the continuum model must be discretized, for instance with a finite element scheme. One may wonder for which reason it is necessary to use the intermediate step involving homogenization and why one should loose all phenomenological control on the third step, i.e. discretization, which is usually performed based only on purely computational considerations. For this reason we have decided to completely refrain from the introduction of any continuum model for the microstructure which we consider in this paper and we preferred to introduce and study directly a discrete (Lagrangian) micro-model (see [START_REF] Turco | Isola. King post truss as a motif for internal structure of (meta)material with controlled properties[END_REF]) for a different example of similar discrete model inspired by trusses).

1.3.

A first discussion of the presented novel results. It is decided to avoid facing the pathologies presented by continuum models, bypassing their formulation, and the well-posedness of corresponding boundary value problems. Instead a finite set of Lagrangian parameters is proposed in order to specify the state of considered system. This set of parameters includes, in particular, the positions of the nodes corresponding to the twisting pivots interconnecting the two families of fibers (Fig. 3) constituting the sheet. It is assumed that such pivots are strong enough to ensure that the two beams they are interconnecting have a constant distance. When this assumption is not satisfied the reader is referred to the treatment proposed in Ref. [START_REF] Spagnuolo | Qualitative pivot damage analysis in aluminum printed pantographic sheets: Numerics and experiments[END_REF], which can be possibly modified to generalize the description presented hereafter. However, the mechanical case is sought in which the beams interconnecting the twisting pivots are designed to be long enough so that a non negligible bending energy is stored between the closest pairs of pivots. It is worth noting that a twisting pivot at the micro-level will produce, at macro-level, a shear stiffness.

In order to account for this additional deformation mechanism, which is not accounted for in Refs. [START_REF] Dell'isola | Large deformations of planar extensible beams and pantographic lattices: Heuristic homogenisation, experimental and numerical examples of equilibrium[END_REF][START_REF] Turco | Hencky-type discrete model for pantographic structures: numerical comparison with second gradient continuum models[END_REF], the discrete Piola-Hencky model is enriched by adding extra degrees of freedom. The discrete kinematic description is enhanced by distinguishing two types of nodes, namely, the first one, which was already considered, interconnects the two families of pantographic fibers while the second one simply interconnects two subsegments (Fig. 1) of the same fiber. The first family of nodes is where the inter-fiber interactions are located, while the second one is where the extra bending energy relative to the fibers is assumed to be concentrated. Therefore the overall considered set of kinematic parameters must also include the positions of these intermediate nodes, which are allowing for the bending of the fibers between two twisting pivots.

The Lagrangian method, the procedure for the search of minima of the total energy, and the step-wise construction of equilibrium configurations, described by one or more parameters, which was already exploited in Refs. [START_REF] Turco | Hencky-type discrete model for pantographic structures: numerical comparison with second gradient continuum models[END_REF][START_REF] Turco | Quantitative analysis of deformation mechanism in pantographic substructures: experiments and modeling[END_REF], allow for the determination of a large family of equilibrium shapes for perfect pantographic sheets. In particular, the Vietnam long neck 1 shape, which is experimentally observed in standard extension bias tests, is predicted with a very short computing times. Further, with a very limited set of constitutive parameters, experimental observations are predicted (e.g.quantitatively describing equilibrium shapes, reaction forces via application of Castigliano theorem, and the first rupture mechanism [START_REF] Turco | Fiber rupture in sheared planar pantographic sheets: numerical and experimental evidence[END_REF]).

The reason for which the aforementioned agreement between numerical predictions and experimental evidences had to be expected can be explained as follows. The numerical code is based on a discrete model that is directly constructed and closely inspired by the understanding of the deformation mechanism(s) of the studied pantographic sheets. There is no need for complex constructions but efficient numerical algorithms. This gives a further example of a theory driven analysis of complex phenomenology, whose study based on data driven approaches would produce an intractable and probably useless amount of computing processes.

After the state-of-the-art presented in this structured Introduction, we discuss the structural design of pantographic sheet having pivots without strain energy in Sec. 2. Successively, in Sec. 3, we report and discuss the results of an extension bias test. We tackle the presentation of the enhanced Piola-Hencky discrete model of the pantographic sheet in Sec. 4 whereas we report some insight about the in-house numerical code to attack the nonlinear system of equilibrium equations in Sec. 5 and the calibration of the numerical model in Sec. 6 by using the 1 see https://en.wikipedia.org/wiki/Neck/ring results of the extension bias test. In Section 3.1 we try to bridge experiments and numerical simulations introducing the tool of digital image correlation. Indeed, up to now, the experimental study of pantographic sheets and the comparison of experimental evidence with numerical predictions has been performed by means of simple tools and by focusing on the measurements of few (particularly significant) kinematical or dynamical parameters. However more systematic and accurate analysis of deformed shapes is required and the tools supplied by digital image correlation (DIC), see [START_REF] Sutton | Image correlation for shape, motion and deformation measurements: basic concepts, theory and applications[END_REF][START_REF] Hild | Digital image correlation. Optical methods for solid mechanics: a fullfield approach[END_REF], seem to be suitable to supply it. Some really suggestive images are presented whose interpretation (from quantitative and qualitative point of view) is supporting strongly the modeling choices presented in this paper. Finally, we close the paper reporting some remarks and future challenges.

Structural design: pivots without strain energy

Up to now, the pantographic sheets, which were fabricated with additive printing, necessarily had two families of fibers interconnected with small cylindrical pivots undergoing rotations at the expense of some non-negligible strain energy. This circumstance limits the possibility of structural optimization, which could require the choice of pivots having variable stiffness in different locations. Conversely, pivots without strain energy could lead to a change of the distribution of strain energy and its partition, e.g. into bending and extensional energies, thereby improving the damage resiliency of such structures.

The resolution of the printing process was 0.1 mm and it was exploited in order to leave un-sintered some grains of polyamide between the two adjacent solidified male and female parts of the printed pivots (Figs. 2 and3). The pantographic sheet shown in Figs. 1-3 is made of a series of parallel micro-beams linked to the orthogonal ones by means of the introduced pivots. Cross-sections of the micro-beams are rectangular with the exception of the parts near the pivots, which are designed to be linked with the orthogonal family of micro-beams by means of specifically tailored the pivots.

Extension bias tests: experimental evidence

The aforementioned first printed specimens (overall dimensions the specimen are 213.4 mm × 71.13 mm) were subjected to extension bias tests and the evolution of their deformed shape recorded by acquiring image series. For the first test, an additional preparation was followed to probe the feasibility of digital image correlation in the present context. For the second test, the grips were left unspeckled and no DIC analysis was run.

3.1. Digital image correlation analysis. In the sequel a preliminary image analysis is performed to quantify global deformation modes in an extension bias test. Digital image correlation (DIC [START_REF] Sutton | Image correlation for shape, motion and deformation measurements: basic concepts, theory and applications[END_REF][START_REF] Hild | Digital image correlation. Optical methods for solid mechanics: a fullfield approach[END_REF]) will be used to measure displacement fields during the reported experiment. To the authors' best knowledge, such technique has never been used in this context. A series of 42 pictures were acquired. In the present case, grey level images are reconstructed by a binning of two process to account for Bayer's filter. They are subsequently cropped so that their final definition is 1147 × 261 pixels. For DIC purposes, the upper and lower parts of the grips were speckled with black and white paints (Fig. 5). In addition, the pivots of the pantograph were marked in black and a red background was used in order to create high contrast with the white pantograph. The analysis of pantograph motions can be performed at different scales. The most obvious, yet the most difficult, would be to have a geometric description that would fit the actual pantograph shape. The corresponding kinematics could then be tailored thanks, for instance, to generalized beam theories. DIC analyses may then be possible as was shown for steel beams [START_REF] Hild | Displacement measurement technique for beam kinematics[END_REF]. In the present case, it would mean having an explicit description of all the fibers that compose the pantograph. This route will not be followed hereafter. Another path consists in measuring macroscopic motions via so-called local or global DIC [START_REF] Hild | Comparison of local and global approaches to digital image correlation[END_REF]. Since the final aim of such analyses will be their comparison with numerical simulations, the displacement fields will be expressed in the same language, namely, with finite element bases. Global DIC with 3-noded triangles with linear interpolation functions is considered (i.e. T3-DIC). Because the mesh is not made compatible with the pantograph mesostructure, mechanical regularization will be used to enable the DIC code to converge (i.e. RT3-DIC [START_REF] Tomičević | Mechanics-aided digital image correlation[END_REF]). To avoid any significant bias, the regularization length is identical to the element length (i.e. 25 pixels in the present case).

Figure 6(a) shows the region of interest (ROI) that extends over the whole pantograph and part of the speckled grips. This setting allows the RT3-DIC code to converge first in the grips and then in the central part of the pantograph where most of the deformation takes place. In the present case, the convergence criterion was set to 10 -5 pixel for the norm of the mean displacement corrections. This very low value could be achieved thanks to the regularization strategy used herein. Figure 7 shows the results of the RT3-DIC analysis of the 5th loading step. First, the grey level residuals should be checked, see Figure 7(a). They correspond to grey level differences between the picture in the reference configuration and that of the deformed configuration corrected by the measured displacement field. If perfect match were achieved, this difference should only contain acquisition noise. This is not observed in the present case since there still are some stigmata of the pantograph geometry. It confirms that a continuous kinematic basis made of 3-noded triangles cannot fully capture the present kinematics. However, the deviations remain very small. Consequently the displacement fields shown in Figure 7(b-c) are deemed trustworthy at the macroscale. The transverse u x displacement field shows that there is a huge contraction in comparison with the longitudinal u y motions. This is due to the geometry of the pantograph. This observation translates into the nominal strain components reported in Figure 7(d-f). In the present case, they are computed from the deformation gradient, which is constant per element. Thus the pixelisation of the strain maps is due to the underlying mesh used in the RT3-DIC calculations. The central part of the pantograph contracts more in the transverse direction, ε xx , than its longitudinal expansion, ε yy . This phenomenon is accompanied with more moderate shear ε xy . In Figure 8 the same fields are shown for the 10th loading step. The grey level residuals only slightly degrade in comparison with the previous loading step, see Fig. 8(a). The pattern of the transverse and longitudinal displacement fields is very similar with higher overall levels. The same observation applies for the strain fields. The deformation mode remains unchanged. It is worth observing that all strain levels in the grips remain close to zero, which is to be expected. This is a further validation of the present results.

The last load level prior to damage inception, i.e. first beam failure, is reported in Fig. 9. In that case the grey level residuals, see Figure 9(a), are significantly higher than the previous two cases (see Figures 7(a) and 8(a)). The chosen kinematics is no longer able to properly describe, in a continuous way, the actual kinematics of the central region of the pantograph. For the parts of the pantograph closer to the grips and the grips themselves, the registration quality is significantly better. This result validates the choice of including part of the speckled grips in the analysis. For this last step the highly deformed region has grown toward both ends of the pantograph, which can be understood by the fact that when struts touch each other, the deformation mechanism moves away from these zones. The DIC results reported herein show that such analyses can be run on pantographs. Prior to damage inception, longitudinal nominal strains of the order of 50% and transverse strains as high as -90% were achieved. Such levels could be measured thanks to the sample procedure followed herein, namely, the grips of the pantograph were speckled. In the future, it is likely that speckled pantographs should be analyzed in order to make the DIC analyses easier. It was also shown that a continuous displacement basis was only partly capturing the actual kinematics of the experiment. This calls for more advanced registration techniques that would be consistent of the particular geometry of the pantographs and their specific kinematics.

3.2. Second test. The peculiarity of the observed behavior was the very low force level, wellbelow the sensitivity of the load cell, that was needed to deform the specimen up to very large elongations (of the order of 30 % of the initial specimen length. Figure 10 shows the measured displacement-force response. The displacement is made non-dimensional, with λ varying between 0 and 1, when the applied displacement of the shorter side of the pantographic sheet along the direction of the longer reaches 60 mm. The displacement-force plot is more like a band than a line, whose thickness is about 1 N in the first part (approximately until λ < 0.7). The reaction forces were less than the load cell uncertainty (which is about 2 N) for elongations up to 14 %. This calls for more resolved load cells [START_REF] Turco | Quantitative analysis of deformation mechanism in pantographic substructures: experiments and modeling[END_REF]. The deformed specimens assumed a very peculiar appearance, which which will be referred to as Vietnam long neck shape (see Figs. 5(b) and 11). In these shapes the fibers undergo finite bending between the interconnecting pivots and, ultimately, the first rupture occurs in a fiber which is extremely elongated, close to the corners between the long and short sides of the specimens (Figs. 5(c)).

The pattern of damage and final failure is very similar to those observed when the pivots supply a significant amount of macro-shear strain energy [START_REF] Turco | Fiber rupture in sheared planar pantographic sheets: numerical and experimental evidence[END_REF], which is sometimes more than 60 %-80 % of the total strain energy [START_REF] Turco | Pantographic lattices with nonorthogonal fibres: experiments and their numerical simulations[END_REF]. Therefore it is concluded that the pantographic microstructure determines such a pattern.

Enhanced Piola-Hencky discrete model

Pantographic sheets were modelled [START_REF] Turco | Hencky-type discrete model for pantographic structures: numerical comparison with second gradient continuum models[END_REF], and evaluated for several loading conditions [START_REF] Turco | Large deformations induced in planar pantographic sheets by loads applied on fibers: experimental validation of a discrete Lagrangian model[END_REF][START_REF] Turco | Non-standard coupled extensional and bending bias tests for planar pantographic lattices. Part I: numerical simulations[END_REF][START_REF] Turco | Non-standard coupled extensional and bending bias tests for planar pantographic lattices. Part II: comparison with experimental evidence[END_REF] with a finite dimensional Lagrangian system, consisting of nodes interconnected two by two by extensional springs, and by rotational springs storing energy when two adjacent segments, interconnecting close nodes, change their relative angle. These nodes are called interconnecting nodes as they establish a mechanical interaction between the two families of fibers constituting the pantographic sheet.

The model introduced in Ref. [START_REF] Turco | Hencky-type discrete model for pantographic structures: numerical comparison with second gradient continuum models[END_REF] presents some difficulties if the described experimental evidences are to be predicted. There is a significant amount of strain energy that is stored via bending of the fibers between the interconnecting nodes that cannot be described by the rotational springs [START_REF] Alibert | Extensional elastica in large deformation as γ-limit of a discrete 1D mechanical system[END_REF]. Additional nodes, which are called intermediate nodes, are considered. They only connect two segments of the same fiber, without coupling it with the transverse family of fibers. As a result, an enhanced Piola-Hencky discrete model is obtained.

The enhanced Piola-Hencky discrete model is formulated for improving the description of the strain energy related to the bending mode of each micro-beam. From experiments it is observed that the beam segment of fibers do bend between two nodes and that it is not possible to properly describe the corresponding bending energy solely with an elastic joint in correspondence with the pivot positions. The simplest way to address this issue is to introduce additional elastic joints between those already used. In this way one avoids the introduction of beam elements of varying sophistications [START_REF] Bilotta | Performance of a high-continuity finite element in three-dimensional elasticity[END_REF][START_REF] Greco | An efficient blended mixed B-spline formulation for removing membrane locking in plane curved Kirchhoff rods[END_REF][START_REF] Greco | B-Spline interpolation of Kirchhoff-Love space rods[END_REF][START_REF] Cazzani | Isogeometric analysis of plane curved beams[END_REF] that can add computational complexities.

Figure 12 shows three elastic joints, depicted by P j-1 , P j and P j+1 in the reference position (on the top), and by p j-1 , p j and p j+1 in the current position (on the bottom). The bending strain energy contribution related to the elastic joint P j (distinguished with a black bullet) reads

E b = b f (β) , ( 1
)
where b is the bending stiffness of the elastic joints, and β the angle between the rigid links p j-1p j and p j -p j+1 . The cosine of the angle β, which expresses the bending strain, is computed by means of Carnot's theorem making use of the coordinates of p j-1 , p j and p j+1 . An enhanced representation of the bending strain energy is obtained by further adding elastic joints between those already mentioned, see Fig. 13, where these additional joints are represented as open circles. The bending contributions related to the elastic joints both between P j-1 and P j , and between P j and P j+1 have to be considered. In Fig. 13 elastic joints between P j-1 -P j and P j -P j+1 ) are shown, whose contributions are computed by using Eq. ( 1). In Fig. 13, γ is the angle, in the current configuration, between the rigid links sharing the considered additional elastic joints. The relationship between the assumed bending strain energy and the constitutive parameter b of the elastic joint, the stiffness of the rotational spring, is determined as the solution to a simple structural problem. The single degree-of-freedom problem is illustrated in Fig. 14. A rigid bar with length ℓ is loaded on the right end with a shearing force F. The displacement on the left end is completely constrained whereas the rotation is resisted by the rotational spring with stiffness b. Let us consider three types of bending strain energy E b differentiated by subscripts 1, 2 or 3

E b 1 = 1 2 b 1 ϕ 2 , ( 2
)
E b 2 = 1 2 b 2 ϑ 2 , ϑ = 2 tan ϕ 2 , ( 3
)
E b 3 = b 3 (1 -cos ϕ) . (4)
Equations ( 2) and ( 4) utilize the Lagrangian parameter ϕ to describe the displacements, while that in Eq. ( 3) uses the Lagrangian parameter ϑ, which is linked to ϕ the change of variables reported in Eq. ( 3). The first expression is widely used in structural mechanics for modeling linear elastic behavior known to be applicable for small rotations. For finite rotation problems, the second [START_REF] Eremeyev | Shell-like Structures[END_REF] and third [START_REF] Dell'isola | Large deformations of planar extensible beams and pantographic lattices: Heuristic homogenisation, experimental and numerical examples of equilibrium[END_REF]) forms have been considered. A simple power series expansion shows that the forms adopted for finite rotations reduce to first order expressions for quadratic terms, while they differ with each other in the higher order terms of ϕ. For all cases, the external work W reads W = Fℓ sin ϕ .

(

) 5
The stationarity condition of the potential energy gives the equilibrium equations for each type of assumed strain energy and, consequently, the relationship between the non-dimensional force Fℓ/b i and the Lagrangian parameter ϕ

Fℓ b 1 = ϕ cos ϕ , ( 6
)
Fℓ b 2 = 2 sin ϕ 2 cos ϕ cos 3 ϕ 2 , ( 7
)
Fℓ b 3 = sin ϕ cos ϕ . ( 8
)
Figure 15 reports the plots of the non-dimensional force Fℓ/b i (i = 1, 2 and 3) assuming that the rotation 0 ≤ 2ϕ/π < 1. The three plots show that when large rotations occur, the non-dimensional force can be significantly different for the assumed models. In particular, the second form represents a stiffening behavior with deformation, while the third form shows a gradual softening behavior. However, it is possible to minimize the difference in the results from the three predictions by selecting a (different) most suitable value for the stiffness parameter b i . The expression most appropriate for modeling the bending energy for the present problem is not immediately apparent. It is clear, however, that for the considered pantographic sheets modeled by the enhanced Piola-Hencky discrete model, which considers thousands of rotational degrees of freedom, the rotations can take values that are relatively small, typically close to the boundaries, as well as those that are rather large, typically in the central part of the sheet. For the sake of convenience, and considering that the stiffness parameter b i will be determined in an empirical manner, the third expression is selected in this work. The strain energy is completed by adding the extensional contribution of each micro-beam

E a = 1 2 a (ℓ -ℓ 0 ) 2 , ( 9
)
where ℓ and ℓ 0 are the length of the micro-beam in the current and reference configuration, respectively, and a its axial stiffness.

For the considered pantographic structure, the two types of strain energies (i.e.rotational and extensional) described above are the only quantities necessary to reconstruct the complete equilibrium path, that is to find the set of couples {λ i , q i } in which λ i is a scalar parameter that represents the applied displacement andq i the vector that collects the parameters which completely define the corresponding deformation.

Numerical code details

The above described enhanced Piola-Hencky discrete model of the pantographic sheet was implemented an in-house MatLab code able to reconstruct the complete equilibrium path using a step-by-step strategy (i.e. find an unknown equilibrium point on the basis of the knowledge of the previous known equilibrium point). The energy associated with the pantographic sheet can then be obtained simply by summing all elementary contributions. If the Lagrangian parameters, which describe the actual configuration, are collected in the vector q then the total strain energy becomes E = E(q). The gradient and the Hessian of the strain energy give the structural response vector s(q) and the tangent stiffness matrix K(q), respectively.

For cases without applied forces (as considered in the experiment presented here), the action is the applied displacement on the smaller side of the sheet along the direction of the longer side and the nonlinear system of equations that expresses the equilibrium conditions reads

s(q) = 0 . (10)
Vector q, which collects the Lagrangian parameters used to describe the actual configuration, is rearranged as

q =        u λ ū        , (11)
where the vector u is the collection of the free Lagrangian parameters, and λ ū of those prescribed. The non-dimensional displacement parameter λ then regulates the step-by-step incrementation of the assigned displacements.

The solution to the nonlinear system of equations ( 10) is obtained, in the framework of a stepby-step strategy, by an algorithm based upon Newton's scheme. If the couple (λ 1 , u 1 ) belongs to the equilibrium path, in other words if the associated residual is zero (i.e. its norm is less than a fixed threshold), a nearby equilibrium point (λ 2 , u 2 ) is computed by means of the iterative scheme based on the first-order approximation of the residual defined as

∆u := u 2 -u 1 = -K-1 s (u 1 , (λ 1 + ∆λ) ū) , (12)
where, in order to make the notation compact, ∆λ = λ 2λ 1 and

s (u 1 , (λ 1 + ∆λ) ū) = s               u 1 (λ 1 + ∆λ) ū               . (13)
Furthermore, the stiffness matrix K is the part of the tangent stiffness matrix, computed in u 1 , related to the unknowns collected in u.

The drawback of Newton's scheme, see Eq. ( 12), is related to the singularity of the tangent stiffness matrix. An arc-length strategy, see [START_REF] Riks | An incremental approach to the solution of snapping and buckling problems[END_REF], is implemented. The equilibrium path is described by using as parameter its arc-length. A detailed description of an efficient implementation of the arc-length iterative scheme can be found in [START_REF] Turco | Elasto-plastic analysis of Kirchhoff plates by high simplicity finite elements[END_REF]. Here we only describe necessary computational details of the present code. A simple way to introduce some form of orthogonality conditions is to start from the residual expression

r u 1 + ∆u + u, (λ 1 + ∆λ + λ) ū ≈ s (u 1 + ∆u, (λ 1 + ∆λ) ū) + K        u λ ū        = 0 , (14)
where λ and u define the needed corrections to the extrapolations λ 1 + ∆λ and u 1 + ∆u. As before, K is the part of the tangent stiffness matrix, now computed in u 1 + ∆u, related to the free unknowns. A general form of the constraint equation necessary to make the number of equations equal to the unknowns reads

∆u T C u + γ∆λ λ = 0 , (15)
where C and γ are a matrix and a scalar, respectively, that can be chosen to, e.g., improve the convergence or make simpler the necessary calculations. The condition ( 15) is general and can be interpreted as a representation of orthogonality between extrapolation (∆λ, ∆u) and correction ( λ, u).

From Equation ( 14) the correction u on the displacement are obtained

u = -K-1 s u 1 + ∆u, (λ 1 + ∆λ + λ) ū , (16)
Using Eq. ( 15), the correction λ becomes

λ = ∆u T C K-1 s (u 1 + ∆u, (λ 1 + ∆λ) ū) γ∆λ -∆u T C K-1 ∆s , (17)
where

∆s = s u + ∆u, (λ + ∆λ + λ) ū -s (u + ∆u, (λ + ∆λ) ū) λ . (18)
Equations ( 16) and ( 17) are general but can be particularized to obtain very simple expressions and computational efficiency. Choosing as (γ, C) the couple (0, K) very simple expressions for the corrections are obtained. The corresponding expression for λ is given by λ

= - ∆u T s (u 1 + ∆u, (λ 1 + ∆λ) ū) ∆u T ∆s . (19)
The algorithm outlined above requires the definition of the first extrapolation in each step. The straightforward way is its evaluation using the results of the previous steps. Using the quantities already defined

∆λ = m(λ 1 -λ 0 ) , ∆u = m(u 1 -u 0 ) , (20)
having used the adaptive coefficient m to modify the arc-length during the step-by-step procedure. In practice, the adaptive coefficient serves to increase the step-length in the parts of the equilibrium path that are practically linear, and decrease the step-length in the parts of the equilibrium path that are highly nonlinear. A simple and efficient formula for the adaptive coefficient m, see [START_REF] Clarke | A study of incremental-iterative strategies for non-linear analyses[END_REF], is

m = 1 - r l -n l r l + n l , (21)
which modifies the arc-length on the basis of r l , the actual loops, and n l , the needed loops, necessary to achieve the convergence in the current step. For the first step m = 1 and the value of ∆λ is fixed. In this way, the arc-length is implicitly defined by computing the corresponding value of ∆u.

Calibration of the numerical model

The determination of the constitutive parameters introduced in the enhanced Piola-Hencky model starting from the geometry and the mechanical properties of the pantographic specimens is not straightforward [START_REF] Abali | An energy-based method to determine material constants in nonlinear rheology with applications[END_REF]. The issue concerning micro-macro identification processes has been addressed from many points of view [START_REF] Alibert | Second-gradient continua as homogenized limit of pantographic microstructured plates: a rigorous proof[END_REF] and we refrain from studying it here. Instead the introduced micro-constitutive parameters are calibrated in order to fit the experimental data at the macro-level. To this end a parametric study was performed wherein some simple estimates of the stiffness parameters are obtained, and then varied over a range suggested by these first estimates. Given the DIC results reported in Section 3.1, another could be followed in the future, provided the measurements are performed at the micro-level.

The procedure starts by considering the geometrical data of the technical drawing reported in Fig. 2. From these a simple first estimate of the stiffness parameters of the discrete model is obtained using the results given by the solution of De Saint Venant problem. In practice, for a straight beam of initial length L, cross-sectional area A, and inertia J, the stiffnesses a and b are estimated as

a = YA ℓ , b =c Y J ℓ , (22)
where Y is the Young's modulus=, L/ℓ the number of segments whose sum is the length of the beam, and c = 1ℓ/L a coefficient accounting for the actual number of bending elastic joints. This coefficients rapidly goes to one when the number of segments increases (this hypothesis is not made herein). Geometric data taken from the technical drawing reported in Fig. 2 may be considered accurate of within 0.1 mm. Conversely, the value of the Young's modulus, which affects both stiffnesses a and b, has to be evaluated carefully. Technical references on the polyamide used in the 3D printing process suggest a value of 1700 MPa. 2 As another estimate that can provide an improved fit, the results of the numerical simulation are rescaled simultaneously modifying the stiffness parameters a and b with an additional condition, for example that R(λ = 1) be the same for the computed and measured global reactions. The Young's modulus Y was determined minimizing the discrepancy, through a suitable measure of chosen quantities, between numerical and physical models. A leasts squares approach has been selected [START_REF] Turco | Tools for the numerical solution of inverse problems in structural mechanics: review and research perspectives[END_REF] δ

= x -x 2 , ( 23
)
where δ is the measure of the discrepancy, in this case the L 2 -norm of the differences between x and x, that is between the predicted, by the numerical model, and the measured, from the physical experiment, data. There are many possible choices (e.g.the discrepancy between the configuration of the pivots relative to the same value of the non-dimensional applied displacement λ, or the discrepancy between the predicted and measured global reaction(s) for various levels of λ).

The best identified values, using as additional condition that R(λ = 1) be the same for the numerical simulation and physical experiment, for the stiffnesses are a = 307.1 N/mm and b = 62.23 Nmm, which correspond to an identified value of Y = 950 MPa. The numerical results relative to this set of stiffnesses a and b are reported in Figs. 17 and 18. In Fig. 17 the history of deformation is reported for four selected valued of the non-dimensional applied displacement λ (0.25, 0.5, 0.75 and 1) whereas in Fig. 18 the comparison between the computed and measured reactions versus λ (a) and the strain energy split in axial and bending parts are shown when λ increases (b).

The most important differences between the numerical simulation and the experiment are:

(1) the pronounced shrinkage in the middle part of the specimen mostly for high values of λ; (2) the experimental curve R(λ) may be approximatively considered as bilinear formed by two straight lines; the same could be said for the computed curve with the transition point shifted forward and a steeper second segment.

Guided by the estimates described above a parametric study of the stiffness set was performed by computing the global reaction R when λ varies between 0 and 1. For the sake of clarity, only the curves obtained by varying the axial stiffness a from 50 to 500 N/mm and the bending stiffness b from 20 Nmm to 200 Nmm are reported (Fig. 19).

From the previous results, it is concluded that:

(1) a reduction of the slope of the second part of the R(λ) curve is obtained increasing b;

(2) the value R(λ = 1) is largely influenced by a and much less from b;

(3) the slope of the first part of R(λ) curve is constant and the curves are shifted upward increasing b and much less increasing a.

These observations suggest a further refinement on the stiffnesses a and b, using as best fit criterion that of minimizing the discrepancy between the whole curve, keeping the value R(λ = 1), gives the values 190 N/mm and 180 Nmm, respectively. Figures 20 and21 show the results obtained by using this set of stiffnesses. Figure 20 reports the history of the deformation for λ varying in the set {0.25, 0.5, 0.75, 1} (color map shows the density of the energy level on the micro-beam). Figure 21 compares the computed and measured R(λ) responses, and the strain energy contributions, split into extensional E a and flexural E b parts, respectively. This last fitting procedure seems rather reliable. The coefficients that were obtained by fitting the shapes and force/displacement relationship for selected elongations continue to be effective in predicting the experimental results outside the range of fitting as well. The procedure could be further refined by introducing some shear resistance in the pivots, thereby introducing some spatial variability.

Concluding remarks and future challenges

A discrete Lagrangian model of the considered mechanical system was introduced based on the physical understanding of its properties. No continuum model was considered or constructed as the result of a homogenization step. The most suitable continuum model should include the dependence of strain energy on the second gradient of displacement, and this requires the introduction of more complex modeling concepts and constitutive parameters. The lack of micro-scale torsional (and macro-scale shear) strain energy in the interconnecting pivots is the ultimate cause of the Vietnam long neck shapes that are observed both in experiments and numerical simulations. These shapes are also the source of many conceptual and numerical difficulties whose rigorous mathematical treatment requires, when using continuum models [START_REF] Giorgio | Continuum modelling of pantographic sheets for out-of-plane bifurcation and vibrational analysis[END_REF][START_REF] Bertram | Gradient materials with internal constraints[END_REF][START_REF] Scerrato | Three-dimensional instabilities of pantographic sheets with parabolic lattices: numerical investigations[END_REF], the use of advanced functional concepts [START_REF] Eremeyev | Linear pantographic sheets: Existence and uniqueness of weak solutions[END_REF]. There is another reason for which a discrete model was preferred herein. There are many material systems in which the discrete approach is preferable over a continuum formulation. The presented discrete model and the second gradient continuum Ansatz are not able to describe the effects of compression when fibers actually touch each other in very large deformation regimes [START_REF] Andreaus | Experimental and numerical investigations of the responses of a cantilever beam possibly contacting a deformable and dissipative obstacle under harmonic excitation[END_REF]. In order to handle such phenomenology some strongly nonlinear nodal interactions need to be incorporated, which in discrete models make energetically difficult the too close placement of adjacent fibers, or exponentially increasing deformation potentials in the continuum model. Finally, some future challenges, which will be tackled in the next future, are listed:

(1) the numerical model developed herein considered the elastic regime. Repeated loading/unloading steps should be analyzed and extra constitutive parameters will be required [START_REF] Altenbach | On the constitutive equations of viscoelastic micropolar plates and shells of differential type[END_REF]; (2) the experimental evidence compared with numerical predictions involved very few data (e.g. reaction forces, some global elongations); it has been observed a posteriori that the deformed and predicted shapes are rather close. This methodology, however reliable and effective, does not allow for a complete use of all available data. It was shown that digital image correlation could be used to quantify the whole experimental kinematics. A complete comparison between measured and computed displacement fields is to be performed for validation and calibration purposes. (3) the experimental set up has high performances. However, its elongation is insufficient to always achieve failure. Changing the geometry of the pantographic substructure (e.g. the fiber angle [START_REF] Turco | Pantographic lattices with nonorthogonal fibres: experiments and their numerical simulations[END_REF]) may lead to failure of specimens in the available range of applied displacements, say, in shear tests [START_REF] Turco | Fiber rupture in sheared planar pantographic sheets: numerical and experimental evidence[END_REF]. Similarly the beam spacing or their sections can be altered; (4) in this paper the stiffnesses of the enhanced Piola-Hencky model were assumed to be invariable over the whole pantographic sheet. This hypothesis may be relaxed with no increasing computational cost. That allows the study of pantographic sheets with defects, see [START_REF] Turco | Pantographic structures presenting statistically distributed defects: numerical investigations of the effects on deformation fields[END_REF], or the optimization of the stiffnesses of micro-beams and pivots; (5) the extension of this study to dynamic cases could give some suggestions for the use of pantographic sheets for mitigating the effects of shocks and vibrations [START_REF] Vescovo | Dynamic problems for metamaterials: review of existing models and ideas for further research[END_REF][START_REF] Engelbrecht | Reflections on mathematical models of deformation waves in elastic microstructured solids[END_REF][START_REF] Berezovski | Interfaces in micromorphic materials: wave transmission and reflection with numerical simulations[END_REF][START_REF] Giorgio | Isola. Dynamics of 1D nonlinear pantographic continua[END_REF];

Figure 3 .

 3 Figure 3. 3D rendering of cells of a pantographic sheet with pivots without strain energy (courtesy of M. Golaszewski).

Figure 4 .

 4 Figure 4. Micro-beam of a 3D printed polyamide pantographic sheet.

Figure 5 .

 5 Figure 5. Grey level images of the pantograph in the reference configuration (a), last analyzed deformed configuration (b) and at failure (c).

Figure 6 .

 6 Figure 6. Region of interest (red box) analyzed via RT3-DIC (a) and finite element mesh with 25-pixel T3 elements (b).

Figure 7 .

 7 Figure 7. Correlation results for the 5th loading step: grey level residuals (a), u x (b) and u y (c) displacement fields (expressed in pixels); ε xx (d), ε xy (e), ε yy (f) nominal strain fields.

Figure 8 .Figure 9 .

 89 Figure 8. Correlation results for the 10th loading step: grey level residuals (a), u x (b) and u y (c) displacement fields (expressed in pixels); ε xx (d), ε xy (e), ε yy (f) nominal strain fields.

Figure 10 .

 10 Figure 10. Extension bias test. Measured force R vs. nondimensional displacement λ.

(a) λ = 0 ( 1 Figure 11 .

 0111 Figure 11. Extension bias test. Deformed shape for five values (0, 0.25, 0.5, 0.75 and 1) of the non-dimensional displacement λ.

Figure 12 .

 12 Figure 12. Kinematics of a discrete micro-beam. Reference (on the top-left) and current (on the bottom-right) configurations.

Figure 13 .

 13 Figure 13. Kinematics of the enhanced approximation of a micro-beam. Reference (on the top-left) and current (on the bottom-right) configurations.

Figure 14 .

 14 Figure 14. Cantilever beam loaded with a shearing force as a simple one degree of freedom problem.

1 E b 2 E b 3 Figure 15 .

 12315 Figure 15. Non-dimensional force Fℓ/b i versus rotation ϕ/ϕ max .

Figure 16 .

 16 Figure 16. Extension bias test. Global reaction R computed (in red) and measured (in black) vs. non-dimensional applied displacement λ for Y = 1700 MPa.

5 (d) λ = 1 Figure 17 .

 5117 Figure 17. Extension bias test. Deformation history for four values (0.25, 0.5, 0.75 and 1) of the non-dimensional applied displacement λ for a = 307.1 N/mm and b = 62.23 Nmm (the color map shows the strain energy density on each micro-beam).

Figure 18 .

 18 Figure 18. Extension bias test. Global reaction (in black: measured force) R (a) and contributions to the strain energy E of the axial E a and bending E b parts (b) vs. non-dimensional applied displacement λ when a = 307.1 N/mm and b = 62.23 Nmm.

Figure 19 .

 19 Figure 19. Global reaction R vs. non-dimensional applied displacement λ varying the stiffnesses a and b.

= 1 Figure 20 .

 120 Figure 20. Extension bias test. Deformation history for four values (0.25, 0.5, 0.75 and 1) of the non-dimensional applied displacement λ for a = 190 N/mm and b = 180 Nmm (the color map shows the strain energy density on each microbeam).

Figure 21 .

 21 Figure 21. Extension bias test. Global reactions (in black: measured force) R (a) and contributions to the strain energy E of the axial E a and bending E b parts (b) vs. non-dimensional applied displacement λ when a = 190 N/mm and b = 180 Nmm.

Figure 23 .

 23 Figure 23. Numerical simulation of the bias extension test. Global reaction (in black: measured force) R (a) and contributions to the strain energy E of the axial E a and bending E b parts (b) vs. the non-dimensional applied displacement λ for a = 148 N/mm and b = 280 Nmm .

see https://www.shapeways.com/rrstatic/material docs/mds-strongflex.pdf.

(6) an adaptation of the discrete model presented herein could be profitably used to model wrinkling mechanics of specifically designed specimens [START_REF] Harrison | Towards comprehensive characterisation and modelling of the forming and wrinkling mechanics of engineering fabrics[END_REF]; [START_REF] Hild | Digital image correlation. Optical methods for solid mechanics: a fullfield approach[END_REF] the enhanced Piola-Hencky discrete model could be considered in the study granular media [START_REF] Challamel | Stability of non-conservative elastic structures under additional kinematics constraints[END_REF][START_REF] Misra | Granular micromechanics based micromorphic model predicts frequency band gaps[END_REF][START_REF] Misra | Grain-and macro-scale kinematics for granular micromechanics based small deformation micromorphic continuum model[END_REF][START_REF] Misra | Elastic behavior of 2D grain packing modeled as micromorphic media based on granular micromechanics[END_REF]], or foams with Pott's model [START_REF] De Masi | Potts models in the continuum uniqueness and exponential decay in the restricted ensembles[END_REF][START_REF] De Masi | Coexistence of ordered and disordered phases in Potts models in the continuum[END_REF][START_REF] Grimmett | Correlation inequalities for the Potts model[END_REF].