Claude Boutin

Francesco Dell'isola

Ivan Giorgio

Luca Placidi

Linear pantographic sheets: Asymptotic micro-macro models identification

In this paper we consider linear pantographic sheets which in their natural configuration are constituted by two orthogonal arrays of straight fibers interconnected by internal pivots. We introduce a continuous model by means of a micro-macro identification procedure based on the asymptotic homogenization method of discrete media. The rescaling of the mechanical properties and of the deformation measures is calibrated in order to complies with the specific kinematics imposed by the quasi-inextensibility of the fibers together with the large pantographic deformability. The obtained high order continuum model shows interesting and exotic features, related to its extreme anisotropy but also to the sub-coercivity of its deformation energy. Some first numerical simulations are presented, showing that the model can account for experimental uncommon phenomena occuring in pantographic sheets. The paper focuses on the precise analysis and the understanding of the effective behaviour based on a well-calibration of the extension and bending phenomena arising at the local scale. In an upcoming work the analysis will be extended to oblique arrays, to some analytical solutions to proposed equations and to some further applications.

Introduction

In the study of generalised (multiscale or architectured) continua, two different and complementary approaches can be identified.

1. A phenomenological approach can applied at the macroscopic scale: it gives a general framework and in general it is based on variational principles and methods ( [START_REF] Mindlin | Second gradient of strain and surface-tension in linear elasticity[END_REF][START_REF] Toupin | Theories of elasticity with couple-stress[END_REF][START_REF] Germain | The method of virtual power in continuum mechanics. Part 2: Microstructure[END_REF][START_REF] Gouin | Variational principle involving the stress tensor in elastodynamics[END_REF][START_REF]Principes variationnels en fluide compressible et en magnétodynamique des fluides[END_REF][START_REF] Dell'isola | Variational principles are a powerful tool also for formulating field theories, volume 535 of Chapter in[END_REF][START_REF] Dell'isola | Boundary conditions at fluid-permeable interfaces in porous media: A variational approach[END_REF][START_REF] Lekszycki | Application of variational methods in analysis and synthesis of viscoelastic continuous systems[END_REF]); however, using this approach it is not possible to make the influence of microstructure on macrophenomena to become explicit. In other words, a model is established but the microstructured material to which it could be applied is unknown.

2. Another approach for which the scale change is made possible through an homogenization. This method has been developed in a vast literature (see e.g. [START_REF] Boutin | Rayleigh scattering in elastic composite materials[END_REF][START_REF] Pideri | A second gradient material resulting from the homogenization of an heterogeneous linear elastic medium[END_REF][START_REF] Bensoussan | Asymptotic analysis for periodic structures[END_REF][START_REF] Boutin | Generalized inner bending continua for linear fiber reinforced materials[END_REF][START_REF] Auriault | Homogenization of coupled phenomena in heterogenous media[END_REF][START_REF] Allaire | Homogenization and two-scale convergence[END_REF][START_REF] Alibert | Second-gradient continua as homogenized limit of pantographic microstructured plates: a rigorous proof[END_REF][START_REF] Dell'isola | Pantographic 2D sheets: Discussion of some numerical investigations and potential applications[END_REF][START_REF] Dell'isola | Plane bias extension test for a continuum with two inextensible families of fibers: a variational treatment with Lagrange Multipliers and a perturbation solution[END_REF]) but it has not been applied too often to the case of multiscale materials in which a great contrast in physical properties appears at micro-level (see [START_REF] Soubestre | Non-local dynamic behavior of linear fiber reinforced materials[END_REF][START_REF] Boutin | Generalized beams and continua. dynamics of reticulated structures[END_REF]). Moreover homogenization methods treat more particular cases and therefore they are not supposed to cover all situations. However a careful homogenization technique gives the possibility to identify the micro-morphologies that lead to the macroscopic behavior characteristic of generalized continua.

Recently the attention has been attracted by a particular class of microstructure: that which produces so called pantographic continua ( [START_REF] Dell'isola | Designing a light fabric metamaterial being highly macroscopically tough under directional extension: first experimental evidence[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][START_REF] Dell'isola | Pantographic 2D sheets: Discussion of some numerical investigations and potential applications[END_REF][START_REF] Dell'isola | Plane bias extension test for a continuum with two inextensible families of fibers: a variational treatment with Lagrange Multipliers and a perturbation solution[END_REF][START_REF] Giorgio | Buckling modes in pantographic lattices[END_REF][START_REF] Madeo | Wave propagation in pantographic 2D lattices with internal discontinuities [lainelevi kahemõõtmelises sisemiste katkevustega pantograafilises võres[END_REF]). This kind of structure is inspired by several natural examples, indeed, some biological tissues present fibers that can be modeled in a similar way (see e.g. [START_REF] Melnik | Dynamic fiber reorientation in a fiber-reinforced hyperelastic material[END_REF][START_REF] Federico | Elasticity and permeability of porous fibre-reinforced materials under large deformations[END_REF][START_REF] Grillo | Remodelling in statistically oriented fibrereinforced materials and biological tissues[END_REF] for some possible applications). Their interest was initially related to the possibility to prove the existence of purely second gradient continua [START_REF] Alibert | Truss modular beams with deformation energy depending on higher displacement gradients[END_REF], but subsequently their practical applicability has been proven for woven fabrics and some 'ad hoc' designed metamaterials (see e.g. [START_REF] Vescovo | Dynamic problems for metamaterials: review of existing models and ideas for further research[END_REF][START_REF] Bîrsan | Deformation analysis of functionally graded beams by the direct approach[END_REF][START_REF] Eremeyev | Existence theorems in the linear theory of micropolar shells[END_REF][START_REF] Altenbach | Micropolar shells as two-dimensional generalized continua models[END_REF][START_REF] Dell'isola | Designing a light fabric metamaterial being highly macroscopically tough under directional extension: first experimental evidence[END_REF][START_REF] Eremeyev | On effective properties of materials at the nano-and microscales considering surface effects[END_REF]).

On the other hand some theoretical interesting results were presented by [START_REF] Chesnais | Wave propagation and non-local effects in periodic frame materials: Generalized continuum mechanics[END_REF][START_REF] Boutin | Generalized beams and continua. dynamics of reticulated structures[END_REF][START_REF] Boutin | Generalized inner bending continua for linear fiber reinforced materials[END_REF]. In these last papers some suitable homogenization methods were introduced to study the dynamics of periodic beam structures and other evidence has been presented concerning the need to introduce second gradient continua when high contrast of mechanical properties is present at micro-level. In all considered cases, directional (anisotropic) materials with high contrast in properties between shear and extension are studied. This paper, differently to what done in [START_REF] Rahali | Homogenization à la Piola produces second gradient continuum models for linear pantographic lattices[END_REF], considers the case of extensible fibers and gives a more solid foundation to and generalizes the heuristic results presented in [START_REF] Placidi | Identification of two-dimensional pantographic structure via a linear D4 orthotropic second gradient elastic model[END_REF], based on an accurate analysis of the different and relative order of magnitude of involved physical phenomena.

In this paper we will use a micro-macro asymptotic identification method and obtain the macroscopic equilibrium equations for pantographic lattices in the neighborhood of a reference configuration. The analysis of obtained equation is started and some equilibrium problems are solved by means of numerical simulations.

In a forthcoming paper we will treat the generalized case of reference configuration constituted by two oblique arrays of straight beams, some semi-analytical solutions and some further application.

Hypotheses and notations

Let us consider a periodic pantograph network of fibers (which we also call pantographic lattice or pantographic sheet) formed of two families of continuous fibers arranged perpendicularly and along the axes x and y. The fibers oriented along x are identical, as well as those oriented along y. However, the two families may differ from one another. These fibers are connected by perfect pivots, with an axis perpendicular to the {x, y} plane. The fibers oriented along x are spaced periodically by the length y . Those oriented along y, are spaced by x . This defines the rectangular mesh, which is constituted by those two elements that are the portions of the orthogonal fibers that cross a pivot. Each pivot is referenced by two integers n x and m y , which are simply its discrete coordinates along the axes x and y, respectively. We will model each segment of fiber between two consequent nodes as a beam. The four beam elements connected to the pivot {n x , m y } are denoted by [n x -1, n x ] and [n x , n x + 1] for the two elements oriented along x , or by [m y -1, m y ] and [m y , m y + 1] for the two elements oriented along y. The understanding of the reader is help by fig. 1.

The parameters of the beam elements are their Young modulus E j , the area A j of their sections and the moments of inertia I j of their sections, with j = x, y. The dimensions of the sections A j are assumed small with respect to the lengths j . Accordingly, the behavior of the inter-pivot elements can be effectively modeled by the Euler beam model. It is further assumed that the geometrical and the mechanical parameters of both types of beams are of the same order of magnitude: to be more precise and in particular, the period is characterized by the length =

x y in such a way that O( x ) = O( y ) = . We will need in the sequel to introduce also Note that the specificity of the adopted assumption of beam behaviour is that the length of the arrays is not considered (i.e. the fibers are not treated as beams) but instead, this is the topology of the pantographic lattice that leads to a local beam description.

We study in the {x, y} plane the quasi-static small deformations of the lattice in the neighborhood of its initial equilibrium position where the lattice is periodic, having a rectangular period. It is also assumed that the spatial variations of the deformation fields and/or of the placement fields occur at large-scale and have a characteristic variation length L that is large compared to the size of period : in other words to have a relevant variation of the deformation fields the space variables must have an increment larger relatively to .

These considerations naturally introduce the small parameter to be used for specifying scale separation ε = /L 1

The macroscopic description of the system, valid at the dominant order, is reached for ε → 0.

3 Homogenization Method: multiscale asymptotic heuristic approach.

Let us construct the global behavior from the local behavior of beam elements connected by pivots. To do this, the problem is first discretized exactly and then it is converted into a continuum macroscopic formulation by an asymptotic homogenization procedure.

Using the balance laws of an Euler beam, the contact actions (normal and shear forces and moment) at the ends of each beam element are explicitly expressed in terms of the kinematic variables (displacements and rotations) evaluated at the same ends. The balance of each element is thereby assured. We will consider the cases in which the beam elements between the nodes deform in a quasi-static regime. In other words our treatment will be applicable when, in considered phenomena, there is a (quasi-)instantaneous equilibrium, at the level of periodic cell, of the interconnected beams. Sometimes this assumption is referred to as: the assumption of local microscopic instantaneous equilibrium. This approach is relevant when dealing with phenomena of slow time evolution. As a consequence of the exact discretization, to specify the conditions of the global equilibrium conditions, it is necessary and sufficient to focus on the equilibrium of each one of the pivots, which will play the role of material points of the homogenized continuum. It is therefore expressed, in each pivot, the balance of forces and moments applied by the four elements therein connected (belonging to the two orthogonal fibers that intersect at each pivot). This gives an exact representation of the original problem in a discrete form of finite difference type, the variables being the kinematic variables and the actions, of each element, evaluated at the pivot-locations. One has to explicitly remark that in the set of the aforementioned four elements one can distinguish two pairs (parallel in the reference configuration) which are indeed part of the same fiber: moving from one to the other elements belonging to the same fiber, the displacements and rotations are continuous (see figure 1, x-fiber and y-fiber).

Remark also that the action of a pivot on the two continuous fibers which the pivot itself is interconnecting is modeled here as concentrated (in a point) force and couple: of course, by action and reaction principle, the action of the pivot on the fiber of one array is opposite to the action exerted by the same pivot on the fiber of the other array.

The passage from the exact discrete formulation to the macroscopic continuous description, valid at the dominant order, is performed as follows [START_REF] Caillerie | Thin elastic and periodic plates[END_REF]. We assume that the overall behavior of the system can be described by a set of macroscopic fields, the generic element of which is denoted by Q(x, y).

Discrete variables q (nx,my) at the pivots are considered as the values at these points of continuous functions Q(x, y).

This procedure has been described in general already in the works by Piola [START_REF] Dell'isola | The complete works of Gabrio Piola: Volume I -Commented english translation[END_REF] (he seems to us to be among the first scientists having introduced such a heuristic method of homogenization (see [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] Rahali | Homogenization à la Piola produces second gradient continuum models for linear pantographic lattices[END_REF])): however he applied it specifically only to the case of fluids.

We therefore set q (nx,ny) = Q(x n , y m );

x n = n x x , y m = m y y
Consistent with the hypothesis of separation of scale, these continuous functions vary at large scale as O(L). Consequently, the increments of the inter-pivot distance O( ) can be expressed by the Taylor expansions of macroscopic fields. Thus, the finite difference of the discrete formulation is converted into series by introducing successive gradients of the macroscopic functions. As the distances between the pivots are constant due to the periodicity, the terms of such series are of the type, for example,

k x ∂ k Q/∂x k = O( k x /L k )O(Q) = ε k O(Q)
, with a multiplication factor. They therefore involve the orders of magnitude in powers of ε. As we are interested in situations where ε → 0, we need the following specification for the consistency of the orders of magnitude, i.e. we need to represent the functions Q(x, y) in the form of asymptotic expansions of the type:

Q(x, y) = ∞ k=0 ε k Q (k) (x, y); O(Q (k) ) = O(Q (0) )
These asymptotic expansions are to be used in the equations of equilibrium expressed by means of Taylor expansions. The dimensionless parameters that emerge from this formulation must be weighed in power of ε to translate correctly the dominant mechanisms in the studied system [START_REF] Hans | Dynamics of discrete framed structures: a unified homogenized description[END_REF][START_REF] Boutin | Homogenisation of periodic discrete medium: Application to dynamics of framed structures[END_REF]. This normalization of the balance equations preserves the same local physics during the transition to the limit ε → 0. Consequently, the macroscopic model, in the limit, preserves -at the dominant orderthe same local physics than that prevailing in the real system where the scale ratio /L is taking small but finite values.

After aforementioned normalization, we obtain a series of equilibrium conditions in terms of the powers of ε, which may be solved term by term. The macroscopic description, at the dominant order, is constituted by the first nontrivial differential system on the macroscopic variables.

Micro-Modeling of pantographic lattices

In the literature many attention has been paid to lattices of beams interconnected by clamping constraints or to trusses. The reader is referred to [START_REF] Németh | Bielastic web of links: A discrete model of Csonka's beam[END_REF][START_REF] Liew | Improved nonlinear plastic hinge analysis of space frame structures[END_REF][START_REF] Noor | Continuum models for beam-and platelike lattice structures[END_REF] and references there cited. However, the existence of so called "floppy modes" at micro-level in the pantographic structures requires the most attentive consideration. To be more precise: what we call a floppy mode, is a deformation of the micro-structure to which is associated a vanishing energy (for more details see [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] and references there cited). An intuitive consequence of the existence of floppy modes is that for the effective medium the standard condition of coerciveness is not verified and therefore this concept needs to be modified or generalized. Actually one can decompose the space of deformations into a coercive subspace plus the space of floppy modes and the role of boundary condition becomes more determinant in well-posedness problems. These microscopic floppy modes were exploited there to prove that, in general, the class of first gradient continua (those introduced by Cauchy and usually considered in continuum mechanics, see e.g. [START_REF] Dell'isola | Higher-gradient continua: The legacy of Piola, Mindlin, Sedov and Toupin and some future research perspectives[END_REF][START_REF] Dell'isola | The postulations á la D'Alembert and á la Cauchy for higher gradient continuum theories are equivalent: a review of existing results[END_REF]) is not sufficient to describe, at a macro-level, all conceivable physical systems. In particular, when there is a high contrast (see [START_REF] Camar-Eddine | Non-local interactions resulting from the homogenization of a linear diffusive medium[END_REF][START_REF] Pideri | A second gradient material resulting from the homogenization of an heterogeneous linear elastic medium[END_REF][START_REF] Hans | Dynamics of discrete framed structures: a unified homogenized description[END_REF][START_REF] Boutin | Generalized inner bending continua for linear fiber reinforced materials[END_REF][START_REF] Cecchi | Heterogeneous elastic solids: A mixed homogenization-rigidification technique[END_REF]) in the physical properties at micro-level then it may become necessary at macro level to introduce higher order continua (see [START_REF] Dell'isola | How contact interactions may depend on the shape of Cauchy cuts in Nth gradient continua: approach "à la D'Alembert[END_REF][START_REF] Dell'isola | Cauchy tetrahedron argument applied to higher contact interactions[END_REF]).

The analysis which we present now adapts to the case of pantographic lattices the studies presented in [START_REF] Boutin | Generalized beams and continua. dynamics of reticulated structures[END_REF][START_REF] Boutin | Generalized inner bending continua for linear fiber reinforced materials[END_REF][START_REF] Hans | Dynamics of discrete framed structures: a unified homogenized description[END_REF][START_REF] Soubestre | Non-local dynamic behavior of linear fiber reinforced materials[END_REF] to take into account its behavior, which can be regarded as being somehow exotic.

We start by using the framework of local microscopic instantaneous equilibrium to formulate a meso-level model where the lattice is described as a set of nodes (the pivots) interconnected by beams element.

Characterization of the mechanical behavior of a beam element

Let the section of the considered straight (in the reference configuration) beam be A, its moment of inertia being I, and let the material constituting it be elastic, isotropic and with Young's modulus E. We denote (referring them to the beam planar reference configuration) by v the axial displacement, by w the transverse displacement, by N the normal force, by T the shear force and by M the bending moment (we will be using the French convention for the orientation of axes). These contact actions are defined as that of the part s < 0 on part s > 0, where s designates the abscissa along the axis of the beam. Consider a portion of the beam between two points B and C spaced apart by a distance that is large enough compared with the size, i.e. √ A, of the beam section. The Euler beam theory can therefore be used. Accordingly, the rotation of the section is related to the transverse displacement via the relation θ(s) = dw (s) ds . The constitutive equations of the beam are the following:

N (s) = EA dv ds ; M (s) = -EI d 2 w ds 2 (1)
and, in the quasi-static regime, the equilibrium equations in differential form are expressed by:

dN ds = 0; dT ds = 0; dM ds + T = 0 (2)
Suppose that for the point B (respectively C) the displacements and the rotation are v B , w B and θ B , ( respectively v C , w C and θ C ). The forces and moment at B and at C are determined using the equations of beams. They are expressed as a function of the kinematic variables as follows

N B = N C = EA (v B -v C ) (3)
T B = T C = 12EI 3 w B -w C + 2 (θ B + θ C ) (4)
M B = 6EI 2 w B -w C + 3 (2θ B + θ C ) (5)
M C = - 6EI 2 w B -w C + 3 (θ B + 2θ C ) (6)
Expressions of N and T, M involve the axial EA and bending 12EI

R = 12EI 3 EA = 12I 2 A = O( A 2 ) 1 (7)
This strong stiffness contrast plays an essential role in the functioning of the system. Thereafter we take into account explicitly that the aspect ratio of the elements is √ A/ = O(ε) which leads to:

R = 12EI 3 EA = O(ε 2 ) (8)

Discrete kinematic variables and equilibrium at pivots

Because of the operating principle of an internal pivot, the ends of the four elements connected to it undergo the same displacement u x (u y ) along the axis x (y), but rotations of those elements belonging to fibers with distinct orientations are not identical: the coinciding ends of the two elements oriented along x (respectively y) undergo the same rotation θ x (respectively θ y ). Thus, each pivot {n x , m y } is described by four kinematic variables u x(nx,my) , u y(nx,my) , θ x(nx,my) and θ y(nx,my) .

The equilibrium at a pivot results: (i) in the balances of force (exerted on the pivot) along x and y, and (ii) in the continuity of the both moment fields arising in the beam elements oriented either along x or along y.

These four equations expressed at the pivot {n x , m y } take the following form:

-T C [my-1,my] + T B [my,my+1] + N C [nx-1,nx] -N B [nx,nx+1] = 0 Balance of force along x (9) T C [nx-1,nx] -T B [nx,nx+1] + N C [my-1,my] -N B [my,my+1] = 0 Balance of force along y (10) M C [nx-1,nx] -M B [nx,nx+1] = 0
Balance of Moments for the elements along x ( 11)

M C [my-1,my] -M B [my,my+1] = 0
Balance of Moments for the elements along y (

here we have denoted by the symbol Q D [p,q] the value of the field Q at the extremity D = C or D = B of the beam element connecting the node p and the node q , where p and q are two consecutive ones in both x or y direction.

The component u x (respectively u y ) of the pivot displacement is: (i) the axial displacement of the ends of the beam element oriented along x (respectively y), and (ii) the opposite transverse displacement (respectively direct) of the ends of the beam element along y (respectively along x). The change of sign results from different orientations of the global frame and of the local frame of the y-oriented fibers.

Thus, by substituting in equations [START_REF] Auriault | Homogenization of coupled phenomena in heterogenous media[END_REF][START_REF] Beirão Da Veiga | A family of C0 finite elements for Kirchhoff plates II: Numerical results[END_REF][START_REF] Bensoussan | Asymptotic analysis for periodic structures[END_REF][START_REF] Bîrsan | Deformation analysis of functionally graded beams by the direct approach[END_REF] the forces by their expressions in terms of the displacement fields [START_REF] Alibert | Truss modular beams with deformation energy depending on higher displacement gradients[END_REF][START_REF] Allaire | Homogenization and two-scale convergence[END_REF][START_REF] Altenbach | Micropolar shells as two-dimensional generalized continua models[END_REF][START_REF] Andreaus | Identification of multiple open and fatigue cracks in beam-like structures using wavelets on deflection signals[END_REF], the balance of force along x is obtained as follows

12E y I y 3 y (u x(nx,my-1) -2u x(nx,my) + u x(nx,my+1) ) + y 2 (-θ y(nx,my-1) + θ y(nx,my+1) ) + E x A x x u x(nx-1,my) -2u x(nx,my) + u x(nx+1,my) = 0. ( 13
)
The continuity at the nodes of bending moments of the elements oriented along y (after the simplification by

2EyIy y ) reads 3 y (u x(nx,my-1) -u x(nx,my+1) ) -(θ y(nx,my-1) + 4θ y(nx,my) + θ y(nx,my+1) ) = 0, (14)
the balance of force along y reads

12E x I x 3 x (u y(nx-1,my) -2u y(nx,my) + u y(nx+1,my) ) + x 2 (-θ x(nx-1,my) + θ x(nx+1,my) ) + E y A y y u y(nx,my-1) -2u y(nx,my) + u y(nx,my+1) = 0, (15)
and finally the continuity at the nodes of bending moments of the elements oriented along x (after the simplification by 2ExIx

x ) reads

3 x (-u y(nx-1,my) + u y(nx+1,my) ) -(θ x(nx-1,my) + 4θ x(nx,my) + θ x(nx+1,my) ) = 0 (16)
These four equations are split into two independent groups of equations: the first one (13-14) couples the variables u x and θ y ; while the second (15-16) couples the variables u y et θ x . Thus, it is sufficient to treat the equations [START_REF] Boutin | Rayleigh scattering in elastic composite materials[END_REF][START_REF] Boutin | Homogenisation of periodic discrete medium: Application to dynamics of framed structures[END_REF], being the results of equations [START_REF] Boutin | Generalized beams and continua. dynamics of reticulated structures[END_REF][START_REF] Boutin | Generalized inner bending continua for linear fiber reinforced materials[END_REF] easily deduced by changing the roles of axes x and y.

Continuous formulation and asymptotic expansions

Let us introduce continuous kinematic descriptors (denoted by uppercase letters) coinciding with the discrete kinematic variables of the pivots {n x , m y } with coordinates x n = n x and y m = m y :

u x(nx,my) = U x (x n , y m ); θ y(nx,my) = Θ y (x n , y m ) (17)
u y(nx,my) = U y (x n , y m ); θ x(nx,my) = Θ x (x n , y m ) (18)
and use Taylor series expansions to express the terms of finite difference equations [START_REF] Boutin | Rayleigh scattering in elastic composite materials[END_REF][START_REF] Boutin | Homogenisation of periodic discrete medium: Application to dynamics of framed structures[END_REF]. By introducing the dimensionless variables x = x/L and y = y/L, we have (the reader is referred for comparison to [START_REF] Dell'isola | The complete works of Gabrio Piola: Volume I -Commented english translation[END_REF] and [START_REF] Carcaterra | Macroscopic description of microscopically strongly inhomogenous systems: A mathematical basis for the synthesis of higher gradients metamaterials[END_REF]):

u x(nx,my-1) -2u x(nx,my) + u x(nx,my+1) = 2 y ∂ 2 U x ∂y 2 (x n , y m ) + 2 4! 4 y ∂ 4 U x ∂y 4 (x n , y m ) + O( 6 y ∂ 6 U y ∂y 6 ) = ε 2 2 y ∂ 2 U x ∂y 2 + ε 4 4 y 2 4! ∂ 4 U x ∂y 4 + ε 6 6 y 2 6! ∂ 6 U x ∂y 6 + O(ε 8 )(19) u x(nx-1,my) -2u x(nx,my) + u x(nx+1,my) = ε 2 2 x ∂ 2 U x ∂x 2 + ε 4 4 x 2 4! ∂ 4 U x ∂x 4 + ε 6 6 x 2 6! ∂ 6 U x ∂x 6 + O(ε 8 )(20) similarly -θ y(nx,my-1) + θ y(nx,my+1) = 2ε y ∂Θ y ∂y + ε 3 3 y 2 3! ∂ 3 Θ y ∂y 3 + ε 5 5 y 2 5! ∂ 5 Θ y ∂y 5 + O(ε 7 ) (21) -u x(nx,my-1) + u x(nx,my+1) = 2ε y ∂U x ∂y + ε 3 3 y 2 3! ∂ 3 U x ∂y 3 + ε 5 5 y 2 5! ∂ 5 u x ∂y 5 + O(ε 7 ) (22)
and finally

θ y(nx,my-1) + 4θ y(nx,my) + θ y(nx,my+1) = 6Θ y + ε 2 2 y ∂ 2 Θ y ∂y 2 + ε 4 4 y 2 4! ∂ 4 Θ y ∂y 4 + ε 6 6 y 2 6! ∂ 6 Θ y ∂y 6 + O(ε 8 )(23)
By construction, the coefficients of the powers expansion for ε in the equations [START_REF] Carassale | Non-linear discrete models for the stochastic analysis of cables in turbulent wind[END_REF][START_REF] Carcaterra | Macroscopic description of microscopically strongly inhomogenous systems: A mathematical basis for the synthesis of higher gradients metamaterials[END_REF][START_REF]Principes variationnels en fluide compressible et en magnétodynamique des fluides[END_REF][START_REF] Cazzani | Isogeometric analysis of plane-curved beams[END_REF] are of the same the dominant order but also they contain terms of lower order. Therefore, to really order the relative weight of the different addends, it is necessary to introduce now the asymptotic expansions of the variables U x , U y , Θ x and Θ y . This is essential to effectively separate the powers exponents and to ensure the coherence of the passage to the limit ε → 0. It should be noted that consecutive terms of Taylor expansions are systematically offset of ε 2 . It is therefore sufficient to introduce the developments in the even powers of ε. Consequently, we are looking for fields U x , U y , Θ x , Θ y in the following generic form:

U x = U (0) x + ε 2 U (2) x + ε 4 U (4) x + O(ε 6 U (6)
x ) Thereafter, we will denote with a tilde the correction terms that are physically observable, e.g., U (4) x = ε 4 U (4) x ; Θ (4) x = ε 4 Θ (4) x Referring the developments in power of ε 2 in equations (20-23), we get:

u x(nx,my-1) -2u x(nx,my) + u x(nx,my+1) = ε 2 2 y ∂ 2 U (0) x ∂y 2 + ε 4 4 y ∂ 2 U (2) x ∂y 2 + 2 4! ∂ 4 U (0) x ∂y 4 + O(ε 6 ) (24)
u x(nx-1,my) -2u x(nx,my) + u x(nx+1,my) = ε 2 2 x ∂ 2 U (0) x ∂x 2 + ε 4 4 x ∂ 2 U (2) x ∂x 2 + 2 4! ∂ 4 U (0) y ∂x 4 + O(ε 6 ) (25) -θ y(nx,my-1) + θ y(nx,my+1) = ε y 2 ∂Θ (0) y ∂y + ε 3 3 y 2 ∂Θ (2) y ∂y + 2 3! ∂ 3 Θ (0) y ∂y 3 + ε 5 5 y 2 ∂Θ (4) y ∂y + 2 3! ∂ 3 Θ (0) y ∂y 3 + 2 5! ∂ 5 Θ (0) y ∂y 5 + O(ε 7 ) (26) -u x(nx,my-1) + u x(nx,my+1) = ε y 2 ∂U (0) x ∂y + ε 3 3 y 2 ∂U (2) x ∂y + 2 3! ∂ 3 U (0) x ∂y 3 + ε 5 5 y 2 ∂U (4) x ∂y + 2 3! ∂ 3 U (2) x ∂y 3 + 2 5! ∂ 5 U (0) x ∂y 5 + O(ε 7 ) (27)
θ y(nx,my-1) + 4θ y(nx,my) + θ y(nx,my-1) = 6Θ (0)

y + ε 2 2 y Θ (2) y + ∂ 2 Θ (0) y ∂y 2 + ε 4 4 y Θ (4) y + ∂ 2 Θ (2) y ∂y 2 + 2 4! ∂ 4 Θ (0) y ∂y 4 + O(ε 6 ) (28)
The calculations presented in this subsection provides an accurate transformation of the finite differences into successive derivatives. This step is essential to get the continuous asymptotic model valid in the limit ε → 0.

Asymptotic macroscopic model

In the expressions (24-28) do appear the macroscopic continuous fields and their macroscopic derivatives. By substituting them in the equilibrium equations [START_REF] Boutin | Rayleigh scattering in elastic composite materials[END_REF][START_REF] Boutin | Homogenisation of periodic discrete medium: Application to dynamics of framed structures[END_REF] we may obtain a macroscopic continuous formulation of the equilibrium of pivots (equilibrium of force along x and equilibrium of moment for the elements along y). To make explicit which are the appearing powers of ε, it is convenient to write the obtained equations in the non-dimensional variables x and y , where L is the reference length. As, by hypothesis, x and y are of order ε with respect to L, we have

x = x = ε x L et y = ε y L.
∂U (0) x ∂y + L y Θ (0) y + O(ε 6 ) = 0 ( 29
)
and the equation ( 13) of the equilibrium of forces along x gives:

R x ε 2 2 y ∂ ∂y ∂U (0) x ∂y + L y Θ (0) y + ε 2 2 y ∂ ∂y ∂U (2) x ∂y + L y Θ (2) y + ∂ 3 ∂y 3 2 4! ∂U (0) x ∂y + 1 3! L y Θ (0) y + O(ε 4 ) +ε 2 2 x ∂ 2 U (0) x ∂x 2 + ε 2 2 x ∂ 2 U (2) x ∂x 2 + 2 4! ∂ 4 U (0) x ∂x 4 + ε 4 4 x ∂ 2 U (4) x ∂x 2 + 2 4! ∂ 4 U (2) x ∂x 4 + 2 6! ∂ 6 U (0) x ∂x 6 + O(ε 6 ) = 0 (30)
where

R x = 12E y I y 3 y x E x A x
The moment equation ( 29) comes in the form of a series in which the mechanical characteristics of beams do not interfere. The convergence of the series when ε → 0 implies that each involved term vanishes. Consequently:

∂U (0) x ∂y + L y Θ (0) y = 0 (31) ∂U (2) x ∂y + L y Θ (2) y = 0 (32) ∂U (4) x ∂y + L y Θ (4) y + 1 3 • 4! ∂ 4 ∂y 4 1 5 ∂U (0) x ∂y -L y Θ (0) y = 0 (33)
This precisely means that, the equilibrium of moments for y fibers requires, with an error being equal to O(ε 4 ), a relationship between their rotation and their transverse gradient, which is expressed in the dimensional fields as follows:

∂U x ∂y + y Θ y = O(ε 4 )
At the dominant order, Θ y can thus be considered a hidden variable which does not emerge in the macroscopic description at the leading order. Only by considering the correction O(ε 4 ), the rotation differs from the transverse gradient, because we have

∂U (4) x ∂y + L y Θ (4) y = - 2 5! ∂ 5 ∂U (0) x ∂y 5
or, in dimensional variables and denoting explicitly the observable corrections U (4)

x = ε 4 U (4)
x ; Θ

x = ε 4 Θ (4) x : ∂ U (4) x ∂y + y Θ (4) y = -4 2 5! ∂ 5 ∂U (0) x ∂y 5 (4)
Let us replace now the results [START_REF] Della Corte | Referential description of the evolution of a 2D swarm of robots interacting with the closer neighbors: Perspectives of continuum modeling via higher gradient continua[END_REF][START_REF] Dell'isola | At the origins and in the vanguard of peridynamics, nonlocal and higher-gradient continuum mechanics: An underestimated and still topical contribution of Gabrio Piola[END_REF] in the balance equation [START_REF] Vescovo | Dynamic problems for metamaterials: review of existing models and ideas for further research[END_REF]. This leads to

R x ε 2 4 y 2 4! ∂ 4 U (0) x ∂y 4 + O(ε 4 ) = 2 x ∂ 2 U (0) x ∂x 2 + ε 2 2 x ∂ 2 U (2) x ∂x 2 + 2 4! ∂ 4 U (0) x ∂x 4 + O(ε 4 ) (34)
To exploit this equation it is necessary to weigh the effects of bending (left-hand side term) and extension (right-hand side term). These effects are a consequence of both the mechanical properties of considered system and of the nature of admitted kinematics. The mechanical parameters which we choose will introduce the high contrast condition (8) which is expressed by

R x = R x ε 2 . ( 35
)
Regarding the nature of the macroscopic kinematics, we are led thus to distinguish between low or high contrast situations in the axial and the transverse gradient.

Low contrast between axial and transverse gradient of U x

We consider here macroscopic kinematics where axial and transverse components of the displacement gradient U x are of the same order, i.e.,

∂U (0) x ∂x = O ∂U (0) x ∂y (36)
This estimate explicitly means that the axial and transverse variations have as a common evolution characteristic value O(L). This hypothesis is usually considered in the case of an elastic composite medium where the terms of the strain tensor components are assumed of the same order. In this case, we obtain successively the equation (each relative to the orders ε 0 , ε 2 , ε 4 ):

∂ 2 U (0) x ∂x 2 = 0 ∂ 2 U (2) x ∂x 2 + 2 4! ∂ 4 U (0) x ∂x 4 = 0 ∂ 2 U (4) x ∂x 2 + 2 4! ∂ 4 U (2) x ∂x 4 + 2 6! ∂ 6 U (0) x ∂x 6 = ( x ) -2 R x 4 y 2 4! ∂ 4 U (0) x ∂y 4
By simplifying and returning to the dimensional variables and observable correctors we deduce that:

E x A x y ∂ 2 U (0) x ∂x 2 = 0 (37)
E x A x y ∂ 2 U (2) x ∂x 2 = 0 (38) E x A x y ∂ 2 U (4) x ∂x 2 = E y I y x ∂ 4 U (0) x ∂y 4 (39)
The equations [START_REF] Dell'isola | Designing a light fabric metamaterial being highly macroscopically tough under directional extension: first experimental evidence[END_REF][START_REF] Dell'isola | Boundary conditions at fluid-permeable interfaces in porous media: A variational approach[END_REF] mean that the tension of the fibers oriented along x is constant to the accuracy ε 4 . Only by considering the order 4 of the correctors, the tension of the beams varies due to the bending of orthogonal beams, as indicated by [START_REF] Dell'isola | Cauchy tetrahedron argument applied to higher contact interactions[END_REF].

The assumption [START_REF] Dell'isola | Large deformations of planar extensible beams and pantographic lattices: Heuristic homogenisation, experimental and numerical examples of equilibrium[END_REF] obviously can not cover all cases of loading. In particular, it is not predictive if the lattice is subjected to an uniaxial extension in a direction that does not coincide with one of the directions of the fibers' arrays (see Fig. 2). Indeed if one considers the "red" fiber in Fig. 2 (which presents an experimental evidence), it is clear that its state of tension cannot be constant: its tension is not vanishing in the clamping but its clearly vanishes at the free end. To describe these situations it is necessary to change the too restrictive hypothesis [START_REF] Dell'isola | Large deformations of planar extensible beams and pantographic lattices: Heuristic homogenisation, experimental and numerical examples of equilibrium[END_REF] by allowing a strong contrast between the axial and transverse components of the gradient of the macroscopic displacement.

Strong contrast between axial and transverse gradient of U x

Thus we are lead to consider macroscopic kinematics where the transverse gradient ∂U (0)

x ∂y dominates in comparison to the axial gradient ∂U (0)

x ∂x , i.e.

∂U (0) x ∂x ∂U (0) x ∂y (40)
This assumption reflects the fact that the axial characteristic length of variation of U x (denoted L x,a ) is much larger than the transverse (denoted L x,t ) length (see e.g. experimental evidence presented in [START_REF] Dell'isola | Large deformations of planar extensible beams and pantographic lattices: Heuristic homogenisation, experimental and numerical examples of equilibrium[END_REF]). This is a consequence of the high shear deformability of pantographic network, compared to its deformability in the axis of the fibers.

To assure scale separation we are thus lead to consider as reference length the smaller between the two i.e. L = L x,t so that

∂ i U (0) x ∂x i = O(U (0) x ) Lx,a L i = O(U (0) x ) L x,t L x,a i ; ∂ j U (0) x ∂y j = O(U (0) x ) Lx,t L j = O(U (0) x )
Once we consider L x,t = ε 2 L x,a we have

∂U (0) x ∂x = ε 2 O ∂U (0) x ∂y ; (41)
∂ 2 U (0) x ∂x 2 = ε 4 O ∂ 4 U (0) x ∂y 4 (42)
The presence of a contrast in the components of the strain tensor as given by ( 41) is unconventional in elastic composite medium but arises naturally in the case of weakly compressible viscous fluid (where the trace of the strain rate tensor is negligible compared to its deviatoric component) or in beams and plates (where deformations in the section of the beam-or in the thickness of the plate-are negligible). This contrast is present also in pantographic sheets.

Replacing the estimate ( 42) into [START_REF] Dell'isola | Pantographic 2D sheets: Discussion of some numerical investigations and potential applications[END_REF] we get at the leading order

∂ 2 U (0) x ∂x 2 = ( x ) -2 R x 4 y 2 4! ∂ 4 U (0) x ∂y 4
or, returning to the dimensional variables and normalizing by introducing the surface of the periodic cell we have:

E x A x y ∂ 2 U (0) x ∂x 2 = E y I y x ∂ 4 U (0) x ∂y 4 (43)
This equation indicates that the normal force (on left-hand side) varies at the first order in the beams due to the shear force exerted by the orthogonal beams (on the right-hand side). This is made possible because the transverse gradient is of two orders of magnitude higher than that of the extension gradient (see again [START_REF] Dell'isola | Variational principles are a powerful tool also for formulating field theories, volume 535 of Chapter in[END_REF]). Moreover, we note that the equation ( 43) is more general than the equation [START_REF] Dell'isola | Designing a light fabric metamaterial being highly macroscopically tough under directional extension: first experimental evidence[END_REF] and it is needed in the considered mechanical system. Moreover equation [START_REF] Dell'isola | How contact interactions may depend on the shape of Cauchy cuts in Nth gradient continua: approach "à la D'Alembert[END_REF] degenerates to the equation [START_REF] Dell'isola | Designing a light fabric metamaterial being highly macroscopically tough under directional extension: first experimental evidence[END_REF] when

∂U (0) x ∂x ε 2 ∂U (0)
x ∂y . We will use in what follows the description [START_REF] Dell'isola | How contact interactions may depend on the shape of Cauchy cuts in Nth gradient continua: approach "à la D'Alembert[END_REF] which applies to more general kinematics.

Synthesis of obtained results

The above results, derived from the system [START_REF] Boutin | Rayleigh scattering in elastic composite materials[END_REF][START_REF] Boutin | Homogenisation of periodic discrete medium: Application to dynamics of framed structures[END_REF] for the variables U (0)

x and Θ (0) y is transposed by a similar analysis of the system [START_REF] Boutin | Generalized beams and continua. dynamics of reticulated structures[END_REF][START_REF] Boutin | Generalized inner bending continua for linear fiber reinforced materials[END_REF] to the variables U (0) y and Θ (0)

x . The description in small deformations of the orthogonal pantographic lattice is therefore obtained by restricting the analysis at the dominant order (for simplicity, the exponents of order (0) are removed)

E x A x y ∂ 2 U x ∂x 2 = E y I y x ∂ 4 U x ∂y 4 (44)
E y A y x ∂ 2 U y ∂y 2 = E x I x y ∂ 4 U y ∂x 4 (45)
with the addition of the relationships between the transverse gradients and rotations

∂U x ∂y + y x Θ y = 0; ∂U y ∂x - x y Θ x = 0 (46)
Remark that equations (44-45) can be rewritten by introducing two intrinsic characteristic lengths η x and η y as follows

∂ 2 U x ∂x 2 = η 2 x ∂ 4 U x ∂y 4 ; ∂ 2 U y ∂y 2 = η 2 y ∂ 4 U y ∂x 4
where η 2

x :=

y E y I y x E x A x ; η 2 y := x E x I x y E y A y .
It is clear that the physics of the system is governed by these internal intrinsic lengths which differ from the size of the cell.

The model governed by [START_REF] Dell'isola | A two-dimensional gradient-elasticity theory for woven fabrics[END_REF][START_REF] Dieudonné | Pour l'honneur de l'esprit humain[END_REF] is of the type of a conservative generalized continuum medium having deformation energy depending on first and second order gradient of displacement (see [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]). Its evolution ruled by two displacement fields which are independent and uncoupled U x and U y . This particular behavior is within the framework of the second gradient continua where the internal actions are described by a stress symmetric tensor T of order two and by a hyper stress tensor T of third order. To make the identification, we pose [START_REF] Dell'isola | A two-dimensional gradient-elasticity theory for woven fabrics[END_REF][START_REF] Dieudonné | Pour l'honneur de l'esprit humain[END_REF] in the following form where the differential operator DIV denotes the Lagrangian divergency:

DIV(T -DIV(T )) = 0 ( 47
)
and where we used the definitions T := C.(∇U + ∇U t )/2, T := D.∇∇U [START_REF] Federico | Elasticity and permeability of porous fibre-reinforced materials under large deformations[END_REF] in which the symbol '.' denote the repeated index saturation between different order tensors and ' t ' the transposition of second order tensors. The elasticity tensors C and D of order four and six respectively have the particular form defined by their components as follows:

C abcd := E x A x y δ ax δ bx δ cx δ dx + E y A y x
δ ay δ by δ cy δ dy (49)

D abcdef := E y I y x δ ax δ by δ cy δ dx δ ey δ f y + E x I x y δ ay δ bx δ cx δ dy δ ex δ f x (50)
Remark that (separately) the tensors C and D are not coercive. For instance, pure shear deformations do not have any first gradient energetic content, while any affine displacement field has not any second gradient energetic content. This is consistent with the existence, for pantographic sheets, of floppy modes as nullifiers of deformations energy, in addition to the standard rigid motions. Nevertheless, considering both first and second gradient energies together, the system can be qualified as "subcoercive" when suitable boundary conditions are imposed in such a way that floppy modes are excluded.

In the framework of this requirement, the set of admissible displacements is restricted if compared to the one needed in first gradient theory, where only rigid body motions are excluded. We conjecture that in the set of considered admissible displacements the total deformation energy is definite positive and leads to well-posed problems. This seems physically well-grounded and is confirmed by all the performed numerical simulations presented in the last section.

To identify the class of physically meaningful boundary conditions, in the next section we more closely study the structure of the considered energy.

Energy Formulation and boundary conditions for pantographic lattices

The general framework of second gradient continua (as formulated e.g. by [START_REF] Dell'isola | The postulations á la D'Alembert and á la Cauchy for higher gradient continuum theories are equivalent: a review of existing results[END_REF]) enables, by means of the energy formulation and of the divergence theorem, to specify the boundary conditions which can be consistently considered as applicable to them.

Let B be a bi-dimensional medium whose elastic energy depends among the displacement gradient ∇U acting on the second order stress tensor, T, and among the second displacement gradient acting on the third order hyperstress tensor, T . In this case the energy W of B is:

2W = ˆB(T : ∇U + T . . .∇∇U) (51)
Transforming this expression through successive integrations by parts enables to make the link between the variation of internal energy and the energy supplied at the boundary of B (see e.g. [START_REF] Dell'isola | The postulations á la D'Alembert and á la Cauchy for higher gradient continuum theories are equivalent: a review of existing results[END_REF]). As the considered system is bi-dimensional, the boundary consists in a set of regular edges ∂B with unit normal n and wedges ∂∂B consisting in the union of a finite number N of vertices S I , i.e. ∂∂B = ∪S I = {S} . For our calculations we introduce Levi-Civita tensorial notation, keeping trace of covariance and contravariance nature of considered tensorial quantities and we use latin indices for Eulerian components and greek letters for Lagrangian components.

δW = ´B(T : ∇δU + T . . .∇∇δU) = ´B T β i δU i ,β + T αβ i δU i ,αβ = ´B -T α i,α δU i -T αβ i,β δU i ,α + ´B T β i δU i + T αβ i δU i ,α ,β = ´B -T α i,α δU i -T αβ i,β δU i ,α + ´∂B T β i δU i + T αβ i δU i ,α n β = ´B T αβ i,βα -T α i,α δU i + ´∂B (T β i -T βα i,α )δU i + T αβ i δU i ,α n β (52)
Using the equilibrium equation ( 47) the first integral in the RHS vanishes. Furthermore, the last term of the second integral can be decomposed into a tangent and a normal contribution. As the considered system is bi-dimensional, the boundary ∂B includes a set of regular edges. Thus, introducing the projector on the tangent direction of ∂B, P = I -n ⊗ n (hence P γ α = δ γ α -n α n γ and P.P = P ), we have:

δW = ´∂B T β i -T βα i,α n β δU i + T αβ i n α n β δU i ,γ n γ + T αβ i n β P δ α δU i ,γ P γ δ = ´∂B T β i -T βα i,α n β δU i + T αβ i n α n β δU i ,γ n γ + T αβ i n b P δ α δU i ,γ P γ δ -T αβ i n β P δ α ,γ P γ δ δU i = ´∂B T β i -T βα i,α n β -T αβ i n β P δ α ,γ P γ δ δU i + ´∂B T αβ i n α n β δU i ,γ n γ + ´∂∂B T αβ i n β P d α N δ δU i (53)
On the last integral on ∂∂B = {S} the vectors n and N takes the values n ± and N ± defined on the both sides of the discontinuity of the edges. This is also the case for T and P . On each side the vector N is the tangent vector to ∂B which is the outer pointing normal to the border of ∂B, see Figure 4. Thus P + .N + = N + and P -.N -= N -. In a condensed way, we can write (where the symbol refers to the tangent projection of a tensor on ∂B)

´B(T : ∇δU + T . . .∇∇δU) = ´∂B (T -DIV T ) .n -DIV T .δU + (T .n.n) δ dU dn + {S} [T .n.N ].δU (54) with [T .n.N ] = T + .n + .N + + T -.n -.N -.
The previous integration by parts implies that in order to have well-posed problems one has to suitably assign boundary conditions in a specific way (see e.g. [START_REF] Mindlin | Micro-structure in linear elasticity[END_REF]). In particular on any regular point belonging to the set of edges ∂B one can assign:

• either a component U i or its dual quantity T b i -T βα i,α n β -T αβ i n β P δ α ,γ P γ δ .
In the case of rectilinear edge ∂B the normal and the projector are constants and this expression simplifies into

T β i -T βα i,α -T αβ i,γ (δ γ α -n γ n α ) n β .
• either a component of displacement normal gradient, dU i /dn, or its dual quantity T αβ i n a n β .

Moreover on vertices constituting ∂∂B = {S} one can assign:

• either a component U i or its dual quantity [T αβ i n b N α ] which expresses the discontinuity of the wedge on the vertices. (see Fig. 4) Therefore the weak condition for equilibrium is given by

δW = ˆ∂B f ext i δU i + τ ext i δU i ,γ n γ + ˆ∂∂B F ext i δU i (55)
where external actions are given by forces per unit line f ext i but also by concentrated forces on vertices F ext i and by double forces per unit line τ ext i (see [START_REF] Germain | The method of virtual power in continuum mechanics. Part 2: Microstructure[END_REF]). One can apply these results to the particular pantographic lattice considered, where the preferred orthogonal directions x and y coincide with the array of material fibers. In the context of small deformations, the Lagrangian and Eulerian coordinates systems can be identify. In the pantographic case given by ( 49) and ( 50) the deformation energy is given by

2W = ˆB E x A x y (U x ,x ) 2 + E y A y x (U y ,y ) 2 + E y I y x (U x ,yy ) 2 + E x I x y (U y ,xx ) 2 dxdy ( 56
)
so that we have that, as T = ∂W ∂∇U and T = ∂W ∂∇∇U , the only nonvanishing components of stress and hyperstress are:

T x x = E x A x y U x ,x ; T y y = E y A y x U y ,y T xx y = E x I x y U y ,xx ; T yy x = E y I y x U x ,yy .
Thus, for the pantographic sheet, the two types of kinematic and static (more often called natural) boundary conditions which apply to a straight line edge of normal n = n x e x + n y e y take the following form: On the vertices, noting that n and N are orthogonal, the flux dual to the displacement components U x and U y are respectively (T -yy

U x dual of T x x n x + -T yy x,y (1 + n 2 x ) + T yy x,
x n - x n - y -T +yy x n + x n + y ) and (T -xx y n - x n - y -T +xx y n + x n + y ).
Remark that they vanish when the vertex angle is of π/2 and the boundary are parallel to the fibers.

For a better physical insight consider for instance a straight edge oriented along the direction x, then n = e y and the boundary condition simplifies into

U x dual of -T yy
x,y = -

E y I y x U x ,yyy ; dU x /dy dual of T yy x = E y I y x U x ,yy U y dual of T y y = E y A y x U y ,y ; dU y /dy dual of T xx y n 2 x = 0
These boundary conditions on the lattice can be easily understood by recalling standard beam theory. They show that:

• a virtual displacement tangent to the fiber materializing the edge develops energy due to the shear forces in the orthogonal fibers,

• a virtual displacement normal to the edge fiber develops energy due to the tension forces in the orthogonal fibers,

• a rotation of the edge fiber develops energy due to the couple in the orthogonal fibers,

• no energy is developed (in small deformations) by the extension of the orthogonal fibers associated with the couple of the edge fiber.

If the vertex at the end of the edge along x presents an internal angle α, the static quantities dual to the components U x and U y are respectively the weighted couples T +yy x sin(2α)/2 and -T +xx y sin(2α)/2. No energy is developed if the vertex angle is of π/2 with edges oriented along the fibers.

Similar interpretation apply when the straight edge presents an angle with the fiber orientation, however due to the coupling between the forces and couples in different directions, the physical interpretation of the boundary condition becomes more difficult.

In this paper, for the sake of simplicity, we consider only imposed boundary conditions on displacements and displacement gradients and we consider weak form (55) of equilibrium conditions, so that no dual boundary conditions (on forces or double forces) are assigned. In further papers we will consider more general situations.

6 Some equilibrium shapes of linear pantographic sheets: numerical simulations Pantographic sheets have an exotic behavior which is not only characterized by their anisotropy as evidenced by (i) their vanishing resistance to shear deformation, and (ii) their significant resistance to elongation along fibers; but also by their capacity to resist to variations in their so-called "geodesic curvature" (see [START_REF] Dell'isola | A two-dimensional gradient-elasticity theory for woven fabrics[END_REF][START_REF] Steigmann | Mechanical response of fabric sheets to three-dimensional bending, twisting, and stretching[END_REF][START_REF] Giorgio | Pattern formation in the threedimensional deformations of fibered sheets[END_REF][START_REF] Giorgio | Buckling modes in pantographic lattices[END_REF]) i.e. the changes of curvature of material curves induced by inplane displacements. Moreover in their deformation patterns one can observe also the onset of inner boundary layers where bending of constituting beams is concentrated, as suggested by the existence of the intrinsic characteristic lengths η x and η y . The aim of this section is to provides numerical illustrations of the theoretical developments presented in the previous sections. The considered examples may seem purely academic or dictated simply by the taste of investigating mathematical structures (see [START_REF] Dieudonné | Pour l'honneur de l'esprit humain[END_REF]). Although we indeed consider that scientific knowledge is based on the study of exemplary cases (see [START_REF] Russo | The forgotten revolution: how science was born in 300 BC and why it had to be reborn[END_REF][START_REF] Woodcroft | The pneumatics of Hero of Alexandria. The Pneumatics of Hero of Alexandria Publisher: TAYLOR WALTON AND MABERLY[END_REF][START_REF] Heath | A History of Greek Mathematics[END_REF][START_REF] Heath | The Method of Archimedes, recently discovered by Heiberg: A supplement to the Works of Archimedes[END_REF][START_REF] Archimedes | On Floating Bodies[END_REF]). In addition, a potential application of the presented results concern the forming of fiber reinforced composites (see e.g. [START_REF] Cao | Characterization of mechanical behavior of woven fabrics: experimental methods and benchmark results[END_REF][START_REF] Launay | Experimental analysis of the influence of tensions on in plane shear behaviour of woven composite reinforcements[END_REF][START_REF] D'agostino | Continuum and discrete models for structures including (quasi-) inextensible elasticae with a view to the design and modeling of composite reinforcements[END_REF][START_REF] Harrison | Modelling the forming mechanics of engineering fabrics using a mutually constrained pantographic beam and membrane mesh[END_REF][START_REF] Abdiwi | Modelling the shear-tension coupling of woven engineering fabrics[END_REF][START_REF] Nikopour | Concentrated loading of a fibre-reinforced composite plate: Experimental and computational modeling of boundary fixity[END_REF]).

All the presented numerical simulations are obtained by a code elaborated using COMSOL Multiphysics. The homogenized energy introduced in this paper ( 56) is minimized by using the package 'Weak Form PDE' and by introducing standard third order Hermite finite elements. While the used code is surely not optimized for the introduced problem (we believe that the recently developed numerical methods would be more efficient, see e.g. [START_REF] Cazzani | Isogeometric analysis of plane-curved beams[END_REF][START_REF] Cazzani | Constitutive models for strongly curved beams in the frame of isogeometric analysis[END_REF][START_REF] Greco | B-Spline interpolation of Kirchhoff-Love space rods[END_REF][START_REF] Greco | An implicit G1 multi patch B-spline interpolation for Kirchhoff-Love space rod[END_REF][START_REF] Greco | An isogeometric implicit G1 mixed finite element for Kirchhoff space rods[END_REF][START_REF] Turco | A three-dimensional B-spline boundary element[END_REF][START_REF] Greco | Consistent tangent operator for an exact Kirchhoff rod model[END_REF][START_REF] Beirão Da Veiga | A family of C0 finite elements for Kirchhoff plates II: Numerical results[END_REF][START_REF] Della Corte | Referential description of the evolution of a 2D swarm of robots interacting with the closer neighbors: Perspectives of continuum modeling via higher gradient continua[END_REF]) its rate of convergence seems satisfactory for getting preliminary results concerning the behavior of the simplest structures: actually it is based on the introduction of an auxiliary tensor field which appears in the deformation energy and which is equated to the displacement gradient by means of suitable fields of Lagrange Multipliers. Remark also that all presented numerical simulations are really and intrinsically meshindependent, because of the properties of introduced continuum model, where the second gradient of displacement is at the same time modeling the relevant physical properties and supplies a regularizing effect on equilibrium equations.

In the presented simulations we have chosen a lattice made of square cells, so that x = y = and we have imposed that the x and y fibers have identical rectangular sections (having sides a and b) and elastic moduli so that E y = E x , I x = I y = ba 3 12 and A x = A y = ab. As a consequence we have that

η 2 = η 2 x = η 2 y = I x A x = I y A y = a 2 12 .
The numerical values a = 0.9 mm and b = 1.6 mm are that used for the pantographic structures used in the experimental measurements (see [START_REF] Dell'isola | Designing a light fabric metamaterial being highly macroscopically tough under directional extension: first experimental evidence[END_REF]) having rectangular sections. Young's modulus is 1600 MPa. We remark that the Elastica model for beams is applicable in the considered situation as = 4.95 mm and the number of cells is sufficiently large (L = 42.42 × ) to apply the homogenized model.

In the following subsections, we present the numerical simulation of bias tests in different configurations. First, we consider rectangular specimens undergoing standard bias test in extension, but In all the following figures the black lines indicate the local actual orientation of the material fibers (which are orthogonal in the reference configuration), and the deformed shape is displayed together with the map of the stored energy density. All the calculations are performed in the framework of small deformations, however, for a better insight, the deformation is magnified in the figures.

Bias tests on rectangular pantographic sheet

We start by considering a pantographic sheet having a rectangular initial shape with the long side three times longer than the short one.

In the standard extension bias test the short sides are rigidly displaced in the direction of the long side. The Figure 5a shows the small deformation pattern. Remark that the deformation energy is concentrated along the material lines constituted by some fibers and that, meanwhile the fibers are extensible in the present case, the structure of the deformed shape is similar to the one described by Pipkin [START_REF] Pipkin | Some developments in the theory of inextensible networks[END_REF][START_REF] Pipkin | Plane traction problems for inextensible networks[END_REF] in the case of inextensible fibers. Notice also that the distribution of deformation is strongly non homogeneous while respecting the symmetry of the loading. The simulation enables to identify several zones with different kinematics. The partitioning of the specimen can be describe as follows:

-The clamping zones that consists in two "isosceles triangles" having as basis the short sides.

These areas do not sustain any significant deformation, neither in extension nor in shear, so that the first gradient description (37) applies. The vanishing stored energy results here from the quasi-rigid body motion of these regions;

-The deformed zones outside of the aforementioned triangles, in which large shear deformation arise. In this highly sheared domain, the behavior is governed by the second gradient description [START_REF] Dell'isola | A two-dimensional gradient-elasticity theory for woven fabrics[END_REF][START_REF] Dieudonné | Pour l'honneur de l'esprit humain[END_REF]. One distinguishes a central and four lateral sub-domains of vanishing stored energy delimited by transition zones. The different sub-domains correspond to the occurence of local "floppy" modes. Indeed, because of the specific boundary condition, "floppy" modes on the whole specimen are forbidden. Nevertheless, the minimum of energy is attained by activating local floppy modes, far from the boundary conditions. This results in large parts of the body where the deformation energy is very close to be vanishing;

-The transition zones between the different domains (of quasi-floppy modes or quasi-monolithic type) consist in inner boundary layers where the bending of the fibers is concentrated to accommodate the different kinematics that prevails in the two regions in contact. Such layers, that concentrate the elastic energy, take place along material lines constituted by fibers and are characterized by large gradient of shear deformation. These specific features are the direct signature of the second gradient effects: a standard Cauchy continuum description (i.e. simple gradient description as [START_REF] Dell'isola | Designing a light fabric metamaterial being highly macroscopically tough under directional extension: first experimental evidence[END_REF]) would lead to an homogeneous deformation pattern and avoid the development of shear bands.

Furthermore, the qualitative and quantitative comparison of the numerical simulations plotted on Figure 5a and the experimental data presented in ( [START_REF] Dell'isola | Designing a light fabric metamaterial being highly macroscopically tough under directional extension: first experimental evidence[END_REF]) clearly argue in favor of the second gradient description [START_REF] Dell'isola | A two-dimensional gradient-elasticity theory for woven fabrics[END_REF][START_REF] Dieudonné | Pour l'honneur de l'esprit humain[END_REF] for two reasons. First, the experimental deformation pattern is obviously non homogeneous with inner boundary layers whose the structure presents a similar geometry as that observed in the simulation. Second, by making a best fit of the experimental data ( [START_REF] Dell'isola | Designing a light fabric metamaterial being highly macroscopically tough under directional extension: first experimental evidence[END_REF]), the effective parameters of the second gradient continuum model has been identified. It happens that this "blindy" procedure supply exactly the same values of the effective parameters than the one calculated from the micro-macro upscaling procedure, once the geometric and mechanic properties of the cell beams of the sheet tested experimentally are taken.

In addition simulations of shear bias tests with uniform lateral displacement imposed on the top side, and bending bias tests with rotational displacement imposed on the top side (corresponding to a rigid body rotation centered in the middle of the specimen) have been done. The results are displayed on Figures 5b and5c. Similar general comments as done for the extension test still apply, however the geometry of the shear bands and the energy distribution is modified. In particular, the partitioning of the specimen submitted to bending shows only one sub-domain of quasi-floppy mode which is confined between two quasi-monolithic zones.

Bias test on circular pantographic sheets

Consider now the same bias test as the previous one, except that the rectangular pantographic sheet is replaced by a sheet of circular initial shape. The imposed deformation are obtained by clamping one circular arc and impose a rigid displacement of the opposed one.

In fig. 6a) we consider the extension imposed by a relative rigid translation of the two arcs in the direction of the common bisecting diameter. Note the great similarity of the deformation pattern obtained with rectangular and circular sheets, when focusing on the internal rectangular domain considered in the standard bias test. This means that, independently of rectangular or circular geometry of the sheet, the structure of the deformation pattern is kept (almost) unchanged when identical kinematic boundary conditions are imposed, Remark also that the largest lateral dimension of the circular sheet enable the inner boundary layers to extend and intersect. Thus the partitioning of the specimen is complemented by the appearance of two additional lateral sub-domains.

These observations indicate that in the considered cases, the structure of the inner boundary layer (hence of the partitioning) mostly results from the geometry of the boundary conditions, while the extension of the shear bands depends on the geometry of the whole body. Furthermore, in the circular case, the activation of the deformation outside of the rectangle implies additional energy. Therefore, to reach the same displacement imposed at the boundaries, a largest force is required in the circular case than in the rectangular case .

In fig. 6b) the relative rigid displacement is in the direction orthogonal to the bisecting diameter and imposes a shear to the specimen. The resulting shear deformation pattern is similar to the one activated in extension 6a) but with different orientations of the inner boundary layers. The comparison with the case of rectangular sheet lead to the similar comments as above.

In fig. 6c) the relative displacement of the two arcs is obtained by fixing one arc and by rotating the second one with respect the center of the circle. Remark that the partitioning of the specimen is similar as in the rectangular case but here the boundary layers invade the whole body.

Bias test on initially circular pantographic sheets with central holes

In the same body considered in the previous subsection a square hole (of 14×14 cells) is now carved in its central part. The initial orientation of the hole relatively to the fibers is varied from zero (i.e. the sides of the square are along the fibers) to π/8 and π/4 (i.e. the edges of the square are along and orthogonal to the extension displacement). These different cases enable to investigate the effect of the hole on the deformation pattern and on the onset of inner boundary layers.

One may expected, that when a hole is carved in a sub-domain corresponding to quasi-floppy modes, its influence should be negligible since in both case, the hole or floppy modes deformation occurs with no energy expense. However, if a hole intersects the deformation boundary layers in which the energy would be localized in absence of hole, then the energy distribution is necessarily modified and so does the deformation patterns. The colors indicate the levels of stored energy density.

q 0 q 0 q 0 This is what is observed in the simulations. The edges of the carved hole are such that when oriented along the fibers, the hole almost belong to the central floppy mode sub-domain that arises in the non-carved sheet, while when rotated, two corners of the hole may cross the boundary layers of the intact sheet. In extension tests presented on Figs. 7a and 7c the intersection of the carved hole with the deformation boundary layers of the intact body is relatively small and the deformation patterns are weakly disturbed. This is not the case in Figs. 7b where the hole cut the material lines i.e. fibers in which the energy would be concentrated in intact body. As a consequence the deformation response and the partitioning of the specimen change significantly. The same trends are observed for the bending test: weak perturbations arise in Fig. 8a, but large and dramatic change appears in Fig. 8b and Fig. 8c respectively.

Conclusions

Pantographic sheets belong to the specific class of architectured materials whose mechanical behavior is characterized by:

• very high contrast of extensional and bending stiffness at microscopic level (see equation ( 35));

• very high contrasted gradients of displacement in the axial and transverse direction at the macro level (see equation ( 40));

• a microstructure which produces a discretely oriented, orthotropic material exhibiting an extreme anisotropy leading to the presence of two preferred material directions having very high extensional stiffness; this circumstance causes the onset of internal boundary layers where gradients of deformation may arise.

The standard (Cauchy or first gradient) continuum models were conceived under some implicit assumptions which do not allow for the description of all just listed extreme mechanical properties [START_REF] Dell'isola | At the origins and in the vanguard of peridynamics, nonlocal and higher-gradient continuum mechanics: An underestimated and still topical contribution of Gabrio Piola[END_REF].

For this reason, it is needed to introduce, for pantographic sheets, a generalized continuum model, by reconsidering, at the very beginning, the standard modeling procedure. Indeed a consistent model should allow for, in particular, the description of deformation gradient concentrations induced, for instance, by imposed boundary displacements. To be driven in the construction of the most suitable model we use the asymptotic homogenization method extensively presented in ( [START_REF] Boutin | Rayleigh scattering in elastic composite materials[END_REF]). It allows for the rigorous construction of a description able to encompass all mentioned atypical properties.

Even if we limit ourselves to the case of linearized models in statics (and in absence of body forces) we still get a nonstandard second gradient continuum model. Moreover, the advantage of used micro-macro upscaling procedure is that it allows us to determine firstly the atypical structure of the constitutive and balance equations and secondly all relative constitutive parameters (i.e. those specifying first gradient and second gradient terms in deformation energy). This is done explicitly in terms of the microscopic properties of the elementary cells constituting the pantographic sheet. Furthermore, the presented micro-macro identification provides a design rule for pantographic bidimensional continua.

The main advances provided by this paper compared to some previous works on pantographic sheets lie in the three following points:

• The homogenization method leads to a clear and rigorous micro-macro identification procedure.

Hence, the structure of the homogenized description and its parameters are explicitly deduced form the cell. This differ from the "a priori" variational approaches where a macroscopic description is proposed but its applicability to specific micro-structure can only be postulated. Beside, compared to direct numerical simulations at the microscale [START_REF] Dell'isola | Pantographic 2D sheets: Discussion of some numerical investigations and potential applications[END_REF] the established model presents a great advantage, in terms of computational cost but also of in terms of physical understanding of the actual behaviour.

• The description enables to account for the low extensibility and the high bending deformability. It thus enable to overcome the drawback of the studies that consider inextensible fibers [START_REF] Dell'isola | Plane bias extension test for a continuum with two inextensible families of fibers: a variational treatment with Lagrange Multipliers and a perturbation solution[END_REF].

• The physical insight of pantographic systems is improved by a simple interpretation of the macroscopic description: the tension in a fiber varies because of the shear forces transferred through the pivots by the orthogonal fibers. Even if the paper focuses on small deformations, this mechanism remains the essentially the same (meanwhile complicated by the non orthogonality of the fibers) at large deformation [START_REF] Dell'isola | Large deformations of planar extensible beams and pantographic lattices: Heuristic homogenisation, experimental and numerical examples of equilibrium[END_REF].

Note that we only consider small in plane deformations. However the same homogenization method can be extended to investigate buckling and out plane deformation as discussed in [START_REF] Giorgio | Buckling modes in pantographic lattices[END_REF][START_REF] Giorgio | Pattern formation in the threedimensional deformations of fibered sheets[END_REF].

The obtained models can be framed in the context of generalized continuum theories and immediately allows for the explicit determination of deformation energy and related boundary conditions. The atypical features of described continuum model for pantographic sheets are reflected by its atypical mathematical properties. That is:

• the PDEs determining the equilibrium configurations involve second and forth order partial derivatives terms which may have a comparable order of magnitude; in the studied case of orthogonal fibers, the displacement fields along the fibers result to be governed by uncoupled PDEs;

• the deformation energy is not coercive in the standard sense. However, using specific boundary conditions, consistently established from the found expression for deformation energy, we establish the definite positiveness of the deformation energy functional, which simultaneously involves first and second order of displacement gradients; we conjecture that such formulated problems are well-posed.

The corresponding atypical physical properties of pantographic sheets consist in the following circumstances:

• the extensional forces along a fiber are balanced by the shear forces due to the bending of its transverse fibers (see equations ( 44), ( 45))

• the balance of couples is separately valid for the two arrays of orthogonal fibers, so that the fiber rotation is proportional to the axial gradient of transverse displacement, see formula [START_REF] Eremeyev | On effective properties of materials at the nano-and microscales considering surface effects[END_REF].

The numerical simulations which we have presented are aimed to illustrate the performances of obtained model. Indeed pantographic sheets show very peculiar deformation patterns exhibiting regions of concentrated deformation energy. In particular when applied to rectangular specimen the model:

• allows for the determination, already in the linear case, of the region where the deformation energy is localized without any further a priori assumptions;

• shows that the diffusion patterns of deformation inside the specimen differ notably from those shown in first gradient material. In particular, they are canalized in a way determined by the material symmetry and boundary conditions, while their thickness is determined by the characteristic length specified by the competition between the first and second gradient stiffnesses;

• the described features are confirmed by considering specimen with circular shape also when including rectangular holes.

All these features and predictions match at least qualitatively the experimental evidence e.g. [START_REF] Dell'isola | Designing a light fabric metamaterial being highly macroscopically tough under directional extension: first experimental evidence[END_REF]. It then appears that pantographic sheet is an archetypical oriented material in which second gradient effect plays a crucial role and for which it is possible to get a close description of the underlying actual physical mechanism.

There are many worthy issues that may be investigated concerning the complex structure here analyzed:

• linearized analysis of pantographic structures with non-orthogonal and uneven fibers (this case is also relevant as it may give a tangent model for large deformations) and the determination of some analytical solutions of found PDEs;

• the comparison of presented theory with further experimental evidence for possible improvement in the idealized model;

• the exploration of the application of the model to the mechanics of fabrics, the analysis of the dissipation that can arise from the interaction of the fibers or from internal friction (see e.g. [START_REF] Nadler | A model for frictional slip in woven fabrics[END_REF]) and the extension to three-dimensional deformations;

• buckling phenomena which can produce wrinkling (see e.g. [START_REF] Giorgio | Buckling modes in pantographic lattices[END_REF][START_REF] Rizzi | The effects of warping on the postbuckling behaviour of thin-walled structures[END_REF][START_REF] Rizzi | Initial postbuckling behavior of thin-walled frames under mode interaction[END_REF][START_REF] Gabriele | On the imperfection sensitivity of thin-walled frames[END_REF][START_REF] Carassale | Non-linear discrete models for the stochastic analysis of cables in turbulent wind[END_REF][START_REF] Piccardo | Some research perspectives in galloping phenomena: critical conditions and post-critical behavior[END_REF]) and damage detection and its evolution (see e.g. [START_REF] Yang | Higher-order continuum theory applied to fracture simulation of nanoscale intergranular glassy film[END_REF][START_REF] Misra | Micromechanical model for viscoelastic materials undergoing damage[END_REF][START_REF] Andreaus | Identification of multiple open and fatigue cracks in beam-like structures using wavelets on deflection signals[END_REF][START_REF] Placidi | A variational approach for a nonlinear 1-dimensional second gradient continuum damage model[END_REF][START_REF] Placidi | A variational approach for a nonlinear one-dimensional damage-elasto-plastic secondgradient continuum model[END_REF]).

Finally, as future and more challenging lines of development, considering biological applications, it is conceivable a morphoelastic pantographic continuum which can be employed as a 'smart' prosthesis (see e.g. [START_REF] Mcmahon | Nonlinear morphoelastic plates II: exodus to buckled states[END_REF]). In this context, it is possible also to equip the pantographic structure with transducers that can act both as sensors or actuators in order to induce some functional adaptation capabilities (see e.g. [START_REF] D'annibale | Linear stability of piezoelectric-controlled discrete mechanical systems under nonconservative positional forces[END_REF][START_REF] D'annibale | On the failure of the 'similar piezoelectric control' in preventing loss of stability by nonconservative positional forces[END_REF][START_REF] Andreaus | Effect of electrical uncertainties on resonant piezoelectric shunting[END_REF]).

Figure 1 :

 1 Figure 1: Topology of the pantographic lattice. The pivots are equally distributed on the whole sheet. The black points are indicating specific monitoring points during bias test experimentations.

3

  rigidities, respectively. The beams slenderness hypothesis √ A implies that such rigidities differ significantly. In fact, considering beams of regular section, say for instance rectangular section of length's sides a and b with b = O(a) = O( √ A), then A = ab, I = ba 3 /12 and 12I/A = a 2 = O(A), and consequently the rigidity ratio R is

  Thus by limiting ourself to the infinitesimals O(ε 6 ), the continuity (14) of moments of the elements oriented along y gives, after grouping different terms,

Figure 2 :

 2 Figure 2: An example of standard bias extension test. Courtesy by Tomasz Lekszycki, Marek Pawlikowski and Roman Grygoruk.

Figure 3 :

 3 Figure 3: An example of generalized bias test with imposing a shear displacement.

Figure 4 :

 4 Figure 4: Sketch of vectors n and N on the boundaries.

Figure 5 :

 5 Figure 5: Equilibrium shapes of rectangular specimens submitted to extension (left, a) shear (center, b) and bending (right, c) displacement. The colours indicate the levels of stored energy density.

Figure 6 :

 6 Figure 6: Equilibrium shapes of circular specimens submitted to extension (left, a) shear (center, b) and bending (right, c) displacement. The colours indicate the levels of stored energy density.

Figure 7 :

 7 Figure 7: Equilibrium shapes for extension bias test of a circular sample with a squared hole: a) whose edges are oriented along fibers; b) rotated of an angle of π/8; c) rotated of an angle of π/4.The colors indicate the levels of stored energy density.

Figure 8 :

 8 Figure 8: Equilibrium shapes for bending test of a circular sample with a squared hole: a) whose edges are oriented along fibers; b) rotated of an angle of π/8; c) rotated of an angle of π/4. The colours indicate the levels of stored energy density.

  x n x n y n y ; dU x /dn dual of T yy x n 2 y U y dual of T y y n y + -T xx y,x (1 + n 2 y ) + T xx y,y n x n y n x ; dU y /dn dual of T xx y n 2

	x

Acknowledgements

Claude Boutin gratefully acknowledges the Università di Roma La Sapienza for a visiting grant. The authors thank Tomasz Lekszycki, Marek Pawlikowski and Roman Grygoruk for having made available the fig. 2.