parallel to one of the two coordinate axes. Each beam in one family is connected by means of an internal pivot to all the beams belonging to the other family and intersecting it (see Fig. 1). Pantographic lattices considered here are slightly different from those considered in [START_REF] Alibert | Truss modular beams with deformation energy depending on higher displacement gradients[END_REF]: indeed there, the topological connection between two arrays of pantographic beams was obtained superimposing the central pantographic pivots (which we denoted with the letter B in the following Fig. 1). The comparison between the pantographic lattice studied here and that which has been previously studied will be the object of future work.

Unfortunately, the pioneering paper by [START_REF] Maxwell | On the calculation of the equilibrium and stiffness of frames[END_REF] has limited the consideration of engineers to beam lattices statically (over-)determined and including only clamping constraints, or to simple trusses. Actually, Maxwell only wanted to supply to the designers a useful mathematical and computational tool and limited his attention to a class of structures simply because their "calculation" via analytical methods was viable. One cannot attribute to Maxwell any intention to limit the class of structures of interest but simply the wish to allow for a rigorous design based on solid mathematical grounds in cases closed form calculations are possible, in an époque when this was the only possible form of computation. Most likely beyond Maxwell's intentions, his algorithmic procedure imposed a paradigm in engineering design: only truss structures or structures with clamped internal constraints between beams are commonly considered in engineering design. This paradigm has no logical nor (by now) "practical" ground: it was actually imposed by the limitations induced by the demand of studying the deformation problems using: (i) the analytical closed form solutions and (ii) as calculating devices simply with paper and pencil. However, this design paradigm has been incorporated (without any logical or algorithmic reason) in the most modern tools of numerical calculation: in many of them it is unreasonably difficult to impose internal constraints between beams, except when this constraint is a clamp. As a consequence, the attention of those researchers looking for homogenized field models (and demonstrating their results using heuristic -as done in the present paper -or by more rigorous identification methods -see [START_REF] Alibert | Second-gradient continua as homogenized limit of pantographic microstructured plates: a rigorous proof[END_REF] for a rigorous result in a similar context) has been limited to the "traditional" class of beam lattices, which has indeed a natural continuum counterpart, namely the standard Cauchy first gradient model. In the beautiful book by [START_REF] Kuznetsov | Underconstrained structural systems[END_REF], underconstrained structural systems are courageously studied against the prescriptions of the aforementioned paradigm: the deep analysis of the author leads to a clear understanding of the "practical" importance of systems which have been unreasonably neglected in engineering sciences. Indeed, not considering the possibility of their use has considerably reduced the potential performances of structural systems: pantographic structures we consider here are underconstrained and show very interesting mechanical features, among which we remark the fact that their continuum counterpart is given by second gradient models.

Classical homogenization of such structures with a discrete topology towards an equivalent Cauchy continuum has been an active research field since many years, and several methods have been developed over the last two decades to find an effective continuum [START_REF] Bornet | Homogénéisation en mécanique des matériaux 1[END_REF][START_REF] Chesnais | Effects of the local resonance on the wave propagation in periodic frame structures: generalized Newtonian mechanics[END_REF][START_REF] Reis | Homogénéisation automatique de milieux discrets périodiques[END_REF][START_REF] Reis | Modélisation de matériaux continus par un treillis répétitive en vue de la simulation de certains cas de plasticité[END_REF][START_REF] Sanchez-Hubert | Introduction aux méthodes asymptotiques et à l'homogénéisation -application à la mécanique des milieux continus[END_REF]) amongst these works, many of them treat the homogenized behavior of elastic and plastic materials reinforced by rigid fibers [START_REF] Bellieud | Homogenization of a soft elastic material reinforced by fibers[END_REF][START_REF] Contrafatto | An enriched finite element for crack opening and rebar slip in reinforced concrete members[END_REF][START_REF] El Jarroudi | Asymptotic behaviour of a cylindrical elastic structure periodically reinforced along identical fibres[END_REF][START_REF] Luongo | Mode localization in dynamics and buckling of linear imperfect continuous structures[END_REF][START_REF] Sili | Homogenization of an elastic medium reinforced by anisotropic fibers[END_REF][START_REF] Turco | A strategy to identify exciting forces acting on structures[END_REF], also in the context of the mechanical study of living and growing tissues (see [START_REF] Federico | An energetic approach to the analysis of anisotropic hyperelastic materials[END_REF]Grillo et al., 2009a;Grillo, Wittum, Giaquinta, & Mi ćunovi ć, 2009b;[START_REF] Lekszycki | Modelling of bone adaptation based on an optimal response hypothesis[END_REF]. The classical homogenization techniques nevertheless encounter limitations when the wavelength of the loading or deformation field becomes comparable with the typical microstructure size; it is a well established fact in the literature [START_REF] Buechner | Size effects in the elasticity and viscoelasticity of bone[END_REF][START_REF] Harrigan | Limitations of the continuum assumption in cancellous bone[END_REF]) that the classical Cauchy theory does not allow for the correct prediction of the mechanical response of bone at sufficiently small scale levels and/or when specific loading conditions are adopted (Andreaus, Giorgio, & Lekszycki, 2014b;Andreaus, Giorgio, & Madeo, 2014c), requiring the improvement of these theories by incorporating additional intrinsic parameters and internal length scales to correlate the microstructure with the macrostructure. Indeed, first gradient Cauchy theories need to be improved by incorporating additional intrinsic parameters and internal length scales to correlate the microstructure with its macro-modeling via suitably introduced macro-fields (see e.g. [START_REF] Auffray | Analytical continuum mechanics a la Hamilton-Piola: least action principle for second gradient continua and capillary fluids[END_REF]dell'Isola, Andreaus, & Placidi, 2015a): an important field of investigation in this context concerns the introduction of second gradient models in the description of damage, see [START_REF] Aifantis | On the role of gradients in the localization of deformation and fracture[END_REF] for a general perspective and [START_REF] Misra | Granular micromechanics model for damage and plasticity of cementitious materials based upon thermomechanics[END_REF], [START_REF] Misra | Thermomechanics-based nonlinear rate-dependent coupled damage-plasticity granular micromechanics model[END_REF], Placidi (2015a), [START_REF] Placidi | Gedanken experiments for the determination of two-dimensional linear second gradient elasticity coefficients[END_REF], [START_REF] Rinaldi | A rational model for 2D disordered lattices under uniaxial loading[END_REF], [START_REF] Rinaldi | A rational model for 2D disordered lattices under uniaxial loading[END_REF]Placidi (2014), Scerrato, Giorgio, Della Corte, Madeo, and[START_REF] Scerrato | A micro-structural model for dissipation phenomena in the concrete[END_REF] for more recent results in this direction. Moreover, if one replaces the pivots in the nodes of the lattice so as to have contact interaction between the two families of beams, it is clear that a damaged pantographic structure could also be modeled by means of piecewise linear systems (Andreaus, Chiaia, & Placidi, 2013b). It has also to be pointed out that a homogenized macroscopic model of a pantographic structure is intrinsically anisotropic, with the anisotropic behavior depending on the distribution of orientation of the fibers, in the same way as it occurs in polycrystalline materials [START_REF] Placidi | Thermodynamics of polycrystalline materials treated by the theory of mixtures with continuous diversity[END_REF][START_REF] Placidi | On the role of grain growth, recrystallization and polygonization in a continuum theory for anisotropic ice sheets[END_REF][START_REF] Seddik | Application of a continuum-mechanical model for the flow of anisotropic polar ice to the EDML core, Antarctica[END_REF].

The objective of homogenization towards generalized continua considering either additional degrees of freedom (like the Cosserat medium as advocated by the two Cosserat brothers in 1986 and[START_REF] Cosserat | At the origins and in the vanguard of peri-dynamics, non-local and higher gradient continuum mechanics. An underestimated and still topical contribution of Gabrio Piola[END_REF]the micromorphic medium;[START_REF] Cosserat | La théorie de l'élasticité[END_REF][START_REF] Cosserat | At the origins and in the vanguard of peri-dynamics, non-local and higher gradient continuum mechanics. An underestimated and still topical contribution of Gabrio Piola[END_REF][START_REF] Eringen | Nonlinear theory of simple microelastic solids[END_REF] or additional higher order gradients (like the second order gradient continuum; [START_REF] Mindlin | On first strain gradient theories in linear elasticity[END_REF][START_REF] Toupin | Elastic materials with couple stresses[END_REF] is to remedy to these limitations and to extend the range of validity of the continuum approach beyond the strict assumption of the scale separation, short range micro-interactions or absence of high contrast of physical and geometrical properties at micro-level (see, [START_REF] Camar-Eddine | Determination of the closure of the set of elasticity functionals[END_REF][START_REF] Pideri | A second gradient material resulting from the homogenization of a heterogeneous linear elastic medium[END_REF][START_REF] Trinh | The role of the fluctuation field in higher order homogenization[END_REF][START_REF] Trinh | Evaluation of generalized continuum substitution models for heterogeneous materials[END_REF]. As mentioned in (Forest, 2006), the method to be used must be able to account for the effect of the morphology and distribution of phases on the material response and to predict scale effects (size effects of the constituents and impact of strong gradients of the loading); several methods have been developed towards this objective in the literature (see [START_REF] Alibert | Truss modular beams with deformation energy depending on higher displacement gradients[END_REF]Dos Reis & Ganghoffer, 2012a;[START_REF] Forest | Mechanics of generalized continua: construction by homogenization[END_REF][START_REF] Forest | Homogenization methods and the mechanics of generalized continua -part 2[END_REF][START_REF] Hirschberger | Computational multiscale modelling of heterogeneous material layers[END_REF][START_REF] Seppecher | Linear elastic trusses leading to continua with exotic mechanical interactions[END_REF].

The Cosserat continuum presents several limitations to the study of the localization of deformations in presence of shear gradients, and a second order grade continuum is more suitable for this purpose (for the asymptotic methods needed in this context see e.g. [START_REF] Luongo | Mode localization by structural imperfections in one-dimensional continuous systems[END_REF][START_REF] Luongo | Mode localization in dynamics and buckling of linear imperfect continuous structures[END_REF]. Several authors have turned their attention towards the development of rigorous homogenization procedures showing the ability to account for the heterogeneous nature of the material at the micro-level using a second gradient macroscopic constitutive law [START_REF] Camar-Eddine | Determination of the closure of the set of elasticity functionals[END_REF][START_REF] Pideri | A second gradient material resulting from the homogenization of a heterogeneous linear elastic medium[END_REF]; works in this direction consider both linear elastic materials [START_REF] Mindlin | Second gradient of strain and surface tension in linear elasticity[END_REF] and non-linear materials [START_REF] Bardenhagen | Derivation of higher order gradient continuum theories in 2, 3-D nonlinear elasticity from periodic lattice models[END_REF][START_REF] Mareno | Global continuation in second-gradient nonlinear elasticity[END_REF]. Higher order homogenization schemes have been built for architectured materials in both the linear and nonlinear regimes in the work of Trinh (2011). [START_REF] Auffray | Comportement des matériaux cellulaires; élaboration, caractérisation et modélisation prédictive des propriétés Ph[END_REF] focused on the mechanical response of honeycombs, accounting for scale effects being present. The recent work of El Jarroudi (2013) treats the homogenization of nonlinear elastic material in contact with a set of more rigid nonlinear elastic fibers periodically distributed within the structure. In these papers (see also [START_REF] Alibert | Truss modular beams with deformation energy depending on higher displacement gradients[END_REF], it is proven that some specific, highly organized, microstructures show a macroscopic behavior which, at macro-level, cannot be described with first gradient models. Unfortunately, it is not still available a general method allowing for the algorithmic classification of a given micro-structure determining the most appropriate macro model suitable for its description.

In the present work, we shall construct an equivalent second order gradient continuum for the pantograph structure considered in the linearized regime and in the neighborhood of a reference configuration where the beams are straight. These are exactly the assumptions accepted in [START_REF] Placidi | A variational approach for a nonlinear one-dimensional damage-elasto-plastic second-gradient continuum model[END_REF], where the aim is complementary to the one of the present paper: indeed, while here we aim to identify the macroscopic constitutive equations via the specification of the microstructure of the system, in the aforementioned paper, this identification is planned via the experimental evidence obtained in a series of well-designed tests.

In future developments, we shall address the problem of finding the equivalent continuum for pantographic lattices in the regime of large deformations. We explicitly remark here that the 2D continuum which we determine is the most suitable model for pantographic lattices, and is more general than standard plate models, as it involves second gradient of in-plane displacements: these 2D continua were already investigated in the literature (see e.g. [START_REF] Altenbach | Linear theory of shells taking into account surface stresses[END_REF][START_REF] Pagnini | Reliability analysis of wind-excited structures[END_REF][START_REF] Kuznetsov | Underconstrained structural systems[END_REF]dell'Isola & Steigmann, 2015;[START_REF] Giorgio | Pattern formation in the three-dimensional deformations of fibered sheets[END_REF][START_REF] Pietraszkiewicz | Extended non-linear relations of elastic shells undergoing phase transitions[END_REF]Steigmann & dell'Isola, 2015) and they present very interesting mathematical features.

Thereby, we shall extend the first order homogenization schemes recently developed for the determination of the effective mechanical properties of periodical lattices considered as Cauchy or micropolar continua (Dos Reis & Ganghoffer, 2012a) towards second order gradient continua. Following the ideas presented in recent contributions devoted to the 1D case [START_REF] Carcaterra | Macroscopic description of microscopically strongly inhomogenous systems: A mathematical basis for the synthesis of higher gradients metamaterials[END_REF][START_REF] Elnady | Construction of second order gradient continuous media by the discrete asymptotic homogenization method[END_REF], we shall deal with 2D structures in the present work. The developed framework based on the principle of virtual powers allows incorporating the local microstructural effects via the consideration of the second order displacement gradient.

The similarities of the methods we present here and the asymptotic expansion and perturbative methods used, in slightly different context, by (see e.g. [START_REF] Luongo | Analytical and numerical approaches to nonlinear galloping of internally resonant suspended cables[END_REF][START_REF] Rinaldi | A rational model for 2D disordered lattices under uniaxial loading[END_REF] are remarkable: we expect that multiscale methods will need to be used explicitly when lattices with multiple different length scales will be considered, thereby leading to homogenized higher gradient continua. In this regards, it has to be noticed that higher gradient continua require enhanced integration schemes to avoid some numerical problems as locking or spurious modes (see e.g. [START_REF] Cazzani | On the dynamics of a beam partially supported by an elastic foundation: an exact solution-set[END_REF][START_REF] Cazzani | On some mixed finite element methods for plane membrane problems[END_REF][START_REF] Cazzani | Numerical aspects of coupling strongly frequency-dependent soil-foundation models with structural finite elements in the timedomain[END_REF]Cazzani, Malagù, & Turco, 2014a;[START_REF] Garusi | An unsymmetric stress formulation for Reissner-Mindlin plates: a simple and locking-free rectangular element[END_REF][START_REF] Reccia | FEM-DEM modeling for out-of-plane loaded Masonry panels: a limit analysis approach[END_REF]: in the present instance, we use some packages from COMSOL adapted by introducing some fictitious microstructural tensor fields equated to displacement gradients via Lagrange multipliers. This approach can be safely made more effective by introducing isogeometric integration schemes (see e.g. [START_REF] Greco | B-Spline interpolation of Kirchhoff-Love space rods[END_REF], 2014;[START_REF] Greco | A procedure for the static analysis of cable structures following elastic catenary theory[END_REF]Cazzani et al., 2014aCazzani et al., , 2014b[START_REF] Cazzani | Constitutive models for strongly curved beams in the frame of isogeometric analysis[END_REF][START_REF] Turco | A three-dimensional B-spline boundary element[END_REF]. It has to be remarked that the most peculiar feature of the obtained homogenized energy consists in the existence of so called floppy modes, e.g. displacements which are not rigid and which correspond to vanishing deformation energy. Indeed, it is easy to check that the found homogenized energy admits as floppy modes all first order polynomials with respect to the Lagrangian coordinates. Also very peculiar is its dynamics behavior: it seems suggestive that pantographic structures, once suitably modified by the introduction of a set of microscopic oscillators, for instance by considering generalized beams as those studied in [START_REF] Piccardo | A complete dynamic approach to the Generalized Beam Theory cross-section analysis including extension and shear modes[END_REF], may favor the energy trapping processes described in [START_REF] Alibert | Truss modular beams with deformation energy depending on higher displacement gradients[END_REF], [START_REF] Carcaterra | Theoretical foundation of apparent damping and irreversible energy exchange in linear conservative dynamical systems[END_REF] and [START_REF] Koc | Energy sinks: vibration absorption by an optimal set of un-damped oscillators[END_REF]. The final aim should be the design of novel metamaterials (see e.g. [START_REF] Boutin | Homogenisation of periodic discrete medium: application to dynamics of framed structures[END_REF][START_REF] Chesnais | Effects of the local resonance on the wave propagation in periodic frame structures: generalized Newtonian mechanics[END_REF][START_REF] Placidi | Reflection and transmission of plane waves at surfaces carrying material properties and embedded in secondgradient materials[END_REF].

The paper is organized as follows: the effective strain energy density per unit surface of the pantograph is derived based on an analytical method in Section 2. As a main aspect, a general methodology is presented in Section 3 to compute the effective mechanical properties of beam lattices up to the second order gradient of the displacement. We shall expose into details (Section 3) the technical steps required to construct the equivalent second gradient continuum, thereby highlighting the forces dual to the first and second order kinematic variables in the sense of the virtual power of internal forces. The proposed method is a variant of similar homogenization schemes already developed for second gradient continuum [START_REF] Kouznetsova | Multi-scale constitutive modelling of heterogeneous materials with a gradient-enhanced computational homogenization scheme[END_REF], with the main difference that the topology of the initial medium is discrete. In Section 4, numerical simulations are performed, exemplifying some features of the obtained second gradient continuum. We conclude this work by a summary of the main results and a few perspectives (Section 5).

Analytical derivation of the effective energy of the pantograph

The pantograph is well-known, for it is a typical example of a structure exhibiting no first order elasticity, but instead its elasticity appears directly at second order when it is viewed as a continuum, thus it is a second order grade continuum. The goal of the work is to derive the effective continuum behavior of the pantograph through its continuum energy.

The unit cell of the pantograph is represented on Fig. 1, together with the kinematic angular variables. The kinematic assumption is as follows: both beams, linking together the pair of nodes (i -1, j -1) -(i + 1, j + 1) and (i + 1, j -1) -(i -1, j + 1), are inextensible and articulated by a pivot at the central node (i, j). This unit is repeated by periodicity in the plane in order to replace the discrete structure by an effective continuum body with homogenized properties at a mesoscopic scale of description. The displacement and rotation are continuous at all nodes.

The deformation energy of the pantograph structure viewed as a continuum is evaluated based on two different approaches: the first method is analytical and relies on the work of [START_REF] Askar | A structural model of a micropolar continuum[END_REF]. The second method is more rigorous and relies on the discrete asymptotic homogenization technique to build the effective second order gradient continuum based on asymptotic expansions.

Analytical expression of the energy

The general expression of the internal deformation energy associated to a beam element (ij, i+1 -j+1), Fig. 2, is obtained in a straightforward manner as:

W = 1 2 ES L ξ j+1 i+1 -ξ j i 2 + 1 2 EI L θ j+1 i+1 -θ j i 2 + 1 2 12EI L 3 η j+1 i+1 -η j i - L 2 θ j+1 i+1 + θ j i 2
This energy is made of three terms representing the different forms of strain energy, namely extensional, flexural and shear contributions. The beams of the pantograph are supposed to be inextensible: the linearized strong from of the inextensibility condition writes for a beam element:

l 2 o = l 2 , with l o the initial beam length, that is r i j -r i+1 j+1 2 = Cte = l 2 o
The following relation between the Lagrangian and Eulerian positions holds

r = R + u and R extremity -R origin = l o e b
The variable u denotes the displacement. Thus, it holds that

r j j -r i+1 j+1 2 = R i j + u i j -R i+1 j+1 -u i+1 j+1 2 = l 2 o e b ⇒ -l o e b + ξ i j -ξ i+1 j+1 e b + η i j -η i+1 j+1 e b⊥ 2 = l 2 o e b
Previous condition becomes in the small perturbations framework

l 2 o e b -2l o ξ i j -ξ i+1 j+1 = l 2 o e b
Thus, one obtains the following linearized inextensibility condition:

ξ i j -ξ i+1 j+1 = 0
Consequently, the energy expresses as

W = 1 2 EI L θ i+1 j+1 -θ i j 2 + 1 2 12EI L 3 η i+1 j+1 -η i j - L 2 θ i+1 j+1 + θ i j 2 introducing therein the flexural rigidity K f = EI L b
, the shear rigidity

K c = 12EI L 3 b
and L b = a the beam length.

Both beams on Fig. 1, linking together the pair of nodes (i -1, j -1) -(i + 1, j + 1) and(i + 1, j -1) -(i -1, j + 1) have shear rigidities denoted K c 1 and K c 2 , and flexural rigidities denoted K f 1 and K f 2 ; these two beams are articulated by a pivot at the central node.

The strain energy of the structure is written as follows:

W = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 1 2 K f 1 θ 1 i+1 j+1 -θ 1 i j + 2 + θ 1 i j - -θ 1 i-1 j-1 2 + 1 2 K f 2 θ 2 i+1 j-1 -θ 2 i j + 2 + θ 2 i j - -θ 2 i-1 j+1 2 + 1 2 K c 1 η 1 i+1 j+1 -η 1 i j -a 2 θ 1 i+1 j+1 + θ 1 i j + 2 + η 1 i j -η 1 i-1 j-1 -a 2 θ 1 i j - + θ 1 i-1 j-1 2 + 1 2 K c 2 η 2 i+1 j-1 -η 2 i j -a 2 θ 2 i+1 j-1 + θ 2 i j + 2 + η 2 i j -η 2 i-1 j+1 -a 2 θ 2 i j - + θ 2 i-1 j+1 2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (1)
Here and in the sequel, indices 1 and 2 refer respectively to the first and second beam.

One denotes separately both contributions to the strain energy associated to each beam:

W 1 = ⎡ ⎢ ⎢ ⎣ 1 2 K f 1 θ 1 i+1 j+1 -θ 1 i j + 2 + θ 1 i j - -θ 1 i-1 j-1 2 + 1 2 K c 1 η 1 i+1 j+1 -η 1 i j -a 2 θ 1 i+1 j+1 + θ 1 i j + 2 + η 1 i j -η 1 i-1 j-1 -a 2 θ 1 i j - + θ 1 i-1 j-1 2 ⎤ ⎥ ⎥ ⎦ (2)
For the second beam, it holds similarly

W 2 = ⎡ ⎢ ⎢ ⎣ 1 2 K f 2 θ 2 i+1 j-1 -θ 2 i j + 2 + θ 2 i j - -θ 2 i-1 j+1 2 + 1 2 K c 2 η 2 i+1 j-1 -η 2 i j -a 2 θ 2 i+1 j-1 + θ 2 i j + 2 + η 2 i j -η 2 i-1 j+1 -a 2 θ 2 i j - + θ 2 i-1 j+1 2 ⎤ ⎥ ⎥ ⎦ (3)
We next elaborate the development of each contribution to the total energy:

W 1 = ⎡ ⎢ ⎢ ⎣ 1 2 K f 1 θ 1 i+1 j+1 -θ 1 i j + 2 + θ 1 i j - -θ 1 i-1 j-1 2 + 1 2 K c 1 η 1 i+1 j+1 -η 1 i j -a 2 θ 1 i+1 j+1 + θ 1 i j + 2 + η 1 i j -η 1 i-1 j-1 -a 2 θ 1 i j - + θ 1 i-1 j-1 2 ⎤ ⎥ ⎥ ⎦
Using two different finite difference schemes further delivers the equalities

η 1 i+1 j+1 -η 1 i j = a ∂η 1 ∂s i j + a 2 2 ∂ 2 η 1 ∂s 2 i j and η 1 i j -η 1 i-1 j-1 = a ∂η 1 ∂s i j - a 2 2 ∂ 2 η 1 ∂s 2 i j (4)
with s therein the curvilinear abscissa along the beam, and

θ 1 i+1 j+1 = θ 1 i + j + a ∂θ 1 ∂s i + j , θ 1 i - j = θ 1 i-1 j-1 + a ∂θ 1 ∂s i - j (5)
We use the Bernoulli model including a local rotation, in addition to the displacement acting as the translational d.o.f. The rotation is limited to the zero order (one order less than the displacement), so that previous expressions in (5) entail

θ 1 i+1 j+1 = θ 1 i + j and θ 1 i - j = θ 1 i-1 j-1 (6)
The continuity of the angular variation of the beam at the central node delivers the condition:

θ 1 i + j = θ 1 i - j = θ 1 i j (7)
Thus W 1 takes the form:

W 1 = ⎡ ⎣ 1 2 K c 1 ⎛ ⎝ a ∂η 1 ∂s i j + a 2 2 ∂ 2 η 1 ∂s 2 i j -a θ 1 i j 2 + a ∂η 1 ∂s i j - a 2 2 ∂ 2 η 1 ∂s 2 i j -a θ 1 i j 2 ⎞ ⎠ ⎤ ⎦ (8)
Due to Bernoulli kinematic assumption, it holds that θ 1

i j = ∂η 1 ∂s | i j . Thus, W 1 simplifies to W 1 = K c 1 a 2 2 ∂ 2 η 1 ∂s 2 i j 2
Taking both stiffnesses equal, that is,

K c 1 = K c 2 = K c further gives W 1 = K c a 2 2 ∂ 2 η 1 ∂s 2 i j 2 (9)
One next develops the second component of the energy

W 2 = ⎡ ⎢ ⎢ ⎣ 1 2 K f 2 θ 2 i+1 j-1 -θ 2 i j + 2 + θ 2 i j - -θ 2 i-1 j+1 2 + 1 2 K c 2 η 2 i+1 j-1 -η 2 i j -a 2 θ 2 i+1 j-1 + θ 2 i j + 2 + η 2 i j -η 2 i-1 j+1 -a 2 θ 2 i j - + θ 2 i-1 j+1 2 ⎤ ⎥ ⎥ ⎦ (10)
Using similarly two distinct finite difference schemes delivers

η 2 i+1 j-1 -η 2 i j = a ∂η 2 ∂s i j + a 2 2 ∂ 2 η 2 ∂s 2 i j ; η 2 i j -η 2 i-1 j+1 = a ∂η 2 ∂s i j - a 2 2 ∂ 2 η 2 ∂s 2 i j (11)
with according to Bernoulli model and the continuity of the angular variation of the beam at the central node the following conditions:

θ 2 i+1 j+1 = θ 2 i + j ; θ 2 i - j = θ 2 i-1 j-1 and θ 2 i + j = θ 2 i - j = θ 2 i j (12)
One then obtains after development the strain energy for the second beam

W 2 = K c a 2 2 ∂ 2 η 2 ∂s 2 i j 2 (13)
The total energy of the pantograph over the unit cell writes accordingly

W = a 4 4 K c ⎛ ⎝ ∂ 2 η 1 ∂s 2 i j 2 + ∂ 2 η 2 ∂s 2 i j 2 ⎞ ⎠ (14)
The transition is next made from the local basis (attached to the beam) to the Cartesian basis: the local components η 1 i j and η 2 i j project along the two components [

u 1 v 1 ].
After development of Eq. ( 14), one obtains (details are in the Appendix, expression

Eq. ( 13)):

W = a 4 2 K c cos β ∂ 2 v 1 ∂x 2 + ∂ 2 v 1 ∂y 2 2 + sin β ∂ 2 u 1 ∂x 2 + ∂ 2 u 1 ∂y 2 2 (15)
For the pantograph, the orientation of the first beam is β = π 4 , thus previous expression of the energy becomes This entails the expression of the energy density by dividing by the area of the elementary unit cell, that is S = (

W = a 4 2 K c √ 2 2 ∂ 2 v 1 ∂x 2 + ∂ 2 v 1 ∂y 2 2 + √ 2 2 ∂ 2 u 1 ∂x 2 + ∂ 2 u 1 ∂y 2 2 = a 4 4 K c ∂ 2 v 1 ∂x 2 + ∂ 2 v 1 ∂y 2 2 + ∂ 2 u 1 ∂x 2 + ∂ 2 u 1 ∂y 2 2 (16)
√ 2a) 2 , viz. W/S = a 2 8 K c ∂ 2 v 1 ∂x 2 + ∂ 2 v 1 ∂y 2 2 + ∂ 2 u 1 ∂x 2 + ∂ 2 u 1 ∂y 2 2 (17)
This expression only involves second order gradient contributions; one notes as expected from the geometry of the system the symmetry with respect to the displacements components u 1 and v 1 .

The constitutive law of the second order equivalent medium can easily be obtained by computing the second order derivatives of the postulated strain energy density, viz.

W = W ε i j , K i jk = 1 2 σ i j ε i j + S i jk K i jk , with K i jk = ε jk,i ≡ 1 2 ∂ 2 u j ∂x k ∂x i + ∂ 2 u j ∂x j ∂x i (18)
This entails based on the comparison of ( 17) with ( 18) that the Cauchy stress vanishes, thus a purely second order grade effective continuum is obtained, with the hyperstress components given by

w = 1 2 K f + a 2 4 K c ∂ 2 v 1 ∂x 2 + ∂ 2 v 1 ∂y 2 2 + ∂ 2 u 1 ∂x 2 + ∂ 2 u 1 ∂y 2 2 → S 11 : = ∂w ∂ ∂ 2 u 1 ∂x 2 = K f + a 2 2 K c ; S 22 := ∂w ∂ ∂ 2 v 1 ∂y 2 = K f + a 2 2 K c S 12 : = ∂w ∂ ∂ 2 u 1 ∂y 2 = K f + a 2 2 K c ; S 21 := ∂w ∂ ∂ 2 v 1 ∂x 2 = S 12
The next section is devoted to the presentation of the discrete homogenization technique, which is extended to derive an effective second order grade continuum for planar lattices.

Discrete homogenization up to second order

We shall account for the second order gradient of the displacement and derive the expression of the hyperforces and hyperstress tensor, based on the virtual power form of the equilibrium equations; this shall later on entail the expression of the energy density.

Expressions of forces

We consider an inextensible 2D beam with length l b = εL b , developing internal work under the action of the transverse forces T O , T E and the moments M O , M E ; the subscripts O and E refer to the origin and extremity nodes of the beam (Fig. 3).

From beam theory and using a so called complete Bernoulli model (Dos Reis, 2010; Dos Reis & Ganghoffer, 2012b) with inextensible beams, we will determine the moment expressions accounting for the second order gradient of the displacement; we summarize below the expressions of the transverse forces and moments in vector format: T

T εb O = K b c e b⊥ • (D ε O -D ε E ) + l ε b 2 (φ ε E + φ ε O ) (19)
εb E = -T εb O (20) M εb O = K b c l ε b 6 l ε b (2φ ε O + φ ε E ) + 3 • e b⊥ • (D ε O -D ε E ) (21) M εb E = K b c l ε b 6 l ε b (φ ε O + 2φ ε E ) + 3 • e b⊥ • (D ε O -D ε E ) (22) with D ε O , D ε E the displacement

Asymptotic expansion of the kinematic variables D ε

One considers displacement fields parameterized by the curvilinear abscissa (denoted by s in the sequel):

D = (U,V ) φ = φe 3 It holds in Cartesian coordinates D ε (x, y) = D O (x) + εD 1( x, y) + ε 2 D 2( x, y) + • • • (23)
The same displacement (vector field) writes in curvilinear coordinate along the beam

D ε = D O (s ε ) + εD 1 (s ε ) + ε 2 D 2 (s ε ) + • • • (24)
Specializing previous displacement field at the origin and extremity nodes gives the Taylor series expansions:

D ε O = D O (s ε ) + εD O 1 (s ε ) + ε 2 D O 2 (s ε ) + • • • (25) D ε E = D O (s ε + εL i δ ib ) + εD E 1 (s ε + εL i δ ib ) + ε 2 D E 2 (s ε + εL i δ ib ) + • • • (26) = D O (s ε ) + εL i δ ib ∂D O (s ε ) ∂s i + ε 2 L 2 i δ 2 ib 2 ∂ 2 D O (s ε ) ∂s 2 i + • • • + εD E 1 (s ε ) + ε 2 L i δ ib ∂D E 1 (s ε ) ∂s i + • • • + ε 2 D E 2 (s ε ) + • • • (27)
This entails the relative displacement

(D ε O -D ε E ) = ε D O 1 -D E 1 -L i δ ib ∂D o (s ε ) ∂s i + ε 2 D O 2 -D E 2 -L i δ ib ∂D E 1 (s ε ) ∂s i - L 2 i δ 2 ib 2 ∂ 2 D O (s ε ) ∂s 2 i ( 28
)
with L i the periodicity length (length of the repetitive unit cell, see Fig. 3), δ i the shift factor equal to ±1) for nodes belonging to a neighboring cell, and nil for nodes located inside the considered cell, and i ∈ {1, 2} indicating the considered axis e 1 or e 2 .

The cell is equilibrated at each node, and the efforts are periodical, with a period length equal to the cell width; this means each boundary node n has a mirror node on an opposite edge of the unit cell, as pictured in Fig. 4: one such node belongs to the elementary cell and is labeled n e , whereas the second node belongs to the adjacent cell, and is called n a . The first node type is an origin node for the beam, whereas the second node type is an extremity node.

One notices that the relative displacement expresses versus both first and second order gradients of the displacement field.

Similarly, one writes the asymptotic expansion of the microrotations, restricting presently to the first order:

φ ε = φ O (s ε ) Thus, φ ε O = φ O (s ε ), φ ε E = φ O (s ε ) = ∂D O (s ε ) ∂s i e b⊥
Inserting Eq. ( 28) into Eqs. ( 19)-( 22), one obtains the following expressions of the transverse forces and moments:

T O = K b c ε e b⊥ • D O 1 -D E 1 -L i δ ib ∂D O (s ε ) ∂s i + L b φ O +ε 2 e b⊥ . D O 2 -D E 2 -L i δ ib ∂D E 1 (s ε ) ∂s i - L 2 i δ 2 ib 2 ∂ 2 D O (s ε ) ∂s 2 i (29) T E = -T O (30) M O = K b c L b 6 ε 2 L b (3φO) + 3e b⊥ D O 1 -D E 1 -L i δ ib ∂D O (s ε ) ∂s i +ε 3 3e b⊥ D O 2 -D E 2 -L i δ ib ∂D E 1 (s ε ) ∂s i - L 2 i δ 2 ib 2 ∂ 2 D O (s ε ) ∂s 2 i (31) M E = K b c L b 6 ε 2 L b (3φO) + 3e b⊥ D O 1 -D E 1 -L i δ ib ∂D O (s ε ) ∂s i +ε 3 3e b⊥ D O 2 -D E 2 -L i δ ib ∂D E 1 (s ε ) ∂s i - L 2 i δ 2 ib 2 ∂ 2 D O (s ε ) ∂s 2 i ( 32
)
The following notations are conveniently introduced:

V O1 = D O 1 • e b⊥ (33) V E1 = D E 1 + L i δ ib ∂D O (s ε ) ∂s i • e b⊥ (34) V O2 = D O 2 . e b⊥ (35) V E2 = D E 2 + L i δ ib ∂D E 1 (s ε ) ∂s i + L 2 i δ 2 ib 2 ∂ 2 D O (s ε ) ∂s 2 i .e b⊥ (36)
This leads to the following simpler expressions of the forces and moments:

T O = K b c ε V O1 -V E1 + L b φ O + ε 2 (VO2 -V E2) (37) T E = -T O (38) M O = K b c L b 6 ε 2 L b (3φO) + 3(V O1 -V E1) + ε 3 (3(VO2 -V E2)) (39) M E = M O (40)
We next determine the derivatives of the displacements, expressions

L i δ ib ∂D O (s ε ) ∂s i , L i δ ib ∂D E 1 (s ε ) ∂s i and L 2 i δ 2 ib 2 ∂ 2 D O (s ε ) ∂s 2 i in a Cartesian
basis. After resolution and simplifications, one obtains:

L i δ ib ∂D O (s ε ) ∂s i = L 1 δ 1 cos θ 1 ∂D O ∂x + sin θ 1 ∂D O ∂y + L 2 δ 2 cos θ 2 ∂D O ∂x + sin θ 2 ∂D O ∂y (41) L i δ ib ∂D E 1 (s ε ) ∂s i = L 1 δ 1 cos θ 1 ∂D E 1 ∂x + sin θ 1 ∂D E 1 ∂y + L 2 δ 2 cos θ 2 ∂D E 1 ∂x + sin θ 2 ∂D E 1 ∂y (42) L 2 i δ 2 ib 2 ∂ 2 D O (s ε ) ∂s 2 i = ⎛ ⎝ L 2 1 δ 2 1 2 cos 2 θ 1 ∂ 2 D O ∂x 2 + sin 2 θ 1 ∂ 2 D O ∂y 2 + 2sinθ 1 cosθ 1 ∂ 2 D O ∂x ∂y + L 2 2 δ 2 2 2 cos 2 θ 2 ∂ 2 D O ∂x 2 + sin 2 θ 2 ∂ 2 D O ∂y 2 + 2sinθ 2 cosθ 2 ∂ 2 D O ∂x ∂y ⎞ ⎠ (43)
The trigonometric factors cos θ 1 , sin θ 1 , cos θ 2 and sin θ 2 are respectively the components of the periodicity vectors

Y 1 = [ cos θ 1 sin θ 1 ] (i, j) and Y 2 = [ cos θ 2 sin θ 2 ] (i, j) (

Homogenization step

Writing of the equilibrium equations in virtual power form

The equilibrium equations write in virtual power form as (self-equilibrium in the absence of external forces)

b∈B(b) T b . V = 0 (44)
with .

V a virtual velocity field chosen to vanish on the edges of the domain. The moment equilibrium applied to the lattice nodes writes (Dos Reis & Ganghoffer, 2012b)

v i ∈Z 2 b∈B R M b O • w b O + M b E • w b E = 0 (45)
with w the virtual rate of rotation. This equation will prove useful later on to solve the localization problem for the kinematic unknowns.

The virtual power of internal forces for the whole lattice expresses accordingly as the sum of a product of internal forces by their kinematic dual quantities at each node:

P = b∈B(b) T i . V i (46)
The virtual velocity field is of the same form as the displacement field evaluated in previous section. In the sequel, we select a field having a simple form, reducing to the zero order term in ε.

The virtual velocity field is then written as:

. D ε (P) = . D O (P) + • • • (47)
Vector P therein represents the coordinate vector of the considered point of the global lattice. Previous velocity field in (47) can be decomposed into longitudinal and transverse components; the virtual transverse velocity writes:

V ε (P) = Ḋε (P) • e b⊥ = . V O (P) + • • • (48)
One may decompose the sum in ( 46) as a double summation on the cells and the nodes of a given cell, viz.

P = c∈Z n∈B T i . V i = c∈Z P e ( 49
)
with i = (c, n), Z being the set of all lattice cells, B the set of all nodes a given elementary cell, and P e the virtual power of an elementary cell. One can thus write

P e = n∈B T i . V i (50)
Since the sum of efforts cancel at each lattice node, it only remains the edge efforts for each node since the virtual power of internal forces mutually cancel:

P e = b T O . V O +T E . V E (51)
Since each beam is equilibrated, it holds the equality T b E = -T b O ; this entails the following simplification of the previous virtual power in (51) as:

P e = b T E . V E - . V O (52)
We next expand the relative virtual velocity ( .

V E - . V O ). A Taylor series expansion of (48) gives . V O (s ε ) = . V O (s ε ) + • • • . V E = . V O (s ε + εL i δ ib ) = . V O +εL i δ ib ∂ . V O (s ε ) ∂s i + ε 2 L 2 i δ 2 ib 2 ∂ 2 . V O (s ε ) ∂s 2 i + • • • Thus ( . V E - . V O ) = εL i δ ib ∂ . V O (s ε ) ∂s i + ε 2 L 2 i δ 2 ib 2 ∂ 2 . V O (s ε ) ∂s 2 i ( 53
)
The expression of the virtual external power of an elementary cell in (52) then becomes

P e = b T E εL i δ ib ∂ . V O (s ε ) ∂s i + ε 2 L 2 i δ 2 ib 2 ∂ 2 . V O (s ε ) ∂s 2 i = b -K b c ε V O1 -V E1 + L b φ O + ε 2 (VO2 -V E2) εL i δ ib ∂ . V O (s ε ) ∂s i + ε 2 L 2 i δ 2 ib 2 ∂ 2 . V O (s ε ) ∂s 2 i (54)
In order to further simplify expression (54), one introduces the following notations:

T E = -k b c ε V O1 -V E1 + L b φ O + ε 2 (VO2 -V E2) = εT 1 E + ε 2 T 2 E (55)
thus leading to

P e = b ε 2 T 1 E L i δ ib ∂ . V O (s ε ) ∂s i + ε 3 T 1 E L 2 i δ 2 ib 2 ∂ 2 . V O (s ε ) ∂s 2 i + T 2 E L i δ ib ∂ . V O (s ε ) ∂s i + ε 4 T 2 E L 2 i δ 2 ib 2 ∂ 2 . V O (s ε ) ∂s 2 i (56)
One recalls the expression of the discrete form of the virtual power of forces for the entire lattice: Considering that the elementary cell is a small surface element

dv = gds → ds = 1 g dv,
with s therein the curvilinear coordinate and g the determinant of the Jacobian of the coordinate transformation: , j) , as pictured in Fig. A.1. This entails the continuous form of the virtual power of internal forces:

g = det ∂R[1] ∂s 1 ∂R[1] ∂s 2 ∂R[2] ∂s 1 ∂R[2] ∂s 2

R is the vector function for the coordinate change between the Cartesian and the curvilinear basis, expressed in the form

R = [L 1 s 1 Y 1 [1] + L 2 s 2 Y 2 [1], L 1 s 1 Y 1 [2] + L 2 s 2 Y 2 [2]] (i
lim ε→0 P = 1 g b ε 2 T 1 E L i δ ib ∂ . V O (s ε ) ∂s i + ε 3 T 1 E L 2 i δ 2 ib 2 ∂ 2 . V O (s ε ) ∂s 2 i + T 2 E L i δ ib ∂ . V O (s ε ) ∂s i + ε 4 T 2 E L 2 i δ 2 ib 2 ∂ 2 . V O (s ε ) ∂s 2 i dv (59)
One denotes the apparition of terms in three orders in ε: the zero order term shall lead to the Cauchy stress, the first order to the coupling tensors between first and second order effects, and the second order term to the hyperstress tensor. One restricts in the sequel to periodic structures with central symmetry, thus the coupling tensors vanish.

From Eqs. ( 41) and ( 43), one can develop the expressions

∂ . V O (s ε ) ∂s i , L i δ ib ∂ 2 . V O (s ε ) ∂s 2 i
and insert them back into (59):

lim ε→0 P = ⎡ ⎢ ⎢ ⎢ ⎣ 1 g b ⎛ ⎜ ⎜ ⎜ ⎝ T 1 E L 1 δ 1 cos θ 1 ∂ . V O ∂x + sin θ 1 ∂ . V O ∂y + L 2 δ 2 cos θ 2 ∂ . V O ∂x + sin θ 2 ∂ . V O ∂y +ε 2 ⎛ ⎝ T 2 E ⎛ ⎝ L 2 1 δ 2 1 2 cos 2 θ 1 ∂ 2 . V O ∂x 2 + sin 2 θ 1 ∂ 2 . V O ∂y 2 + 2sinθ 1 cosθ 1 ∂ 2 . V O ∂x ∂y + L 2 2 δ 2 2 2 cos 2 θ 2 ∂ 2 . V O ∂x 2 + sin 2 θ 2 ∂ 2 . V O ∂y 2 + 2sinθ 2 cosθ 2 ∂ 2 . V O ∂x ∂y ⎞ ⎠ ⎞ ⎠ ⎞ ⎟ ⎟ ⎟ ⎠ ⎤ ⎥ ⎥ ⎥ ⎦ dv (60)

Equivalence with a second gradient continuum

The principle of virtual power is involved to establish the equilibrium equations used for the determination of the stress and hyperstress tensors. Let recall the constitutive law for an anisotropic second gradient continuum: σ i j = C i jpq pq + M pqri j K pqr S i jk = M i jkpq pq + A i jkpqr K pqr introducing therein σ i j , S i jk , pq et K pqr the stress, hyperstress, deformation and gradient of deformation tensors respectively.

The tensors C i jpq , A i jkpqr , M pqri j are the first and second order elasticity tensors and the coupling tensor respectively. For periodic uniform structures endowed with central symmetry, the coupling tensor M vanishes.

Let then determine the expressions of the stress and hyperstress tensors, starting from the internal power of stresses for a second order grade continuum in the case of linear elasticity, to highlight the comparison with Eq. ( 60):

P i = σ -S • ∇ • ∇ • . D dv (61)
An expansion of previous expression in index form leads to: Dr e r

P i = σ i j e i ⊗ e j -S
∂x p ∂x q dv = σ • e q • ∂ . D ∂x q - S • e p • e q • ∂ 2 . D ∂x p ∂x q dv (62)
Previous expression can be recast under the form

P i = F q • ∂ . D ∂x q -H pq • ∂ 2 . D ∂x p ∂x q dv (63)
introducing therein the force vector we use here a slightly different convenient notation by underlining tensors which are indexed by additional superscript

F q = σ • e q = σ i j e i ⊗ e j • e q = σ iq e i ( 64
)
and the force hyper-vector

H pq = S • e p • e q = S klm e k ⊗ e l ⊗ e m • e p • e q = S kqp e k (65)
Based on Eqs. ( 64) and ( 65), one can next reconstruct the stress and hyperstress tensors σ and S σ = σ iq e i ⊗ e q = F q ⊗ e q (66) S = S kqp e k ⊗ e q ⊗ e p = H pq ⊗ e q ⊗ e p (67)

Let introduce vectors

T 1 E = T 1 E e b⊥ et T 2 E = T 2 E e b⊥
, in order to formulate Eq. ( 60) under the form

lim ε→0 P = b F q • ∂ . DO ∂x q -H pq • ∂ 2 . DO ∂x p ∂x q dv (68) with q ∈ {1, 2}, F 1 = T 1 E g (L1δ1cosθ1+L2δ2cosθ2) (69) F 2 = T 1 E g (L1δ1sinθ1+L2δ2sinθ2) ( 70
)
with g = L 1 L 2 (cosθ 1 sinθ 2sinθ 1 cosθ 2 ) the determinant of the Jacobian matrix and the pair (p, q) ∈ {(1, 1), (2, 2), (1, 2)}.

There are three combinations (instead of four), since one makes the summation on the combination (1,2) and (2,1) of the crossed derivative: Table 1 Geometrical and mechanical data for the pantographic structure shown in The automatic treatment of the discrete homogenization method is implemented in a MAPLE code, allowing the direct computation of the stress and hyperstress tensors; this will enable the evaluation of the strain energy density. The method allows to treat any centrosymmetric lattice in a systematic manner. One obtains a vanishing Cauchy stress at the first order, viz.

H 11 = H 11 = T 2 E g L 2 1 δ 2 1 cos 2 θ 1 2 + L 2 2 δ 2 2 cos 2 θ 2 2 (71) H 22 = T 2 E g L 2 1 δ 2 1 sin 2 θ 1 2 + L 2 2 δ 2 2 sin 2 θ 2 2 (72) H 12 = T 2 E g L 2 1 δ 2 1 cosθ 1 sinθ 1 +L 2 2 δ 2 2 sinθ 2 cosθ 2 (73)

σ = 0

The hyperstress third order tensor S is obtained as (without accounting for periodicity, thus here the shift factor δ must be removed so we will put it equal to 1 and the periodicity length L i is replaced by the length of each beam L b in Eqs. ( 26)-( 28), then the same steps are followed):

S[1..2, 1..2, 1] = a 2 4 k c ∂ 2 u 1 ∂x 2 + ∂ 2 u 1 ∂y 2 ∂ 2 v 1 ∂x 2 + ∂ 2 v 1 ∂y 2 0 0 , S[1..2, 1..2, 2] = 2 4 k c 0 0 ∂ 2 u 1 ∂x 2 + ∂ 2 u 1 ∂y 2 ∂ 2 v 1 ∂x 2 + ∂ 2 v 1 ∂y 2
This entails the following expression of the strain energy density:

W/S = a 2 8 K c ∂ 2 v 1 ∂x 2 + ∂ 2 v 1 ∂y 2 2 + ∂ 2 u 1 ∂x 2 + ∂ 2 u 1 ∂y 2 2
The same expression of the strain energy density characteristic of a pure second order medium is obtained by both methods.

Determination of the internal lengths

In order to identify the internal lengths associated to the different deformation modes (extension and in-plane shear in the present case), one can rely on the equilibrium equations (in the absence of body forces)

- ∂σ k1 ∂x 1 + ∂σ k2 ∂x 2 + ∂ 2 S k11 ∂x 1 ∂x 1 + ∂ 2 S k12 ∂x 1 ∂x 2 + ∂ 2 S k21 ∂x 2 ∂x 1 + ∂ 2 S k22 ∂x 2 ∂x 2 = 0, k = 1, 2,
inserting therein the obtained expressions of the stress and hyperstress tensors. Alternatively, one can use the obtained strain energy density

W = W ε i j , K i jk = 1 2 σ i j ε i j + S i jk K i jk , with K i jk = ε jk,i
l 2 1 = 2a λ + 2μ
In a similar manner, one elaborates the contribution 1 2

(σ12ε12 + S 122 K 122) = 1 2 {2μ + a 345 ∇ 2 }ε 12
leading to the identification of the internal length in shear l 2 2 = a 345 2μ [START_REF] Papanicolopulos | Continua with microstructure: Part II: second-gradient theory[END_REF]). In the present case of the pantograph, isolating contributions from the strain energy density leads to the successive identification of the following internal lengths in extension and shear: 1 2

(σ11ε11 + S 111 K 111 ) = 1 2 K c 0 + a 2 4 ∇ x ∂u 1 ∂x + ∂ 2 u 1 ∂y 2 ⇒ l 2 xx = a 2 4 1 2 (σ22ε22 + S 222 K 222) = 1 2 K c 0 + a 2 4 ∇ y ∂v 1 ∂y + ∂ 2 v 1 ∂x 2 ⇒ l 2 yy = a 2 4 1 2 (σ12ε12 + S 122 K 122) = 1 2 K c 0 + a 2 4 ∇ y ∂u 1 ∂y + ∂ 2 u 1 ∂x 2 ⇒ l 2 xy = l 2 yx = a 2 4
For more complex structures, the asymptotic homogenization is quite efficient (analytical methods become quite complicated and inapplicable).

Simulations at the structural level based on the equivalent second gradient continuum

In order to regularize the second order gradient continuum, we incorporate small but finite stretches of the two beams of the pantograph structure, leading to an additional first order energy density.

Expression (1) of the strain energy of the structure accordingly becomes: Following the same steps as in Section 2, we obtain the following continuum strain energy density:

W = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 1 2 K f 1 θ 1 i+1 j+1 -θ 1 i j + 2 + θ 1 i j - -θ 1 i-1 2 + 1 2 K f 2 θ 2 i+1 j-1 -θ 2 i j + 2 + θ 2 i j - -θ 2 i-1 j+1 2 + 1 2 K c 1 η 1 i+1 j+1 -η 1 i j -a 2 θ 1 i+1 j+1 + θ 1 i j + 2 + η 1 i j -η 1 i-1 j-1 -a 2 θ 1 i j - + θ 1 i-1 j-1 2 + 1 2 K c 2 η 2 i+1 j-1 -η 2 i j -a 2 θ 2 i+1 j-1 + θ 2 i j + 2 + η 2 i j -η 2 i-1 j+1 -a 2 θ 2 i j - + θ 2 i-1 j+1 2 + 1 2 EI L θ 1 i j + -θ 2 i j + 2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦
W/S = a 2 8 K c ∂ 2 v 1 ∂x 2 + ∂ 2 v 1 ∂y 2 2 + ∂ 2 u 1 ∂x 2 + ∂ 2 u 1 ∂y 2 2 + K f 2a 2 ∂u 1 ∂x + ∂u 1 ∂y 2
In Fig. 5, a picture of a real pantographic structure is shown; it gives a motivation of the micro-macro identification procedure presented in this paper; the bending stiffness of the beams is uniform and has been selected as K b c = 7.85 × 10 -3 N m. It has to be remarked that numerical simulations using micro models based on Cauchy 3D continua require few days of computations on a powerful workstation. By comparison, the simulations performed using the homogenized model constructed in the present paper were obtained with the same workstation and the same software (Comsol Multiphysics) in a few minutes. Note that the 2D continuum obtained with the homogenization procedure assumes the existence of two families of inextensible fibers with no relative slip (for more details about these continuum models, see dell'Isola, Della Corte, Greco, & Luongo, 2015c): therefore, the minimization process which we used to get all shown equilibrium configurations requires the method of Lagrange multipliers. The boundary conditions imposed in all numerical simulations consist in imposed displacement (hard device boundary conditions). The normal derivatives of displacement are left arbitrary and therefore the dual double forces are assumed to vanish.

The geometrical and mechanical data for the pantographic shown in Fig. 5 are summarized in Table 1. In Fig. 6, the 2D equivalent model constructed in previous section is used to obtain the deformed shape in extension of a rectangular specimen having a side three times longer than the other. The same figure shows the actual shape of inextensible fibers and a scale of colors measures the density of the second gradient deformation energy, which has been computed as related to the microscopic bending of inextensible fibers.

It appears from Fig. 7 that the deformation energy is concentrated in some transitions zones and at wedges of the specimen: this reflects the occurrence, in 2D second gradient materials, of contact forces concentrated in points. Fig. 8 shows the nearly piecewise constant shear deformation field in the case of a sheared specimen; we observe the occurrence of five shear regions (denoted by green, blue and red), with relative boundary layers: this circumstance could not be accounted for without a second gradient model.

The surface double traction and the edge tractions may be evaluated from the principle of virtual work [START_REF] Germain | The method of virtual power in continuum mechanics. Part 2: Microstructure[END_REF][START_REF] Papanicolopulos | Continua with microstructure: Part II: second-gradient theory[END_REF], using integration by part in the virtual work of internal forces (suggesting a similar expression of the virtual work of external forces), viz. the equalities

-δW int + δW ext = 0 δW int = V σ i j δε i j + S i jk δK i jk dV δW ext = V f k δu k dV + S P k δu k dS + S R k Dδu k dS + C E k δu k dS
The quantities P k , R k , E k therein are respectively the surface tractions, the surface double tractions and the edge tractions. The body forces f k are balanced by the first order and second order stresses in the equilibrium equation σ jk -S i jk,i , j + f k = 0 Classical derivations lead to the identification of the external forces versus Cauchy stress and hyperstress (the double bracket [[.]] stands for the jump of the enclosed quantity) as P k = n j σ jk -S i jk ,i -D j n i S i jk + (D l n l )n j n i S i jk R k = n i n j S i jk E k = s m e ml j n l n j S i jk with e ml j the components of the permutation tensor, and s m the components of the surface vector. The occurrence of these forces involving contributions which are specific to second gradient continua leads to specific features like boundary layers and concentrated forces, as will be exemplified in the sequel thanks to numerical simulations. Fig. 9 shows the deformation energy field relative to a combined bending and compression test. One can observe the onset of three boundary layers, the major of which involves a series of sub-layers characterized by an increasing level of deformation energy (pictured in red, yellow, green and cyan in Fig. 9), while the third one with lowest intensity cuts the first boundary layer being concentrated in a inextensible fiber orthogonal to the most deformed ones: this situation is better illustrated by Fig. 10.

The previous integrals building up the work of external forces have to be interpreted in the present 2D context (V denotes the pantograph domain, a surface, S a line and C a domain edge); the surface tractions and double tractions are distributed within the observed boundary layers (Figs. 9 and 10), while the concentrated force observed at the top right corner of the tested specimen (Figs. 9 and 10) corresponds to quantity E k expressed in the last expression.

Fig. 11 shows 3D plots of the deformation energy in the same in-plane bending and compression test: it is intended to depict the onset of concentrated contact forces. Fig. 12 shows the nearly piecewise constant pattern of the shear deformation field: one can identify four shear regions (orange, blue and red), separated by the corresponding transition zones with, once more, the regularizing effect of second gradient energy.

Finally, in Figs. 13 and 14, the effect of the presence of two symmetric holes is shown. In Fig. 13, a sort of deformation energy fringes is obtained together with the "dispersive" effect determined by the presence of holes, in the case of pure extension. In the case of compression and bending, the presence of holes exacerbates the intensity of the deformation energy in transition zones and the intensity of concentrated forces (Fig. 14).

Conclusion

In the present work, we have shown that the linearized homogenized model for the pantographic lattice described at microlevel by Fig. 1 must necessarily be a second gradient continuum. Some numerical simulations are finally shown, to illustrate some peculiarities of the obtained continuum models. The developed homogenization technique is general and is capable of treating any lattice, periodical or not. It provides the full set of first and second gradient effective moduli accounting for the full anisotropy of the effective continuum. For more complex structures, the asymptotic homogenization is quite efficient (in comparison, analytical methods become quite complicated and inapplicable).

If we consider the total expression of the rotation (5) including the gradient of the rotation, the energy density w := W/S takes now the more general form

w = 1 2 K f + a 2 4 K c ∂ 2 v 1 ∂x 2 + ∂ 2 v 1 ∂y 2 2 + ∂ 2 u 1 ∂x 2 + ∂ 2 u 1 ∂y 2 2
This constitutes an extension of previous Bernoulli beam model towards a second gradient continuum in which the fibers curvatures related to second order displacement gradients are present: as underlying beams forming the pantographic lattice are Euler beams, these gradients determine the rotation of beams fibers. In this sense, the considered model supplies an example of a constrained Cosserat continuum.

The well-posedness results available in the literature do not apply immediately to the energy we have determined via our homogenization procedure: it seems that some investigations about its mathematical nature may be needed.

However, the numerical simulations presented in this paper show that there is a reasonable large set of boundary conditions for which a suitable formulated integration scheme actually do converge. Few additional comments are here in order about the results of our simulations: the model is capable to describe (i) the onset of zones where deformation and deformation gradients are concentrated; (ii) the establishment of different segregated deformation zones (or phases); and the global effect of the presence of holes on the overall deformation pattern. All these features indicate that the modeling effort presented in this paper deserve further efforts.

One has to remark that the model which we obtained here via a micro-macro identification procedure has been heuristically introduced by Rivlin and Pipkin, already in the decades 1950-1960: for a review and a reformulation of the results of last authors, the reader shall refer to dell 'Isola et al. (2015c). Although the assumption of inextensible fibers has not a general validity, the results obtained in Section 4 are an important step towards the understanding of the mechanical behavior of the considered systems: we are confident that the generalization to the case of extensible, but stiff fibers will be a relatively easy task.

There is no unique way to generate second order grade continuum models from such discrete homogenization of lattices, since the effective continuum is clearly dependent on the initial static model at the beam level. Another alternative to the complete Bernoulli beam model which is more appropriate for bending dominated lattices (like the pantograph) would be to consider Bernoulli kinematics without rotation (appropriate for tension dominated lattices). Such an alternative is clearly of high interest and will be investigated in future contributions.

The consideration of the large displacements response in such pantographic structures constitutes a natural extension of the present work. It will not be possible, in general, to determine a closed form expression of the homogenized energy for the equivalent continuum; nevertheless, the determination via numerical computations of the effective stiffness coefficients will be helpful, and shall constitute an objective for future investigations. Moreover, in the more general context of large displacements and deformations of thin and light structures, may be of great interest further investigations about external non-conservative interactions as studied in [START_REF] Carassale | Effects of free-stream turbulence and corner shape on the galloping instability of square cylinders[END_REF]Marrè-Brunenghi (2013), Nguyen, Freda, Solari, and[START_REF] Nguyen | Aeroelastic instability and wind-excited response of complex lighting poles and antenna masts[END_REF], [START_REF] Pagnini | Reliability analysis of wind-excited structures[END_REF], [START_REF] Solari | A numerical algorithm for the aerodynamic identification of structures[END_REF]. These phenomena, involving stability issues, can be useful addressed using perturbation methods as e.g. described in [START_REF] Luongo | Mode localization by structural imperfections in one-dimensional continuous systems[END_REF][START_REF] Luongo | Mode localization in dynamics and buckling of linear imperfect continuous structures[END_REF]. In the present situation, one has Y 1 = [ cos θ 1 = 1 sin θ 1 = 0 ] (i, j) and Y 2 = [ cos θ 2 = 0 sin θ 2 = 1 ] (i, j) .

Inserting the present values of the trigonometric factors cos θ 1 , sin θ 1 , cos θ 2 and sin θ 2 , one obtains the transformation of the derivatives:

∂η ∂s (A.12) Thus, the energy writes: .13)

W = a 4 4 K c ⎛ ⎝ ∂ 2 η 1 ∂s 2 i j 2 + ∂ 2 η 2 ∂s 2 i j 2 ⎞ ⎠ = a 4 4 K c cos β ∂ 2 v 1 ∂x 2 + ∂ 2 v 1 ∂y 2 -sin β ∂ 2 u 1 ∂x 2 + ∂ 2 u 1 ∂y 2 2 + cos β ∂ 2 v 1 ∂x 2 + ∂ 2 v 1 ∂y 2 + sin β ∂ 2 u 1 ∂x 2 + ∂ 2 u 1 ∂y 2 2 = a 4 2 K c cos β ∂ 2 v 1 ∂x 2 + ∂ 2 v 1 ∂y 2 2 + sin β ∂ 2 u 1 ∂x 2 + ∂ 2 u 1 ∂y 2 2 (A

Fig. 1 .

 1 Fig. 1. Unit cell of the pantograph structure. The link at node (i,j) is a pivot.

Fig. 2 .

 2 Fig. 2. Kinematics of a beam element.

Fig. 3 .

 3 Fig. 3. Kinematics and statics for a single isolated beam e 3 .

Fig. 4 .

 4 Fig. 4. Elementary cell in a 2D space.

  vectors of the origin and extremity nodes of each beam respectively, l ε b the beam length, and K

  as pictured in Fig. A.1 in the Appendix).

  of this power is obtained, transforming the double summation into the Riemann integral in the domain when ε

.-

  klm e k ⊗ e l ⊗ e m • ∂e p Dr e r -S klm e k ⊗ e l ⊗ e m • S klm e k ⊗ e l ⊗ e m • e p • e q • ∂ 2 .

Fig. 5 .

 5 Fig. 5. Pantographic structure, by courtesy of Tomasz Lekszycki, Marek Pawlikowski and Roman Grygoruk. The rectangle has sides 7 cm × 21 cm. The bending stiffness of the beams is selected to be K b c = 7.85 × 10 -3 N m.

Fig

  in compression and bending Minimum: 7 mm has a linear distribution Maximum: 14 mm Bending stiffness of the beam K b c = 7.85 × 10 -3 N m

Fig. 6 .

 6 Fig. 6. Deformed shape in extension: the scale of colors measures bending energy. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 8 .

 8 Fig. 8. Shear deformation field in the case of sheared specimen.

Fig. 9 .Fig. 10 .

 910 Fig. 9. Bending and compression test: formation of internal boundary layers as determined by the plot of the deformation energy. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)

Fig. 11 .

 11 Fig. 11. The same deformation as in Fig. 5: 3D plots of deformation energy showing deformation concentrations in layers and corners.

Fig. 12 .

 12 Fig. 12. Shear deformation field showing regions of constant shear and transition zones. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)

Fig. 13 .

 13 Fig. 13. Effect of presence of holes on the strain energy field in extension.

Fig. 14 .

 14 Fig. 14. Effect of the presence of holes on deformation energy field in compression (top) and bending (bottom).

Appendix:(

  Fig. A.1. Change of basis.

  v 1sin β u 1) + ∂ ∂y (cos β v 1sin β u 1)

  dx 2 (cos β v 1sin β u 1) + d 2 dy 2 (cos β v 1sin β u 1) (cos β v 1 + sin β u 1) + d 2 dy 2 (cos β v 1 + sin β u 1)

	∂η 2 ∂s	i j	=	∂η 2 ∂x	i j	+	∂η 2 ∂y	i j	=	∂ ∂x	(cos β v 1 + sin β u 1) +	∂ ∂y	(cos β v 1 + sin β u 1)
			=	cos β	∂v 1 ∂x	+ sin β	∂u 1 ∂x	+ cos β	∂v 1 ∂y	+ sin β	∂u 1 ∂y	(A.10)
	∂ 2 η 1 ∂s 2		i j	=	∂ 2 η 1 ∂x 2	i j	+	∂ 2 η 1 ∂y 2	j i	=	d 2
				=	cos β	∂ 2 v 1 ∂x 2 -sin β	∂ 2 u 1 ∂x 2 + cos β	∂ 2 v 1 ∂y 2 -sin β	∂ 2 u 1 ∂y 2	(A.11)
	∂ 2 η 2 ∂s 2		i j	= dx 2 = ∂ 2 η 2 ∂x 2 i j + ∂ 2 η 2 ∂y 2 i j d 2 = cos β ∂ 2 v 1 ∂x 2 + sin β ∂ 2 u 1 ∂x 2 + cos β	∂ 2 v 1 ∂y 2 + sin β	∂ 2 u 1 ∂y 2
		=		cos β	∂v 1 ∂x	-sin β	∂u 1 ∂x	+ cos β	∂v 1 ∂y	-sin β	∂u 1 ∂y	(A.9)