Houssam Abdoul-Anziz

Cédric Bellis

Pierre Seppecher
email: seppecher@imath.fr

Keywords: Periodic homogenization, second gradient, strain-gradient, Gamma-convergence Mathematics Subject Classification: 35B27, 35J30

We determine in the framework of static linear elasticity the homogenized behavior of three-dimensional periodic structures made of welded elastic bars. It has been shown that such structures can be modeled as discrete systems of nodes linked by extensional, flexural/torsional interactions corresponding to frame lattices and that the corresponding homogenized models can be strain-gradient models, i.e., models whose effective elastic energy involves components of the first and the second gradients of the displacement field. However, in the existing models, there is no coupling between the classical strain and the strain-gradient terms in the expression of the effective energy. In the present article, under some assumptions on the positions of the nodes of the unit cell, we show that classical strain and strain-gradient strain terms can be coupled. In order to illustrate this coupling we compute the homogenized energy of a particular structure which we call asymmetrical pantographic structure.

Homogenization of frame lattices leading to second gradient models coupling classical strain and strain-gradient terms Introduction On the one hand, macroscopic behavior of elastic composite materials has been widely studied both in the mechanics as well as in the mathematics literature. From the mathematical point of view, the homogenization of periodic elastic media with moderate contrast is well founded, see, for instance 10,35 . It consists in taking into account the fact that the size of the unit cell is much smaller than the characteristic size L of the whole sample of the material and in passing to the limit when the ratio ε := /L tends to zero. This approach called asymptotic homogenization has been widely used to study conduction or elasticity problems in static or dynamic cases 17,28,30,35 . The effective behavior of the periodic elastic medium is characterized explicitly in terms of a local minimization problem set in a rescaled cell (see formula (3.6) of 6 ).

On the other hand, strain-gradient (or second gradient) materials are expected to have an exotic behavior 24,36 . They are obtained in the literature by using different methods: heuristic homogenization techniques 12,13,25 and variational methods 4,9,15,31 . In articles 4,9,15,31 the second gradient energies obtained involve only the components of the skew-symmetric part of the gradient of the displacement: the limit models obtained enter the framework of the so-called couple-stress models.

computed from the components of the second gradient ∇∇u of the displacement field. As the inverse is also true, strain-gradient materials can be equivalently called second gradient materials.

In classical linear elasticity, the elastic energy density is a non-negative quadratic form * Q(e(u)) = 1 2 e(u) : A : e(u)

of the strain where A is the rigidity tensor (symmetric fourth-order tensor: A ijkl = A klij = A jikl ) and the stress tensor is its differential σ = A : e(u). In linear strain-gradient elasticity the elastic energy density is a non negative quadratic form of the strain and its gradient Q(e(u), ∇e(u)) = 1 2 e(u) : (A : e(u)) + e(u) : (B :

. ∇e(u)) + 1 2 ∇e(u) :

. (C :

. ∇e(u)) where C is a symmetric sixth-order tensor (C ijklmn = C lmnijk = C jiklmn ) and B is a fifth-order tensor (with the symmetry B ijklm = B jiklm = B ijlkm ).

Equilibrium equations under the action of some external body force f are easily recovered from the minimization of the total energy

W (u) := Ω Q(e(u), ∇e(u)) -f • u dx.
This is a lower semi-continuous (for instance for the L 2 topology) functional over the displacement field u and under suitable coercivity assumptions, the existence of a minimizer is ensured. The variational formulation of the minimization problem reads ∀v, Ω (e(u) : A : e(v) + e(u) : B :

. ∇e(v) + e(v) : B :

. ∇e(u) + ∇e(u) : . C :

. ∇e(v)f • v) dx = 0.

Let us define the third and second-order tensors H := B T : e(u) + C :

. ∇e(u) and σ := A : e(u) + B :

. ∇e(u)div(H),

where B T stands for the transposed tensor: (B T ) ijklm := (B) klmij . Integrating by parts, the Euler equation of the minimization problem takes, on the interior of the domain, the usual form div(σ) + f = 0 (3) of equilibrium equation. The fundamental difference with classical elasticity is that, now, the stress tensor σ depends on the second and third-order partial derivatives of u. The boundary conditions are also fundamentally different from these in classical elasticity: the boundary terms which appear when integrating by parts give the natural free boundary conditions associated with the model. In case of a smooth domain, they read

(H • n) • n = 0 and σ • n + div // (H • n) // = 0. ( 4
)
where n stands for the exterior normal to the domain and div // for the surface divergence operator on the boundary † . We refer to 20-23 for the interpretation of the first equation in terms of double-forces, the dependence of surface force on the curvature of the surface (through div // n) as well as the possible presence of edge forces in case of a non smooth boundary (not presented here). In these articles the case when non vanishing forces or double-forces are applied on the boundary are detailed and the dual conditions (in terms of imposed displacement and normal derivative of the displacement) are also considered.

To understand the effect of the different terms and of the boundary conditions on the equilibrium, let us consider, for the time being, a one-dimensional case: a bar made of a strain-gradient material is fixed at its center while a non vanishing displacement u 0 is imposed at its left hand extremity (see Figure 1). Thus the left hand half is expanding while the right half of the bar is left free. In ordinary (first gradient) elasticity the solution is obvious: the left half

(a ⊗ b) ij := a i b j , (A • B) i 1 ...in := A i 1 ...ip j B j i p+1 ...in (A : B) i 1 ...in := A i 1 ...ip j k B k j i p+1 ...in , (A : . B) i 1 ...in := A i 1 .
..,ip j k l B l k j i p+1 ...in . † For any second-order tensor M , notation M // stands for the projection

M // := M -(M • n) ⊗ n
of the bar is subjected to a uniform dilatation while the right half remains at rest (see the case = 0 in Figure 2a). For a strain-gradient material the situation is more complex: the dilation tends to propagate in the free right part of the bar and when a coupling term is present this propagation has its own behavior. Let us perform completely the minimization of the total energy for this one-dimensional example: the domain Ω = (-L, L) is an interval and the total energy reads

L -L Q(u , u ) dx = L -L 1 2 a(u (x)) 2 + bu (x)u (x) + 1 2 c(u (x)) 2 dx
under the constraints u(0) = 0, u(-L) = u 0 . Introducing the intrinsic length := c √ ac-b 2 and adimensional parameter k := b c , the energy is proportional to

1 2 L -L (u (x)) 2 + ( u (x) + ku (x)) 2 dx.
We can see that two terms are in competition: the classical term, proportional to (u ) 2 , tends to minimize the strain while the other one tends, when k = 0 to minimize the variation of the strain, that is to favor a constant strain and when k = 0 to favor an exponentially varying strain. The equilibrium is the solution on (-L, 0) and (0, L) of the linear ordinary differential equation

2 u -(1 + k 2 )u = 0
with the boundary and continuity conditions u

(-L) = u 0 , u (-L) + ku (-L) = 0, 2 u (L) -(1 + k 2 )u (L) = 0, u (L) + ku (L) = 0, u(0 -) = u(0 + ) = 0, u (0 + ) = u (0 -) and u (0 -) + ku (0 -) = u (0 + ) + ku (0 +
). The influence of the parameters and k over the equilibrium are shown in Figure 2. The effect of length is clear.

The expansion imposed to the left half decreases exponentially in the free right half. This length is characteristic of that decreasing. The effect of the adimensional coupling parameter k is more intriguing.

Pantographic structures

Even if the results of 3 give a straightforward procedure for analyzing the strain-gradient effective behavior of periodic frame lattices, it is not obvious to exhibit structures which actually show such a behavior. Indeed one has to conceive structures which are able to propagate the strain and, up to now, very few structures are known to do so. Pantographic structures are among them 4,5,19,27,32-34,36-38 : they lead through homogenization process to interesting strain-gradient models. A way for building them is to consider a periodic array of cells made of six nodes, linked by elastic beams as represented in Figure 3. The "pantographic" beam thus obtained can be used to build 2D or 3D materials by arranging many parallel pantographic beams and linking them in a suitable way (see Figures 45).

Let us be more specific: the considered periodic lattice is made of nodes y (i,j,k),s := ε(y s + it 1 + jt 2 + kt 3 ) where y s stands for the position of node s (s ∈ {1, . . . , 6}) in the prototype rescaled cell, t 1 , t 2 , t 3 are three independent periodicity vectors and i, j and k vary from 1 to ε -1 (which is assumed to be an integer). Fixing

y 1 = 1 4 , 0, 0 , y 2 = 0, - 1 4 , 0 , y 3 = 0, 1 4 , 0 , y 4 = 1 2 , 1 2 , 0 , y 5 = 1 2 , - 1 2 , 0 , y 6 = 3 4 , 0, 0 , (5)
assuming that, inside each cell, only pairs of nodes [START_REF] Abdoul-Anziz | Homogénéisation de composites élastiques périodiques à fort contraste: Conception de métamatériaux de second gradient[END_REF][START_REF] Abdoul-Anziz | Homogenization of periodic graph-based elastic structures[END_REF], [START_REF] Abdoul-Anziz | Homogénéisation de composites élastiques périodiques à fort contraste: Conception de métamatériaux de second gradient[END_REF][START_REF] Abdoul-Anziz | Strain gradient and generalized continua obtained by homogenizing frame lattices[END_REF], [START_REF] Abdoul-Anziz | Homogénéisation de composites élastiques périodiques à fort contraste: Conception de métamatériaux de second gradient[END_REF][START_REF] Alibert | Second-gradient continua as homogenized limit of pantographic microstructured plates: a rigorous proof[END_REF], [START_REF] Abdoul-Anziz | Homogénéisation de composites élastiques périodiques à fort contraste: Conception de métamatériaux de second gradient[END_REF][START_REF] Alibert | Truss modular beams with deformation energy depending on higher displacement gradients[END_REF], [START_REF] Abdoul-Anziz | Homogenization of periodic graph-based elastic structures[END_REF][START_REF] Alibert | Truss modular beams with deformation energy depending on higher displacement gradients[END_REF], [START_REF] Abdoul-Anziz | Strain gradient and generalized continua obtained by homogenizing frame lattices[END_REF][START_REF] Alibert | Second-gradient continua as homogenized limit of pantographic microstructured plates: a rigorous proof[END_REF], [START_REF] Alibert | Second-gradient continua as homogenized limit of pantographic microstructured plates: a rigorous proof[END_REF][START_REF] Allaire | Homogenization and two-scale convergence[END_REF], [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] are interacting, assuming that nodes 4 and 5 of each cell are respectively interacting with nodes 3 and 2 of the next cell following εt 1 (with t 1 = (1, 0, 0)) and assuming moreover that node 6 of each cell is interacting with both nodes 3 and 2 of the next cell following εt 1 , we get a pantographic beam as shown in Figure 6. Assuming moreover that node 1 of each cell is interacting with node 1 of the next cells following εt 2 and εt 3 (with t 2 = (0, √ 3, -1) and t 3 = (0, √ 3, 1)) we get a 3D structure as shown in Figures 45. The procedure described in 3 gives the effective behavior of this structure which can be considered as a microstructured material. The homogenized elastic energy is a quadratic functional of the displacement field u. It reads

E(u) = 1 2 Ω (e 13 (u)) 2 + σ (e 12 (u)) 2 + γ (e 11 (u)) 2 + κ( ∂ 2 u 1 ∂x 2 1 2 + ∂ 2 u 2 ∂x 2 1 2 dx
under the constraints e 22 (u) = e 33 (u) = e 23 (u) = 0. The values of material parameters , σ, γ, κ are explicit in terms of the rigidities of the bars the structure is made of. Note also that the second partial derivatives could be written in terms of the strain-gradient. As noticed in 3 this energy does not couple strain and strain-gradient. For such a coupling, a term like

∂ 2 u 1 ∂x 2 1 - ∂u 1 ∂x 1 2
would be necessary. As, at equilibrium, the continuum tends to minimizes its energy, the displacement field would tend to fulfill the equation

∂ 2 u 1 ∂x 2 1 - ∂u 1 ∂x 1 = 0
that is would tend to increase exponentially in the direction t 1 . How could that be possible ? The answer is natural: modifying the pantograph we just described, making it asymmetrical, must lead to a structure which not only propagates the dilatation but also modifies it from cell to cell. However the attempt to simply modify the coordinates of the nodes given in (5) at order O(1) in the rescaled cell is bound to fail: indeed one would remain in the framework studied in 3 where it is proved that the homogenized energy density is the sum of a quadratic form of the strain-gradient due to the extensional rigidities of the bars and of a quadratic form of the classical strain due to the torsion or bending rigidities of the bars. Hence no coupling can arise. The remark on which we base our study is that, even though it seems a natural idea, there is no reason for assuming that the positions y s of the nodes inside the rescaled cell do not depend on ε. We thus propose to modify the positions y 4 and y 5 of nodes 4 and 5 as

y ε 4 := 1 2 + α ε, 1 2 , 0 and y ε 5 := 1 2 + α ε, - 1 2 , 0 ,
respectively, where α is a parameter which allows us to tune the asymmetry of our structure.

Prepared using sagej.cls

The new pantographic structure is now very slightly asymmetrical as shown in Figure 6 but there is no need to draw the new 3D structure which remains essentially similar to the one represented in Figures 45. Next section is devoted to the study of the effective properties of periodic lattices when the position of the nodes in the rescaled cell is allowed still to depend on ε. Later on, the homogenization result will be applied to the asymmetrical pantographic structure we just described and we will indeed verify that the effective model is a strain-gradient model with coupling.

① ③ ② ⑥ ④ ⑤ εt 1

Statement of the homogenization problem

Considered geometry

We restrict ourselves to 3D structures leading to 3D continuum models. The study of structures leading to 2D or 1D models, that is to membranes, plates and beams can be treated in a very similar way. This has been done in 3 for non-coupling cases.

So we consider in R 3 a periodic structure with three independent vectors of periodicity t 1 , t 2 and t 3 . The domain Ω is the parallelepiped Ω := {x 1 t 1 + x 2 t 2 + x 3 t 3 : (x 1 , x 2 , x 3 ) ∈ [0, 1] 3 }. We consider a discrete lattice made of a finite number of points (called "nodes") repeated periodically following the vectors t i . We use I := (i 1 , i 2 , i 3 ) ∈ Z 3 to label the different cells of the lattice and s ∈ {1, . . . , K} to label the different nodes in each cell. The geometry of the lattice depends on a small parameter ‡ ε and the position of node s of cell I is assumed to be of the form

y ε I,s := ε(y s + εz s + i 1 t 1 + i 2 t 2 + i 3 t 3 ), (6)
where (y s ) and (z s ) are two fixed families of K vectors of R 3 . The considered lattice is made of those nodes which lie inside Ω and is therefore divided in ε -3 cells similar to εΩ.

Note that the role of the parameter ε is twofold: first it measures the ratio of the sizes of the periodic cell and of the considered domain (in that sense, the fact that ε is very small is the standard hypothesis on which homogenization is based); secondly one can see in [START_REF] Allaire | Homogenization and two-scale convergence[END_REF] that it measures the speed of convergence of the nodes positions y s + εz s in the rescaled prototype cell towards their limit positions y s . Here lies the main difference with previous works 2,3 where these rescaled positions were assumed not to depend on ε.

We introduce

y ε I := ε(i 1 t 1 + i 2 t 2 + i 3 t 3 )
as a reference point in the cell I and

I ε := (i 1 , i 2 , i 3 ) ∈ 0, ε -1 3
the set of indices of the cells of the structure. In the sequel, we will use the notation I φ I,s to denote the mean value of any field φ I,s defined at the nodes of the lattice:

I φ I,s := ε 3 I∈I ε φ I,s .
‡ For sake of simplicity, ε -1 is assumed to be an integer.

In order to describe the interactions between the nodes of a generic cell I and the interactions between the nodes of this cell with the nodes of its 26 closest neighbor cells I + p (actually, due to periodicity, only half of them have to be considered), we introduce the set P := {(0, 0, 0), (1, 0, 0), (0, 1, 0), (1, 1, 0), (1, -1, 0), (0, 0, 1), (0, 1, 1), (1, 0, 1), (0, 1, -1), (-1, 0, 1), (1, 1, 1), (1, 1, -1), (-1, 1, 1), (1, -1, 1)} that we identify with {1, . . . , 14}. For any p := (p 1 , p 2 , p 3 ) ∈ P we denote p := p 1 t 1 + p 2 t 2 + p 3 t 3 the corresponding vector so that y ε I+p,s = y ε I,s + εp. For any pair of distinct nodes (y ε I,s , y ε I+p,s ), we introduce the associated rescaled length and direction:

ε p,s,s := ε -1 y ε I+p,s -y ε I,s , τ ε p,s,s := y ε I+p,s -y ε I,s ε ε p,s,s
.

Using ( 6) we get the expression

τ ε p,s,s = y s -y s + p + ε(z s -z s ) y s -y s + p + ε(z s -z s ) .
A Taylor series expansion in ε gives

τ ε p,s,s = τ p,s,s + ερ p,s,s + ε 2 ζ ε p,s,s , (7)
with ζ ε p,s,s ≤ C (C being a constant independent of ε), τ p,s,s := (y s -y s + p)/ p,s,s , ρ p,s,s := 1 p,s,s (Id -τ p,s,s ⊗ τ p,s,s ) • (z s -z s ), (8)
and p,s,s := y sy s + p , where Id stands for the identity tensor.

Mechanical interactions

The nodes of the structure are linked by elastic bars made of a homogeneous isotropic elastic material. These bars are slender cylinders with circular cross-section whose radius is of the form r ε := βε 2 , so that flexural and torsional energies are of order ε 2 compared to the extensional energy. We will not do here the asymptotic study of slender cylinders for which extension, bending and torsional rigidities can be found in many mechanics textbooks, for instance in 26 . We will neither do the integration along the bars (see 2 ) which reduces to a quadratic energy involving rigid motions associated to the two interacting nodes. Instead we will start directly from this quadratic interaction energy. Any node of the structure is endowed with a rigid motion, that is a couple of vectors (U I,s , θ I,s ) (remember that in small deformations, rotations are represented by vectors). Hence the kinematics of the structure is represented by the two families of vectors (i.e. discrete fields) U := (U I,s ) and θ := (θ I,s ).

For any interacting pair of distinct nodes (I, s), (I + p, s ) , we introduce the extension and the global rotation of the bar by setting respectively

(ρ u ) I,p,s,s := U I+p,s -U I,s ε • τ ε p,s,s and (α u ) I,p,s,s := U I+p,s -U I,s ε p,s,s × τ ε p,s,s .
Since the bar linking the nodes (I, s) and (I + p, s ) is elastic, an energy proportional to (ρ u ) 2 I,p,s,s is needed for extending it while the energy needed to bend and twist it is a quadratic form of (α u ) I,p,s,sθ I,s and (α u ) I,p,s,sθ I+p,s .

Let us denote A ⊂ P × {1, . . . , K} 2 the set of indices (p, s, s ) such that nodes (I, s) and (I + p, s ) are linked by a bar. Without loss of generality, we assume that a cell is in interaction only with its closest neighbors. Indeed, the periodic cell can always been chosen large enough for this assumption to be satisfied as illustrated in Figure 7.

We are not interested in disconnected structures. We make, like in 2,3 , the "connectedness assumption" that the bars connect all the nodes of the lattice. More precisely § we assume that, for any p ∈ P and any (s, s ) ∈ {1, . . . , K} 2 , there exist a finite path in the structure which connects the node (I, s) to the node (I + p, s ) that is a finite sequence (s 1 , . . . , s r+1 ) in {1, . . . , K}, (p 1 , . . . , p r ) in P,

• • • • • y 1 • • • • • • • • y 2 • • • • • • • • • • • • • • • • • • • (a) Long-range interactions: 16 cells • y 1 • y 2 • y 3 • y 4 • • • • • y 8 • y 7 • y 6 • y 5 • • • • • • • • • • • • • • • • • • • • (b)
( 1 , . . . , r ) in {-1, 1} such that s 1 = s, s r+1 = s , r i=1 i p i = p, i > 0 =⇒ (p i , s i , s i+1 ) ∈ A and i < 0 =⇒ (p i , s i+1 , s i ) ∈ A.
The global extensional energy reads

E ε (U ) := ε -2 I (p,s,s )∈A a p,s,s 2 (ρ u ) 2 I,p,s,s . (9)
Here

a p,s,s = Yπβε 3 p,s,s (10)
where Y stands for the Young modulus of the material the bars are made of. We assume that the Young modulus is of order ε -3 . This simply means that we have chosen the energy unit in such a way that the positive coefficients a p,s,s are of order one.

The global bending/torsional energy reads

F ε (U , θ) := I (p,s,s )∈A (θ I,s -(α u ) I,p,s,s ) • B p,s,s 2 • (θ I,s -(α u ) I,p,s,s ) + (θ I,s -(α u ) I,p,s,s ) • C p,s,s • (θ I+p,s -(α u ) I,p,s,s ) + (θ I+p,s -(α u ) I,p,s,s ) • D p,s,s 2 • (θ I+p,s -(α u ) I,p,s,s ) . (11)
Here

B p,s,s = D p,s,s = β 2 a p,s,s Id - 3 + 4ν 4 + 4ν τ p,s,s ⊗ τ p,s,s , C p,s,s = β 2 a p,s,s 2 Id - 3 + 2ν 2 + 2ν τ p,s,s ⊗ τ p,s,s ,
where ν is the Poisson ratio of the material the bars are made of. For any (p, s, s ) ∈ A, the block matrix

  B p,s,s C p,s,s C T p,s,s D p,s,s   is positive.
To summarize, the total elastic energy of the considered structure reads

E ε (U , θ) := E ε (U ) + F ε (U , θ). (12)
We associate to any sequence U ε the families of vectors m ε I , v ε I,s and χ ε I,p defined by

m ε I := 1 K K s=1 U ε I,s , v ε I,s := ε -1 (U ε I,s -m ε I ), χ ε I,p := ε -1 (m ε I+p -m ε I ) (13)
and the family of real numbers ω ε I,p,s,s defined by

ω ε I,p,s,s := ε -2 (U ε I+p,s -U ε I,s ) • τ ε p,s,s if (p, s, s ) ∈ A, 0 otherwise. ( 14
)
Using this notation, we can rewrite the two addends E ε (U ε ) and F ε (U ε , θ ε ) of the total elastic energy E ε (U ε , θ ε ) under the following forms:

E ε (U ε ) ≡ Ēε (v ε , χ ε ) = ε -2 I (p,s,s )∈A a p,s,s 2 (v ε I+p,s -v ε I,s + χ ε I,p ) • τ ε p,s,s 2 , (15)
F ε (U ε , θ ε ) ≡ Fε (v ε , χ ε , θ ε ) = I (p,s,s )∈A θ ε I,s -(α ε v,χ ) I,p,s,s • B p,s,s 2 • θ ε I,s -(α ε v,χ ) I,p,s,s + θ ε I,s -(α ε v,χ ) I,p,s,s • C p,s,s • θ ε I+p,s -(α ε v,χ ) I,p,s,s (16)
+ θ ε I+p,s -(α ε v,χ ) I,p,s,s • D p,s,s 2 θ ε I+p,s -(α ε v,χ ) I,p,s,s ,
where (α ε v,χ ) I,p,s,s :=

τ ε p,s,s p,s,s × (v ε I+p,s -v ε I,s + χ ε I,p ).

Asymptotic expansion

In order to guess the limit of E ε as ε tends to zero, let us first compute the limit of E ε (U ε , θ ε ) for particular sequences of displacement fields of the form

U ε I,s := u(y ε I ) + εv s (y ε I ) + ε 2 w s (y ε I ) and θ ε I,s := θ s (y ε I ) (17)
where u, v s , w s , θ s are smooth functions satisfying

K s=1 v s = K s=1 w s = 0.
Remind that y ε I is the reference point in cell I. We still associate to (U ε , θ ε ) the quantities m ε I , v ε I,s , χ ε I,p and ω ε I,p,s,s as in ( 13) and ( 14). We have

m ε I = u(y ε I ), v ε I,s = v s (y ε I ) + εw s (y ε I ), χ ε I,p = ε -1 (u(y ε I+p ) -u(y ε I )) = ∇u(y ε I ) • p + ε 2 ∇∇u(y ε I ) : (p ⊗ p) + O(ε 2 ) and v ε I+p,s -v ε I,s = v s (y ε I+p ) -v s (y ε I ) + ε(w s (y ε I+p ) -w s (y ε I )) = v s (y ε I+p ) -v s (y ε I ) + v s (y ε I ) -v s (y ε I ) + ε w s (y ε I+p ) -w s (y ε I ) = ε∇v s (y ε I ) • p + v s (y ε I ) -v s (y ε I ) + ε w s (y ε I+p ) -w s (y ε I ) + O(ε 2 ). Hence (v ε I+p,s -v ε I,s + χ ε I,p ) • τ ε p,s,s = (∇u(y ε I ) • p + v s (y ε I ) -v s (y ε I )) • τ p,s,s + O(ε)
and the extensional energy E ε (U ε ) cannot remain finite unless the functions u and v s satisfy the constraint

(∇u • p + v s -v s ) • τ p,s,s = 0. ( 18
)
So now we restrict ourselves to functions satisfying [START_REF] Maso | An Introduction to Γ-Convergence[END_REF]. In that case we get

(v ε I+p,s -v ε I,s + χ ε I,p ) • τ ε p,s,s = ε 2 ∇∇u(y ε I ) : (p ⊗ p) + ε∇v s (y ε I ) • p + ε(w s (y ε I ) -w s (y ε I )) • τ p,s,s + ε (∇u(y ε I ) • p + v s (y ε I ) -v s (y ε I )) • ρ p,s,s + O(ε 2
), with ρ p,s,s defined in [START_REF] Barbarosie | Optimization of bodies with locally periodic microstructure by varying the periodicity pattern[END_REF]. By Riemann sum,

lim ε→0 Ēε (v ε , χ ε ) = lim ε→0 ε -2 I (p,s,s )∈A a p,s,s 2 (v ε I+p,s -v ε I,s + χ ε I,p ) • τ ε p,s,s 2 = Ω (p,s,s )∈A a p,s,s 2
1 2 ∇∇u(x) : (p ⊗ p) + ∇v s (x) • p + w s (x) -w s (x) • τ p,s,s + (∇u(x) • p + v s (x) -v s (x)) • ρ p,s,s 2 dx.
Prepared using sagej.cls

In order to write this result in a more synthetic way, let us fix some notation. For any functions u, v s and w s in L 2 (R 3 ), we set in the sense of distributions, for any ¶ (p, s) ∈ P × {1, . . . , K},

(η u ) p,s := ∇u • p, (19)
(ξ u,v ) p,s := ∇v s • p + 1 2 ∇∇u : (p ⊗ p) (20)
and we define

Ē(w, ξ, v, η) := Ω (p,s,s )∈A a p,s,s 2 ((w s -w s + ξ p,s ) • τ p,s,s + (v s -v s + η p,s ) • ρ p,s,s ) 2 dx. (21)
Using this notation, the previous limit reads lim ε→0 Ēε (v ε , χ ε ) = Ē(w, ξ u,v , v, η u ) and constraint ( 18) can be written Ē(v, η u , 0, 0) = 0.

Computing the bending-torsion energy is simpler: we have

(α ε v,χ ) p,s,s = τ ε p,s,s p,s,s × (v ε I+p,s -v ε I,s + χ ε I,p ) = τ p,s,s p,s,s × (v s (y ε I ) -v s (y ε I ) + ∇u(y ε I ) • p) + O(ε).
Therefore, setting (α v,ηu ) p,s,s := τ p,s,s p,s,s

× (v s -v s + (η u ) p,s ) and F (v, η, θ) = Ω (p,s,s )∈A (θ s -(α v,η ) p,s,s ) • B p,s,s 2 • (θ s -(α v,η ) p,s,s ) + (θ s -(α v,η ) p,s,s ) • C p,s,s • (θ s -(α v,η ) p,s,s ) + (θ s -(α v,η ) p,s,s ) • D p,s,s 2 • (θ s -(α v,η ) p,s,s ) dx (22)
we get lim ε→0 F ε (U ε , θ ε ) = F (v, η u , θ).

To conclude, the limit of the elastic energy E ε (U ε , θ ε ) is Ē(w, ξ u,v , v, η u ) + F (v, η u , θ) under the constraint Ē(v, η u , 0, 0) = 0. At this point the reader may consider that the effective behavior of the structure has been ascertained. That is not quite true as sequences of equilibrium solutions have no reason to fulfill the assumption [START_REF] Cioranescu | An Introduction to homogenization[END_REF]. It remains to prove that the result remains true in the general case.

Homogenization framework

To state a homogenization result one has to fix the relationship between sequences of equilibrium displacement fields and their macroscopic limits. In the standard case this relationship is natural: all displacement fields belong to the same functional space (for instance L 2 (Ω)) and their limit is taken is the sense of a suitable topology on this space. We cannot do that here, as, for different values of ε, the displacement field is defined on different sets of nodes. Moreover at fixed ε the displacement field is a finite family of vectors while we expect the effective displacement field to be a continuous function on Ω.

A usual way for transferring the problem on a fixed functional space is to introduce an extension over Ω of the fields which are defined only at the nodes of the structure. As such an extension is somehow arbitrary, we prefer to see the displacement field at fixed ε as a discrete measure and to define its continuous limit as the diffuse limit measure.

Definition 1. Let Z ε = (Z ε I ) I∈I ε be a sequence of families of real numbers. We say that Z ε converges to a measurable function z, and we write Z ε z, when the following weak- * convergence of measures holds:

I Z ε I δ y ε I * z(x) dx ( 23
)
where δ y ε I stands for the Dirac measure at the point y ε I . This definition is extended to vector fields by applying it component-wise.

Remark 1. The convergence [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] means that, for any continuous function ϕ on Ω, the following limit relation holds true:

lim ε→0 I Z ε I ϕ(y ε I ) = Ω z(x)ϕ(x) dx. ( 24
)
When applied to sequences (Z ε ) that satisfy I Z ε I 2 ≤ C where C is a constant independent on ε, this notion of convergence ensures the existence of some function z ∈ L 2 (Ω) such that, up to a sub-sequence, (Z ε ) converges to z (see Lemma 10.1 in 16 ). In view of the relation [START_REF] Eringen | Microcontinuum field theories II: Fluent media[END_REF], we note that we can replace in [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] the Dirac measure δ y ε I by δ y ε I,s or even δ y ε I+p,s . Indeed, for any continuous function ϕ, we have ϕ(y ε I+p,s )ϕ(y ε I,s ) = o(1). Remark 2. When it holds for any s ∈ {1, . . . , K}, the convergence of measures Z ε s z s is closely related to the notion of two-scale convergence developed by Nguetseng 29 and Allaire 6 . The discrete variable s plays the role of the fast variable. In that case, for any convex lower semi-continuous function Φ we have

lim inf ε→0 I 1 K K s=1 Φ(Z ε I,s ) ≥ Ω Φ 1 K K s=1 z s (x) dx. ( 25
)
see Lemma 3.1 in 11 .

Remark 3. If Z ε stands for a family (Z ε s ) K s=1 and if, for any s ∈ {1, . . . , K}, Z ε s z we simply write Z ε z as no confusion can arise. Equation (25) in the previous remark then reduces to lim inf ε→0 I

1 K K s=1 Φ(Z ε I,s ) ≥ Ω Φ (z(x)) dx.
The goal of homogenization is the determination of an elastic problem whose solution is the limit of the equilibrium states when ε tends to zero. A unique homogenized energy is expected to determine the effective elastic problem under the action of diverse external forces. The notion of Γ-convergence is one of the suitable tool for that aim. Definition 2. The sequence of energies E ε is said to Γ-converge to E if i) for any family U ε of displacements at the nodes of the structure,

U ε u =⇒ lim inf ε→0 E ε (U ε ) ≥ E(u)
ii) for any u, there exists a family U ε such that

U ε u and lim sup ε→0 E ε (U ε ) ≤ E(u).
We refer to 14,18 for the definition and properties of Γ-convergence. The two most important properties for our purpose are the following: (i) provided that sequences (U ε ) with bounded energy (E ε (U ε ) < M ) are relatively compact for the chosen convergence, then minimizers of E ε converge to a minimizer of E; (ii) provided that a given external force potential is a continuous perturbation of the elastic energy, one can add this potential to the elastic energy and recover the same potential in the limit energy. Let us assume that external forces F ε I are applied at the nodes y ε I,1 of the structure and are of the type

F ε I = ε 3 F (y ε I,1 )
where F is a continuous vector function. The associated potential satisfies

lim ε→0 I F ε I • U ε I,1 = lim ε→0 I F (y ε I,1 ) • U ε I,1 = lim ε→0 Ω F (x) d( I U ε I,1 δ y ε I (x)) = Ω F (x) • u(x) dx.
as soon as U ε u and therefore is a continuous perturbation. A consequence of the two aforementioned properties is that the equilibrium displacement family U ε under the action of F ε I,1 converges to a minimum of E(u) -Ω F (x) • u(x) dx. Of course, the same also holds for other nodes y ε I,s with s = 1. So Γ-convergence allows to study the homogenization of our structure without taking into account external forces. In that sense it follows closely the intuitive notion of effective behavior.

Note however that the action of more oscillating external forces as well as the action of forces concentrated along the boundary are not encompassed by the present results that must therefore be adapted to these cases. It is also,the case if the displacement is imposed on some nodes (on the boundary or not) of the structure. Therefore the particular problem in the section entitled "strain-gradient models in elasticity" and represented in Figure 1 does not rigorously enter the scope of our result: we have chosen it for illustrative purposes and the reader may be confident that sufficiently concentrated but diffuse forces at the vicinity of the points where the displacement is imposed would lead to solutions close to the ones presented in Figure 2.

Homogenization result

In this section we seek the Γ-limit E of the sequence of functionals E ε defined in [START_REF] Boutin | Microstructural effects in elastic composites[END_REF]. We consider free boundary conditions. In this case, the equilibrium of the structure can be reached only when the applied external forces are balanced and the solution of equilibrium problems is defined up to a global rigid motion. In order to ensure uniqueness of the equilibrium solution, we forbid a global rigid motion by assuming

I U ε I,1 = 0 and I θ ε I,1 = 0. ( 26
)
This condition, which in the sequel is referred to as the "zero mean rigid motion" assumption , is taken into account in the energy by setting E ε (U ε , θ ε ) = +∞ whenever it is not satisfied. We extend the functionals Ē and F defined in ( 21) and ( 22) onto whole L 2 (Ω) by setting Ē = +∞ or F = +∞ whenever the integrands are not square integrable. The homogenized energy corresponding to the microscopic energy E ε is provided by the following theorem.

Theorem 1. The sequence of energies E ε defined in (12) Γ-converges to the energy E defined by

E(u) := inf w,v,θ Ē(w, ξ u,v , v, η u ) + F (v, η u , θ) : Ē(v, η u , 0, 0) = 0, Ω θ 1 (x) dx = 0 ( 27
)
if Ω u(x) dx = 0, and E(u) := +∞ otherwise, where η u , ξ u,v are the quantities defined in [START_REF] Dell'isola | Large deformations of planar extensible beams and pantographic lattices: Heuristic homogenization, experimental and numerical examples of equilibrium[END_REF] and [START_REF] Dell'isola | The relationship between edge contact forces, double forces and interstitial working allowed by the principle of virtual power[END_REF]. Here the infimum is taken on all functions w, v and θ in L 2 (Ω). More precisely we have (i) (Compactness) For all sequences (U ε , θ ε ) with bounded energy (i.e. E ε (U ε , θ ε ) ≤ M ), there exist u and θ such that, up to a sub-sequence, U ε u and θ ε θ.

(ii) (Lower bound) For all sequences (U ε , θ ε ) such that U ε u and θ ε θ, we have

lim inf ε→0 E ε (U ε , θ ε ) ≥ E(u).
(iii) (Upper bound) For any u satisfying E(u) < +∞, there exists a sequence (U ε , θ ε ) such that

U ε u and lim sup ε→0 E ε (U ε , θ ε ) ≤ E(u).
The next subsections are devoted to the proof of this theorem.

Compactness of sequences with bounded energy

In this subsection, we prove point (i) of Theorem 1, that is the relative compactness of sequences (U ε , θ ε ) with bounded energy. We also prove the relative compactness of the associated sequences m ε , v ε , χ ε , ω ε p,s,s defined in ( 13) and ( 14). This compactness result is necessary to ensure that any sequence of minima of E ε (U ε , θ ε ) converges towards a minimum of E(u).

The next lemma is a discrete version of the Korn inequality. Its proof is essentially based on the classical Korn inequality applied to an affine interpolation of the nodes displacements.

Lemma 1. There exists a positive constant C depending only on Ω such that, for any sequence of families (U ε I ) I∈I ε , (θ ε I ) I∈I ε satisfying, for any p ∈ P,

I U ε I = 0, I θ ε I = 0, I U ε I+p -U ε I ε -θ ε I × p 2 < 1, ( 28
)
one has

I θ ε I 2 < C, I U ε I+p -U ε I ε < C and I U ε I 2 < C. ( 29
)
The choice of the node 1 in this assumption is arbitrary. A mean value of the displacements of all node could have been chosen instead. Rotation could have alternatively been forbidden by imposing a condition on the displacements.

Proof. Here C is a positive constant (depending only on Ω) whose value may vary from line to line.

We first divide Ω into ε -3 small parallelepipeds Ω ε I := y ε I + εΩ. Then we divide each Ω ε I into six tetrahedra by considering the six different triplets (p i , p j , p k ) of distinct directions in {t 1 , t 2 , t 3 } and the associated tetrahedra with vertices (y ε I , y ε I+pi , y ε I+pi+pj , y ε I+pi+pj +p k ). On each tetrahedron we define the affine interpolation U ε of U ε by setting for any 0

≤ z ≤ y ≤ x ≤ 1, U ε (y ε I + ε(xp i + yp j + zp k )) := (1 -x)U ε I + (x -y)U ε I+pi + (y -z)U ε I+pi+pj + zU ε I+pi+pj +p k .
Checking that this piecewise affine function belongs to H 1 (Ω) is straightforward. Moreover, there exists a constant C such that

U ε 2 L 2 (Ω) = ε 3 |det(t 1 , t 2 , t 3 )| I∈Iε 1 0 1 0 1 0 |U ε (y ε I + ε(xp i + yp j + zp k ))| 2 dx dy dz ≥ C I U ε I 2 . ( 30
)
In the tetrahedron (y ε I , y ε I+pi , y ε I+pi+pj , y ε I+pi+pj +p k ) we have

ε∇U ε = (U ε I+pi -U ε I ) ⊗ p * i + (U ε I+pi+pj -U ε I+pi ) ⊗ p * j + (U ε I+pi+pj +p k -U ε I+pi+pj ) ⊗ p * k , (31)
where (p * i , p * j , p * k ) is the dual basis corresponding to (p i , p j , p k ). Therefore, denoting by e(U ε ) the symmetric part of ∇U ε , the following equalities hold:

εe(U ε ) : (p i ⊗ p i ) = (U ε I+pi -U ε I -εθ ε I × p i ) • p i , εe(U ε ) : (p j ⊗ p j ) = (U ε I+pi+pj -U ε I+pi -εθ ε I × p j ) • p j , εe(U ε ) : (p k ⊗ p k ) = (U ε I+pi+pj +p k -U ε I+pi+pj -εθ ε I × p k ) • p k , εe(U ε ) : (p i ⊗ p j + p j ⊗ p i ) = (U ε I+pi -U ε I -εθ ε I × p i ) • p j + (U ε I+pi+pj -U ε I+pi -εθ ε I × p j ) • p i , εe(U ε ) : (p i ⊗ p k + p k ⊗ p i ) = (U ε I+pi -U ε I -εθ ε I × p i ) • p k + (U ε I+pi+pj +p k -U ε I+pi+pj -εθ ε I × p k ) • p i , εe(U ε ) : (p j ⊗ p k + p k ⊗ p j ) = (U ε I+pi+pj -U ε I+pi -εθ ε I × p j ) • p k + (U ε I+pi+pj +p k -U ε I+pi+pj -εθ ε I × p k ) • p j
, where we have used the identity (θ ε I × p) • q + (θ ε I × q) • p = 0 with p, q ∈ {p i , p j , p k }. As p i ⊗ p i , p j ⊗ p j , p k ⊗ p k , (p i ⊗ p j + p j ⊗ p i ), (p i ⊗ p k + p k ⊗ p i ), (p j ⊗ p k + p k ⊗ p j ) make a basis for symmetric matrices, then by using [START_REF] Jikov | Homogenization of Differential Operators and Integral Functionals[END_REF] and summing over all parallelepipeds Ω ε I , we obtain that e(U ε ) 2 L 2 (Ω) ≤ C. We can use classical Korn inequality on Ω (see Theorem 2.5 in 30 ): there exist a constant C and a global rigid motion R ε (y) := a ε + B ε • y (where a ε is a constant vector and B ε is a skew-symmetric matrix) such that

U ε -R ε 2 H 1 (Ω) ≤ C. In particular we have ∇U ε -B ε 2 L 2 (Ω) ≤ C. (32)
From ( 31) and ( 32), we get for any direction p i , the inequality

I U ε I+pi -U ε I ε -B ε • p i 2 ≤ C (33)
which, together with assumption (28), gives using Jensen and triangle inequalities,

I B ε • p i -θ ε I • p i 2 ≤ I B ε • p i -θ ε I × p i 2 ≤ C. ( 34
)
We obtain

B ε • p i -I θ ε I • p i 2 ≤ C which, since I θ ε I = 0, reduces to B ε • p i 2 ≤ C.
Going back to [START_REF] Rahali | Homogenization à la Piola produces second gradient continuum models for linear pantographic lattices[END_REF] we obtain I θ ε I × p i 2 ≤ C. These two inequalities, being true for any direction p i , we deduce that

B ε 2 ≤ C and I θ ε I 2 ≤ C. ( 35
)
Using again the triangle inequality and inequality [START_REF] Placidi | A review on 2D models for the description of pantographic fabrics[END_REF], we obtain the second desired bound:

I U ε I+pi -U ε I ε 2 ≤ C. (36)
Prepared using sagej.cls

On the other hand, it is clear that U ε -R ε is the affine interpolation of U ε I -R ε • y ε I . We can then apply inequality (30) to these functions. This gives

C I U ε I -R ε • y ε I 2 ≤ U ε -R ε 2 L 2 (Ω)
and thus

I U ε I -R ε • y ε I 2 ≤ C.
Using again Jensen inequality and the assumption I U ε I = 0, we get R ε • ( I y ε I ) 2 ≤ C and thus R ε 2 ≤ C. This leads to the last desired bound:

I U ε I 2 ≤ C. (37)
This concludes the proof of the lemma.

Lemma 2. Compactness. Let (U ε , θ ε ) be a sequence of displacements satisfying E ε (U ε , θ ε ) ≤ C and with zero mean rigid motion as in [START_REF] Germain | Cours de Mécanique des Milieux Continus. Tome I: Théorie générale[END_REF]. Then the sequences

I U ε I,s 2 , I θ ε I,s 2 , I m ε I 2 , I v ε I,s 2 , I χ ε I,p
2 and I (ω ε I,p,s,s ) 2 are bounded.

Proof. Here C is a constant independent of ε whose value may vary from line to line.

The inequality

E ε (U ε , θ ε ) ≤ C implies that E ε (U ε ) ≤ C and F ε (U ε , θ ε ) ≤ C. The first bound E ε (U ε ) ≤ C leads directly to I (ω ε I,p,s,s ) 2 = ε -2 I U ε I+p,s -U ε I,s ε • τ ε p,s,s 2 ≤ C. (38)
We first consider any (p, s, s

) in A. The second bound F ε (U ε , θ ε ) ≤ C implies I θ ε I,s -τ ε p,s,s × U ε I+p,s -U ε I,s ε p,s,s 2 ≤ C (39)
and

I θ ε I+p,s -τ ε p,s,s × U ε I+p,s -U ε I,s ε p,s,s 2 ≤ C. (40)
Using the triangle inequality, then (39) and (40) yield

I θ ε I+p,s -θ ε I,s 2 ≤ C.
Now let us consider any (p, s, s ) in P × {1, . . . , K} 2 . Owing to the connectedness assumption we consider the path (s i , p i , i ) i=1,...,r connecting nodes (I, s) and (I + p, s ). We set p j := j-1 i=1 i p i . Using the triangle inequality, we get

I θ ε I+ pj ,sj -θ ε I,s 2 ≤ C. (41)
Hence, from (39) or (40), we can deduce

I j U ε I+ pj+1,sj+1 -U ε I+ pj ,sj ε pj ,sj ,sj+1 × τ ε pj ,sj ,sj+1 -θ ε I,s 2 ≤ C. (42)
Let us set U ε I,s,J,s := U ε J,sθ ε I,s × (y ε J,sy ε I,s ) so that

U ε I,s,J+p,s -U ε I,s,J,s = U ε J+p,s -U ε I,s -ε p,s ,s θ ε I,s × τ ε p,s ,s .
Then inequality (42) yields

I U ε I,s,I+ pj+1,sj+1 -U ε I,s,I+ pj ,sj ε pj ,sj ,sj+1 × τ ε pj ,sj ,sj+1 2 ≤ C.
Since (38) also implies that

ε -2 I U ε I,s,I+ pj+1,sj+1 -U ε I,s,I+ pj ,sj • τ ε pj ,sj ,sj+1 2 ≤ C, we get ε -2 I U ε I,s,I+ pj+1,sj+1 -U ε I,s,I+ pj+1,sj 2
< C which leads, still using the triangle inequality, to

ε -2 I U ε I,s,I+p,s -U ε I,s,I,s 2 ≤ C,
or equivalently to

ε -2 I U ε I+p,s -U ε I,s -θ ε I,s × (y ε I+p,s -y ε I,s ) 2 ≤ C. (43)
We focus temporarily on the particular case s = s = 1 which reads

I U ε I+p,1 -U ε I,1 ε -θ ε I,1 × p 2 < C. (44)
Applying Lemma 1 we obtain I θ ε I,1

2 ≤ C. Owing to (41) this bound extends to any s ∈ {1, . . . , K}:

I θ ε I,s 2 ≤ C. (45)
Using inequality (43), we obtain for any (p, s, s ) ∈ P × {1, . . . , K}

U ε I+p,s -U ε I,s 2 ≤ C. 2 , ε -2 I
As Lemma 1 also gives I U ε I,1

2 ≤ C, we get for any s ∈ {1, . . . , K}, by choosing p = 0 and s = 1,

I U ε I,s 2 ≤ C. (47)
Taking the mean value with respect to s in (47) and (46) (with p = 0) we get respectively

I m ε I 2 ≤ C, I v ε I,s 2 ≤ C, (48)
and taking the mean value with respect to s and s in (46), we get

I χ ε I,p 2 ≤ C. ( 49
)
The proof is concluded by collecting ( 38), ( 45), ( 47), ( 48) and (49).

Properties of the double-scale limits

The estimates established in Lemma 2 imply that for any sequence (U ε , θ ε ) with bounded energy, there exist functions θ, u, v, χ p and ω p,s,s in L 2 (Ω) such that, for any s ∈ {1, . . . , K} and up to sub-sequences,

θ ε s θ s , m ε u, v ε s v s , χ ε p χ p and ω ε p,s,s ω p,s,s . (50)
In the following lemma we establish some properties of these limits.

Lemma 3. Let (U ε , θ ε ) be a sequence such that E ε (U ε , θ ε ) ≤ C. We have U ε s u, Ω u dx = 0, Ω θ 1 dx = 0, K s=1 v s = 0 and χ p = ∇u • p. (51)
Moreover there exist some fields w s and λ in L 2 (R 3 ) such that, for any (p, s, s ) ∈ A,

ω p,s,s = w s -w s + ∇(v s + λ) • p + 1 2 ∇∇u : (p ⊗ p) • τ p,s,s + (v s -v s + ∇u • p) • ρ p,s,s . (52)
Prepared using sagej.cls

Proof. The convergence of v ε s v s implies that (U ε sm ε ) 0 and so, for any s, U ε s u. The assumptions I U ε I,1 = 0 and I θ ε I,1 = 0 also hold in the limit and give Ω u dx = 0 and Ω θ 1 dx = 0. The fact that K s=1 v ε I,s = 0 clearly implies that K s=1 v s (x) = 0. For any smooth test field ϕ with compact support in Ω, we have

Ω χ p (x) • ϕ(x) dx = lim ε→0 I ε -1 (m ε I+p -m ε I ) • ϕ(y ε I ) = lim ε→0 I m ε I • ε -1 (ϕ(y ε I-p ) -ϕ(y ε I )) = lim ε→0 I m ε I • (-∇ϕ(y ε I ) • p + O(ε)) = - Ω u(x) • (∇ϕ(x) • p) dx = Ω (∇u(x) • p) • ϕ(x) dx.
We deduce that χ p = ∇u • p. In order to characterize the limit ω p,s,s , we introduce the set D A of families of distributions in H -1 (R

ε -2 (U ε I+p,s -U ε I,s ) • (φ p,s,s (y ε I ) τ ε p,s,s ).
On the other hand, using functions φ p,s,s satisfying (53) and (55) and using the decomposition (7) of τ ε p,s,s , we obtain

I (p,s,s )∈A ε -2 (U ε I+p,s -U ε I,s ) • (φ p,s,s (y ε I ) τ ε p,s,s ) = I (p,s,s )∈A (ε -1 (v ε I+p,s -v ε I,s ) + ε -2 (m ε I+p -m ε I )) • (φ p,s,s (y ε I ) τ ε p,s,s ) = I (p,s,s )∈A (ε -1 (v ε I+p,s -v ε I,s ) + ε -2 (m ε I+p -m ε I )) • (φ p,s,s (y ε I ) τ p,s,s ) + I (p,s,s )∈A (v ε I+p,s -v ε I,s + ε -1 (m ε I+p -m ε I )) • (φ p,s,s (y ε I ) ρ p,s,s ) + I (p,s,s )∈A ε(v ε I+p,s -v ε I,s + ε -1 (m ε I+p -m ε I )) • (φ p,s,s (y ε I ) ζ ε p,s,s ).
It is clear that the third addend in the last sum tends to zero, since ζ ε p,s,s is bounded and

v ε I+p,s -v ε I,s + ε -1 (m ε I+p -m ε I ) tends to v s (x) -v s (x) + ∇u(x) • p as ε → 0.
Let us estimate separately the remaining two addends. To that aim we restrict ourselves to smooth functions φ p,s,s with compact support in Ω. Owing to Remark 1, we get lim ε→0 I (p,s,s )∈A The previous results prove that the distribution

(v ε I+p,s -v ε I,s + ε -1 (m ε I+p -m ε I )) • (φ p,
(ε -1 (v ε I+p,s -v ε I,s ) + ε -2 (m ε I+p -m ε I )) • (φ p,s,s (y ε I ) τ p,s,s ) = (p,s,s )∈A (ε -1 (v ε I+p,s -v ε I,s ) + ε -1 (v ε I,s -v ε I,s ) + ε -2 (m ε I+p -m ε I )) • (φ p,s,s (y ε I ) τ p,s,s ) = (p,s,s )∈A ε -1 (v ε I+p,s -v ε I,s ) • (φ p,s,s (y ε I ) τ p,s,s ) + (p,s,s )∈A ε -2 (m ε I+p -m ε I ) • (φ p,
ω p,s,s -∇v s • p + 1 2 ∇∇u : (p ⊗ p) • τ p,s,s -(v s -v s + ∇u • p) • ρ p,s,s
is orthogonal to all smooth functions in D A ⊥ with compact support in Ω. As such functions are dense in D A ⊥ , there exist some fields w s and λ in L 2 (R 3 ; R 3 ) such that, for any (p, s, s ) ∈ A,

ω p,s,s = w s -w s + ∇(v s + λ) • p + 1 2 ∇∇u : (p ⊗ p) • τ p,s,s + (v s -v s + ∇u • p) • ρ p,s,s .
This concludes the proof of the lemma.

Prepared using sagej.cls

Proof of the homogenization result (i) Lower bound: Let (U ε , θ ε ) be a sequence such that U ε u and θ ε θ. If E ε (U ε , θ ε ) = +∞, there is nothing to prove. We then assume that E ε (U ε , θ ε ) ≤ C. From Lemma 3 we know that the constraints Ω θ 1 (x) dx = 0 and Ω u(x) dx = 0 are satisfied. From the same lemma we know that v ε s converges to some v s and that χ ε p η u . As E ε (U ε ) is bounded, then ε 2 E ε (U ε ) tends to zero. By virtue of Remark 2, we have

0 = lim inf ε→0 (ε 2 E ε (U ε )) = lim inf ε→0 (ε 2 Ēε (v ε , χ ε )) ≥ Ē(v, η u , 0, 0)
which gives the constraint Ē(v, η u , 0, 0) = 0. Definition ( 14) allows us to rewrite E ε (U ε ) as I (p,s,s )∈A a p,s,s 2 (ω ε I,p,s,s ) 2 . Therefore, using (50), (51), (52), Remarks 1 and 2, we get

lim inf ε→0 ( Ēε (v ε , χ ε ) + Fε (v ε , χ ε , θ ε )) ≥ Ē(w, ξ u,v+λ , v, η u ) + F (v, η u , θ) ≥ inf w,v,θ { Ē(w, ξ u,v+λ , v, η u ) + F (v, η u , θ)}. (56)
We obtain the desired inequality by noticing that F (vλ, η u , θ) = F (v, η u , θ) and Ē(v + λ, η u , 0, 0) = Ē(v, η u , 0, 0) = 0.

(ii) Upper bound: By a density argument, it is enough to consider a function u ∈ C ∞ (Ω) such that E(u) < +∞.

In virtue of coercivity and lower semi-continuity of the functionals Ē and F , we can introduce the fields (v, w, θ) belonging to

C ∞ (Ω) such that E(u) = Ē(w, ξ u,v , v, η u ) + F (v, η u , θ), Ē(v, η u , 0, 0) = 0, Ω u(x) dx = 0
and Ω θ 1 (x) dx = 0. As the constraint Ē(v, η u , 0, 0) = 0 is equivalent to equation ( 18), we can apply the results of the asymptotic expansion which state that

lim ε→0 E ε (U ε , θ ε ) = Ē(w, ξ u,v , v, η u ) + F (v, η u , θ) = E(u)
where U ε and θ ε are defined by U ε I,s := u(y ε I ) + εv s (y ε I ) + ε 2 w s (y ε I ) and θ ε I,s := θ s (y ε I ).

(57)

However we must modify such families in order to fulfill the constraints I U ε I,1 = 0 and I θ ε I,1 = 0. This is easy to do: as the functions u, v 1 , w 1 , θ 1 are smooth and satisfy Ω u(x) dx = 0 and Ω θ 1 (x) dx = 0, the quantities

M ε := I u(y ε I ) + εv 1 (y ε I ) + ε 2 w 1 (y ε I ) and t ε := I θ 1 (y ε I )
are O(ε). We now redefine U ε and θ ε by setting

U ε I,s := u(y ε I ) + εv s (y ε I ) + ε 2 w s (y ε I ) -M ε and θ ε I,s := θ s (y ε I ) -t ε . (58)
It is clear that I U ε I,1 = I θ ε I,1 = 0 and that the limits U ε u and θ ε s θ s still hold. As the energy E ε (U ε , θ ε ) is invariant by the addition of a rigid motion we still have lim

ε→0 E ε (U ε , θ ε ) = E(u).

Explicit computation of the homogenized energy

Cell problems

In the homogenized energy we obtained, i.e.,

E(u)

:= inf w,v,θ { Ē(w, ξ u,v , v, η u ) + F (v, η u , θ) : Ē(v, η u , 0, 0) = 0},
one has to compute the minimum with respect to three extra kinematic variables w, v, θ. Let us explain how that can be done.

We first remark that no partial derivative of w appears in the expression of Ē(w, ξ u,v , v, η u ). Thus the minimization with respect to w is a local cell problem in which ξ u,v , v, η u are three vectorial parameters. The linear problem associated with the finite dimensional unknown w can easily be solved. Similarly, no partial derivative of θ appears in the expression of F (v, η u , θ). The minimization with respect to θ is a local cell problem in which η u is a vectorial parameter. Again the linear problem associated with the finite dimensional unknown θ can easily be solved.

Minimizing with respect to v is more subtle. The constraint Ē(v, η u , 0, 0) = 0 imposes v to minimize Ē(v, η u , 0, 0) and this problem can be solved like the two first ones. The unknown v depends linearly on the "parameter" η u (or equivalently on ∇u) up to an element of the kernel of the functional Ē(v, η u , 0, 0). Indeed, this kernel is never trivial as it contains at least the uniform fields (i.e. the families (v s ) which do not depend on s) but it may be much larger. Hence we get a linear operator L such that v = L • ∇u + λ, where λ is any field in (L 2 (Ω; R 3 )) K with values in the kernel of Ē(v, η u , 0, 0). A first consequence is that the constraint Ē(v, η u , 0, 0) = 0 leads to the following constraint on ∇u:

Ē(L • ∇u, η u , 0, 0) = 0

Then we can replace w, θ and v by the solutions of these cell problems. The definition [START_REF] Dell'isola | Large deformations of planar extensible beams and pantographic lattices: Heuristic homogenization, experimental and numerical examples of equilibrium[END_REF] of η u and the dependence of ξ u,v on ∇v through [START_REF] Dell'isola | The relationship between edge contact forces, double forces and interstitial working allowed by the principle of virtual power[END_REF] imply that the homogenized energy can be written as a quadratic functional of ∇u, ∇∇u, λ and ∇λ or equivalently of e(u), ∇e(u), λ and ∇λ:

E(u) = inf λ∈L 2 (Ω) Ω Q(e(u), ∇e(u), λ, ∇λ), ( 59
)
where Q is a non negative quadratic form. An extra kinematic variable λ still remains in the expression (59) but the minimization with respect to this variable cannot be performed locally as its gradient is involved. Indeed, many examples provided in 3 show that the effective behavior of frame lattices may correspond to materials which must be described in terms of extra kinematic variables (e.g. Cosserat, Timoshenko or Reissner models). Then, in general the limit model obtained is both a generalized continuum model and a second gradient model 3 . The variable λ can be interpreted as a microadjustement.

A general algorithm for making explicit the homogenized energy E(u) has been given in 3 for the particular case where z s = 0. The procedure is quite similar in the more general case that we consider here. Only the canonical form of the functional energy Ē has to be modified as the expression of the limit extensional energy Ē obtained here is different from the one obtained in 3 .

Back to the asymmetrical pantographic structure

The asymmetrical structure described in the third section enters our general description of lattices. The family (y i ) is the one defined in [START_REF] Alibert | Truss modular beams with deformation energy depending on higher displacement gradients[END_REF]. Family (z i ) is z 1 = z 2 = z 3 = z 6 = 0, z 4 = z 5 = (α, 0, 0) .

(60)

Considering the mechanical interactions which were described and drawn in Figures are the only non vanishing interaction coefficients.

To fix ideas, we consider the case where all bars are made of the same material with Poisson coefficient ν = 0.3 and have the same circular section with radius ε 2 . Hence the non vanishing coefficients a p,s,r are proportional to -1 p,s,r . A simple choice of energy unit allows us to fix a p,s,r = -1 p,s,r . We have β = 1, thus bending and torsional rigidities are fixed. Applying the procedure described in the previous section is straightforward but needs the help of a computer at least for the determination of bending and torsion parts. Minimization with respect to λ is obtained in this particular case with λ = 0. So no extra kinematic variable is needed and the effective energy (59) reduces to the form [START_REF] Abdoul-Anziz | Homogénéisation de composites élastiques périodiques à fort contraste: Conception de métamatériaux de second gradient[END_REF]. More specifically we get the homogenized elastic energy

Conclusions

In this article we have provided a new rigorous homogenization result for periodic structures made of welded elastic bars. We have shown that allowing the positions of the nodes of the structure inside the rescaled periodic cell to depend on the small parameter ε, may lead to a strain-gradient model in which strain-gradient and classical strain terms are coupled. We have illustrated the astonishing effects of such a coupling on the equilibrium solutions. Note that, for this coupling to arise, the nodes positions usually assumed to be y (i,j,k),s := ε(y s + it 1 + jt 2 + kt 3 ) are modified at the order ε 2 within each cell.

To summarize the results of the present study and those of 1 , let consider a generic modification ε α z s of the positions of the nodes within the periodic lattice.

-If α = 1, then y s is simply replaced by y s + z s . The new geometry can still lead to a homogenized uncoupled strain-gradient model but, in the generic case, the structure will become rigid (see 1 ). -If 1 < α < 2 then, in the generic case, the structure will again become rigid.

-When α = 2, as in the present study, then an effective energy with a coupling between classical strain and strain-gradient terms can be obtained. -A modification with α > 2 will have no effect.

These results may explain how difficult it is to build in practice and characterize experimentally an efficient microarchitectured structure leading to a strain-gradient material. Indeed they show how sensitive is the effective model to the design of the lattice and to manufacturing defects.

Since second gradient effects are very difficult to measure experimentally, in order to get an experimental evidence of them it would be interesting to: (i) seek optimal topologies for frame structures, (ii) extend our study to nonlinear elasticity to take into account geometrical non-linearities that may arise in micro-structures. Some works in these directions can be found in the literature 7,8 but, up to our knowledge, in the context of nonlinear elasticity, second gradient materials have only been obtained heuristically 19 .

Figure 1 .Figure 2 .

 12 Figure 1. Imposing an extension on the left half of a bar.

1 Figure 3 .

 13 Figure 3. Standard pantographic beam

Figure 4 .Figure 5 .

 45 Figure 4. Perspective view of the 3D structure

Figure 6 .

 6 Figure 6. Asymmetrical pantographic structure. Nodes 4 and 5 have been slightly moved on the right from their previous positions (here in red) of Figure 3.

Figure 7 .

 7 Figure 7. Two points of view for describing the same lattice: In figure (a), we have a "small cell" with long-range interactions while in figure (b) we have chosen the cell large enough for our "next cell hypothesis" to be satisfied.

E(u) = 1 2 Ωu 1 ∂x 2 1 - 2 + ∂ 2 u 2 ∂x 2 1 2 dx

 2122 (e 13 (u)) 2 + σ (e 12 (u)) 2 + γ (e 11 (u)) 2 + κ( ∂ 2 8α(e 11 (u)) while the constraint inf v Ē(v, η u , 0, 0) = 0 reads e 22 (u) = e 33 (u) = e 23 (u) = 0. The effective material coefficients are ≈ 4.3, σ ≈ 7.1, γ ≈ 307.5, κ = 12 47 .

  3 ) defined byD A := ψ p,s,s = (w sw s + ∇λ • p) • τ p,s,s such that (p, s, s ) ∈ A, w s ∈ L 2 (R 3 ), λ ∈ L 2 (R 3 ) . orthogonal D A ⊥ isthe set of families (φ p,s,s ) (p,s,s )∈A of functions in H 1 (R 3 ) such that, for all ψ p,s,s ∈ D A , (p,s,s )∈A ψ p,s,s , φ p,s,s = 0. Let us remark that, for any φ ∈ D A ⊥ we have = (w s ) s=1,...,K ∈ (L 2 (R 3 ; R 3 )) K ,

	Its (p,s,s )∈A	(∇φ p,s,s • p) τ p,s,s = 0,	(53)
	and for any w (p,s,s )∈A	((w s -w s ) • τ p,s,s ) φ p,s,s = 0.	(54)
	We extend φ by setting φ p,s,s = 0 whenever (p, s, s ) / ∈ A. Then we can rewrite the equation (54) as
		(p,s,s )	(τ p,s,s φ p,s,s -τ p,s ,s φ p,s ,s ) = 0.	(55)
	On the one hand, we have
	Ω (p,s,s )∈A	ω p,s,s (x) φ p,s,s (x) dx = lim ε→0	I (p,s,s )∈A

  s (x)v s (x) + ∇u(x) • p, (φ p,s,s (x) ρ p,s,s ) .

	s,s (y ε I ) ρ p,s,s )
	=
	(p,s,s )∈A
	Using (54), the first addend becomes
	(p,s,s )∈A

v

  s,s (y ε I ) τ p,s,s ).

	The first term yields
	lim ε→0	I (p,s,s )∈A	ε -1 (v ε I+p,s -v ε I,s ) • (φ p,s,s (y ε I ) τ p,s,s )
							= lim ε→0	I (p,s,s )∈A	v ε I,s • (ε -1 (φ p,s,s (y ε I-p ) -φ p,s,s (y ε I )) τ p,s,s )
							= lim ε→0	I (p,s,s )∈A	v ε I,s • (-∇φ p,s,s (y ε I ) • p) τ p,s,s + O(ε)
							=	(p,s,s )∈A	v s (x), ((-∇φ p,s,s (x) • p) τ p,s,s )
							=	(p,s,s )∈A	∇v s (x) • p, (φ p,s,s (x) τ p,s,s ) .
	For the second term, using (53) we get
	lim ε→0	I (p,s,s )∈A	ε -2 (m ε I+p -m ε I ) • (φ p,s,s (y ε I ) τ p,s,s )
			= lim ε→0	I (p,s,s )∈A	ε -2 m ε I • (φ p,s,s (y ε I-p ) -φ p,s,s (y ε I )) τ p,s,s
			= lim ε→0	I (p,s,s )∈A	m ε I •	-ε -1 ∇φ p,s,s (y ε I ) • p +	1 2	∇∇φ p,s,s (y ε I ) : (p ⊗ p) τ p,s,s + O(ε)
			= lim ε→0	I (p,s,s )∈A	m ε I •	1 2	∇∇φ p,s,s (y ε I ) : (p ⊗ p) τ p,s,s + O(ε)
			=	(p,s,s )∈A	u(x),	1 2	∇∇φ p,s,s (x) • p τ p,s,s
			=	(p,s,s )∈A	1 2	∇∇u(x) : (p ⊗ p), (φ p,s,s (x)τ p,s,s ) .

  4-5-6 implies that coefficients a 1,1,2 , a 1,1,3 , a 1,1,4 , a 1,1,5 , a 1,2,5 , a 1,3,4 , a 1,4,6 , a 1,5,6 , a 2,4,3 , a 2,5,2 , a 2,6,2 , a 2,6,3 , a 3,1,1 , a 6,1,1 , a 9,1,1

Abdoul-Anziz & al

Prepared using sagej.cls

§ For nodes (I, s) and (I + p, s ) lying sufficiently inside Ω, the joining path which is assumed to be independent of I, lies inside Ω and connectedness is thus assured. Some nodes lying near the boundary may remain disconnected from the structure. In that case, such nodes must be deleted otherwise compactness could not be assured.Prepared using sagej.cls

¶ Note that (ηu)p,s actually does not depend on s. It is introduced in this way because the expression of the limit energy is simpler when ηu and ξu,v share the same tensorial nature.

Acknowledgments

The authors thank the French Région Sud and C.N.R.S (Projets Exploratoires Premier Soutien) for their financial support.


Yiyun Chen

Florian Jabbour
email: florianjabbour@mnhn.fr

Andriy Novikov

Wei Wang

Sylvain Gerber

A study of floral shape variation in Delphinieae (Ranunculaceae) using geometric morphometrics on herbarium specimens

Keywords: Aconitum, Delphinium, Flower morphology, Nectar spur, Staphisagria, Taxonomy

The variation in the shape of flowers, the reproductive structures of angiosperms, is generally investigated in a qualitative way, or using multivariate statistical analyses of distance data (i.e. traditional morphometrics). In this study, we evaluate the application of geometric morphometric methods to flowers of herbarium specimens, a material that is still being massively overlooked as far as these methods are concerned. Flowers of Delphinieae (Ranunculaceae) are synorganised: the spurred and nectariferous dorsal petals are nested within the spur or hood of the dorsal sepal. The study of synorganised structures is challenging, and geometric morphometric-oriented investigation may provide working hypotheses for an evo-devo exploration of such hyperorgan formation. We analyze here the floral shape of Delphinieae using geometric morphometric methods and demonstrate the effectiveness of these methods on dry and compressed flowers. We also show the diagnostic value of the dorsal half of the perianth of the Delphinieae flower, and suggest some trends to help design future studies on the evo-devo of spur formation and the co-evolution of Delphinieae flowers with their pollinators.

Introduction

The extent of shape variation in angiosperm reproductive structures is remarkable [START_REF] Endress | Diversity and evolutionary biology of tropical flowers[END_REF][START_REF] Endress | Symmetry in flowers: diversity and evolution[END_REF][START_REF] Sauquet | The ancestral flower of angiosperms and its early diversification[END_REF]. Most botanical studies aiming at quantifying the variation in floral shape have so far employed traditional morphometrics, the multivariate statistical analysis of distance data (e.g. [START_REF] Nagahama | Taxon delimitation in the Andropogon lateralis complex (Poaceae) in southern South America based on morphometrical analyses[END_REF][START_REF] Chartier | The evolution of afro-montane Delphinium (Ranunculaceae): Morphospecies, phylogenetics and biogeography[END_REF], while modern geometric morphometric methods (GMM) have remained surprisingly underused. However, the development of geometric morphometrics since the nineties, based on the spatial relationships among homologous anatomical landmarks, has revolutionized the study of biological shape variation [START_REF] Rohlf | A revolution morphometrics[END_REF]. The main advantage of geometric morphometrics over traditional morphometrics is that the geometry of the biological structure under study (as captured by the relative positions of landmarks) is preserved throughout the analyses, allowing for insightful visualizations and refined interpretations of morphological patterns (Adams, Rohlf, and Slice2004).

Among the botanical studies that have applied geometric morphometrics, some have used landmark data to study shape variation in fresh plant material in two dimensions (e.g. leaves:

Jensen, Ciofani, and Miramontes 2002; or seeds: Jacques and Zhou 2010) or fresh flowers in 2D [START_REF] Gómez | Natural selection on Erysimum mediohispanicum flower shape: insights into the evolution of zygomorphy[END_REF][START_REF] Abdelaziz | Using complementary techniques to distinguish cryptic species: a new Erysimum (Brassicaceae) species from North Africa[END_REF][START_REF] Savriama | Geometric morphometrics of corolla shape: dissecting components of symmetric and asymmetric variation in Erysimum mediohispanicum (Brassicaceae)[END_REF][START_REF] Fernández-Mazuecos | Corolla morphology influences diversification rates in bifid toadflaxes (Linaria sect. Versicolores)[END_REF][START_REF] Gardner | Characterizing floral symmetry in the Core Goodeniaceae with geometric morphometrics[END_REF] or 3D [START_REF] Van Der Niet | Three-dimensional geometric morphometrics for studying floral shape variation[END_REF][START_REF] Wang | Quantifying floral shape variation in 3D using microcomputed tomography: a case study of a hybrid line between actinomorphic and zygomorphic flowers[END_REF]. However, fresh plant material can be difficult to gather and conserve over the time of the study. Here, we evaluated the application of GMM to dry and compressed flowers of herbarium specimens, a material that is still being massively overlooked (but see [START_REF] Shipunov | Geometric morphometrics as a tool for understanding Dactylorhiza (Orchidaceae) diversity in European Russia[END_REF], where GMM were applied to flowers dried for the purpose of the study) compared with dry leaves from herbarium specimens [START_REF] Viscosi | Leaf morphology, taxonomy and geometric morphometrics: a simplified protocol for beginners[END_REF][START_REF] Cope | Plant species identification using digital morphometrics: A review[END_REF].

Synorganisation is defined as the intimate structural connection of two or several neighboring floral organs to form a functioning system or apparatus ("hyperorgan"). It may occur through congenital or postgenital fusion or by coordination of organs without fusion [START_REF] Specht | Flower evolution: the origin and subsequent diversification of the angiosperm flower[END_REF][START_REF] Endress | Synorganisation without organ fusion in the flowers of Geranium robertianum (Geraniaceae) and its not so trivial obdiplostemony[END_REF][START_REF] Endress | Development and evolution of extreme synorganization in angiosperm flowers and diversity: a comparison of Apocynaceae and Orchidaceae[END_REF]. Synorganisation is rare in basal angiosperms, but becomes an important trait in many monocot and eudicot flowers [START_REF] Specht | Flower evolution: the origin and subsequent diversification of the angiosperm flower[END_REF], such as orchids [START_REF] Vogel | Organographie der Blüten kapländischer Ophrydeen[END_REF][START_REF] Rudall | Roles of synorganisation, zygomorphy and heterotopy in floral evolution: the gynostemium and labellum of orchids and other lilioid monocots[END_REF]. It is hypothesized to be a key process in floral diversification [START_REF] Endress | Angiosperm floral evolution: morphological developmental framework[END_REF]). The study of synorganised structures is challenging, and GM-oriented investigation may provide working hypotheses for an evo-devo exploration of such hyperorgan formation.

Although morphologically diverse, the flowers of Delphinieae Warm. (Ranunculaceae) are all based on a single relatively simple body plan. The calyx is petaloid with a dorsal sepal deformed into a pocket (spur or hood). The corolla is reduced to a single petal (nectariferous; in Consolida and Aconitella), a pair of petals (both nectariferous; in Aconitum and Gymnaconitum), or four petals (the two dorsal-most being nectariferous; in Delphinium and Staphisagria), located in the dorsal half of the flower. The nectariferous petals are made of three components: stalk, spur, and labium. The proportions of these three structures vary among the major clades of Delphinieae (e.g. the petal of Aconitum has a long stalk with a short spur and labium, whereas the petal of Delphinium has a very short stalk, a relatively long spur, and a short labium). The nectariferous petals are enclosed into the pocket of the dorsal sepal, which makes the dorsal part of Delphinieae flowers synorganised (Jabbour and Renner 2012a).

A (traditional) morphometric study of Delphinium floral shape using herbarium specimens was first conducted by [START_REF] Chartier | The evolution of afro-montane Delphinium (Ranunculaceae): Morphospecies, phylogenetics and biogeography[END_REF], who measured 13 floral traits of afro-montane species and circumscribed three morphospecies, namely D. dasycaulon, D. leroyi and D. macrocentrum. We analyze here the floral shape of Delphinieae using GMM, in order to 1) assess the usefulness and effectiveness of geometric morphometrics on dry and compressed flowers from herbarium specimens, and 2) evaluate the diagnostic value of the dorsal module of the perianth (dorsal sepal + one dorsal petal) of the Delphinieae flower. Our preliminary study will help design future studies on the evo-devo of spur formation and the co-evolution of Delphinieae flowers with their pollinators.

Materials and Methods

Definition of categories of flower morphologies and taxonomic sampling

Delphinieae flowers can be classified into three morphological types: Aconitum, Delphinium and Consolida. Aconitum L., Delphinium L., and Consolida (DC.) Gray are three genera from the tribe Delphinieae (Fig. 1) that show contrasted floral morphologies, although they share the same floral groundplan (five petaloid sepals, the dorsal one enclosing the spurred petal(s) of the reduced corolla). The other three genera of Delphinieae can be associated with one of these three morphological types, as far as floral morphology is concerned.

In flowers of the type Aconitum, the dorsal sepal is hood-shaped or nightcap-shaped (Jabbour and Renner 2012a). It contains two dorsal petals, being usually hammer shaped with a long grooved stalk and a spurred limb [START_REF] Kosuge | Morphology of the petal in Aconitum[END_REF]. The limb consists of two parts, a labium and a spur. The spur is projected upward inside the upper hollowed sepal, often elongated and curved or coiled, but sometimes short and saccate [START_REF] Kosuge | Morphology of the petal in Aconitum[END_REF].

In flowers of the type Delphinium, the dorsal sepal is spurred and the corolla consists of four petals. The two dorsal nectariferous ones extend into the sepal spur, and the lateral ones, usually flat, cover the entrance of the spur. In flowers of the type Consolida, the sepal is spurred and the corolla consists of a single spurred petal that results from the postgenital fusion of the two dorsal-most petal primordia. The six other petal primordia have an arrested development (Jabbour and Renner 2012a) We follow here the same classification of Delphinieae as in Jabbour and Renner (2012b). The number of species sampled and the total number of species in each infratribal group is given between parentheses in Fig. 1.

Digitization of herbarium specimens and selection of the photographs

The specimens used in this study are from P (Herbier National, Muséum national d'Histoire naturelle, Paris) and KRA (Jagellonian University, Poland) herbaria. All photographs were taken by one of us (AN) based on an experimental protocol using a slab of LEDs that allows to trans-illuminate the herbarium specimens and obverse internal (i.e. hidden) floral structures (Fig. 2). All pictures were taken in a standardized way with a resolution of 3648×2736 pixels, and included a scale.

A total number of 2217 photographs were taken: 983 are from specimens with flowers of the type Aconitum and 1234 are from specimens with flowers of the type Delphinium (including the type Consolida, as Consolida and Aconitella specimens are often classified under Delphinium in the herbarium collections). Selecting a subset of photographs for our study was a two-step process. First, we selected the photographs that could be used for GMs, the socalled "good" photographs. A good photograph should meet two criteria: show clearly the internal structures and show a flower with as few folded structures as possible (Fig. 2).

After this first step, we kept 43% (423) and 38% (469) of the photos of the flowers from the type Aconitum and Delphinium, respectively. The second step consisted in selecting 20 photographs (corresponding to 20 species) for each of the three categories of flowers, choosing the taxa in order for the sampling to be representative of the main clades in the phylogeny of Delphinieae (Fig. 1 and Suppl. Table 1). After sorting the specimens that could be used for GM, our sampling included a single type: voucher R. Gombault 1717 (P00197360), holotype of Consolida gombaultii (J. Thiébaut) Munz.

Geometric morphometrics of dorsal sepal and petal

We focused on the dorsal sepal and one dorsal petal, since these are the structures producing, protecting, and limiting the access to the nectar. The remaining parts of the dried flowers are often too folded to be morphometrically reliable and too dark even when trans-illuminated. As most of the Delphinieae (except Consolida and Aconitella) have a pair of dorsal petals, and as they are both very similar in shape and size (FJ, pers. obs.), we decided to work on only one dorsal petal.

To characterize the shape of these floral structures, we used a set of nine landmarks (three for the dorsal sepal and six for the dorsal petal; Fig. 3 and Suppl. Table 2) and seven series of 20 to 25 semilandmarks to capture the homologous curved outlines of the sepal and the petal between the landmarks [START_REF] Gunz | Semilandmarks: a method for quantifying curves and surfaces[END_REF]. The coordinates of these landmarks and semilandmarks were manually recorded from the sixty photographs using the software tpsDig2 [START_REF] Rohlf | The tps series of software[END_REF]. We decided not to consider the stalk of the petal in the analyses, as this structure is very slender in Aconitum, making it difficult to place the semilandmarks between the landmarks 4 and 5, and 4 and 9.

Analyses at different levels were performed (sepal and petal considered together or separately), but a similar protocol was followed for each of them. Flowers from the Consolida type were not included in the comparison of the dorsal module shape among Delphinieae, as their single petal is not homologous to any of the two dorsal petals from flowers from the Aconitum or Delphinium types. Landmark configurations were aligned using partial generalized Procrustes analysis (GPA; [START_REF] Rohlf | The tps series of software[END_REF]Slice 1990, Dryden and[START_REF] Dryden | Statistical analysis of shape[END_REF] and semilandmarks were slid using the minimum bending energy criterion [START_REF] Bookstein | Morphometric tools for landmark data: geometry and biology[END_REF]) using tpsRelw [START_REF] Rohlf | The tps series of software[END_REF]. Procrustes analysis extracts shape variation by filtering out the effects of size, position and orientation (e.g., [START_REF] Mitteroecker | Advances in geometric morphometrics[END_REF]. Aligned coordinates were then orthogonally projected onto a linear, tangent shape space using the mean shape as point of tangency. All subsequent analyses were carried out on the tangent space coordinates.

The main features of variation were highlighted by principal component analysis (PCA; [START_REF] Jolliffe | Principal Component Analysis and Factor Analysis[END_REF]) and visually depicted by means of Thin-Plate Spline deformation grids [START_REF] Bookstein | Principal warps: Thin-plate splines and the decomposition of deformations[END_REF].

Results

Variation in the shape of the dorsal sepal + dorsal petal in each of the types Delphinium, Aconitum, and Consolida

Within each of the three types of dorsal floral module identified, namely Delphinium, Aconitum and Consolida, the quantitative capture of the shape of the dorsal sepal and of the dorsal petal allows to circumscribe well the different taxonomic subgroups (Fig. 4 and Supp.

Fig. 1). In the Delphinium type, the principal components PC1 and PC2 together explain c. 64% of the shape variation in the dorsal floral module (Fig. 4A). Towards positive scores on PC1, the inner and outer spurs get narrower and the labium of the dorsal sepal covers more closely the labium of the dorsal petal. Along PC2, the inner and outer spurs get narrower and the inner spur tend to occupy more fully the cavity of the sepal spur. The labium of the dorsal sepal gets also better adjusted to the labium of the dorsal petal. The representatives of the genus Staphisagria (samples 1 and 2) are separated from the rest of the samples (23 to 40).

The inner and outer spurs of Staphisagria flowers are wider, and the outer spur does not fit tightly around the inner spur. The rest of the samples belong to a grade of subgenera within the genus Delphinium: subg. Delphinium (samples 23 to 29), subg. Anthriscifolium (30), and subg. Delphinastrum (31 to 40). From the first cited subgenus to the third, the inner and outer spurs get wider and the labium of the dorsal sepal moves back from the labium of the dorsal petal. The sample corresponding to subg. Anthriscifolium (30) is morphologically closer from the samples belonging to subg. Delphinium than from subg. Delphinastrum. In the latter subgenus, samples are largely distributed across PC2, and hence show a high diversity of narrow and wide spurs, the outer one covering more or less tightly the inner one (Fig. 4A).

In the Aconitum type, the principal components PC1 and PC2 together explain c. 66% of the shape variation in the dorsal floral module (Supp. Fig. 1A). Towards positive scores on PC1, the dorsal sepal gets narrower and the limb of the dorsal sepal occupies a more vertical position within the hood. Along PC2, the spur and labium of the dorsal petal get narrower.

The representatives of the subgenus Lycoctonum (3 to 9) are separated from the specimens belonging to subgenus Aconitum (10 to 21). The single specimen of the monotypic genus Gymnaconitum (22) occupies an extreme position in the PCA.

In the Consolida type, the principal components PC1 and PC2 together explain c. 62% of the shape variation in the dorsal floral module (Supp. Fig. 1B). Towards positive scores on PC1, the inner spur tends to occupy less fully the space within the outer spur. Along PC2, both outer and inner spurs get wider. The representatives of Aconitella (57 to 60) are separated from the specimens identified as Consolida s. str. (41 to 56).

Variation in the shape of the dorsal sepal in Delphinieae

In Delphinieae as a whole, the quantitative capture of the shape of the dorsal sepal allows to circumscribe well several taxonomic subgroups (Fig. 4B). PC1 and PC2 together explain c. 72% of the shape variation in the dorsal sepal. Along PC1, the spur gets wider and becomes hood-shaped and the rim of the labium becomes concave. Along PC2, the very wide (hood-shape) spur of the dorsal sepal gets narrower. Representatives of A. subg. Aconitum (10 to 21) form a homogenous group that is different from the rest of Delphinieae. The samples belonging to A. subg. Lycoctonum (3 to 9) and the genus Aconitella (57 to 60) occupy in the PCA an intermediary position between A. subg. Aconitum (10 to 21) and the rest of Delphinieae. The samples of Consolida (41 to 56) cluster together with the samples from the three subgenera of Delphinium (23 to 40). The two specimens of Staphisagria (1 and 2) cluster with a specimen of Consolida (43, C. armeniaca) and all three occupy an extreme position in the PCA (Fig. 4B).

Variation in the shape of the dorsal petal in the Delphinium and Aconitum types

In the species belonging to the Aconitum and Delphinium types, shape variation of the dorsal petal allows to circumscribe well the genera and subgenera (Fig. 4C). PC1 and PC2 together explain c. 73% of the shape variation in the dorsal petal. Along PC1, the middle part of the limb gets bigger and the spur gets down-curved. Along PC2, the labium of the petal gets more compact. The subgenera A. subg. Aconitum and A. subg. Lycoctonum cluster together and are isolated from the rest of the Delphinieae on the PCA. All the genera and subgenera belonging to the Delphinium type form the second cluster. The sampled species of Staphisagria, the sister genus to all the other Delphinieae, occupy an intermediary position between the cluster of Aconitum species and that of Delphinium species (Fig. 4C).

Discussion

Geometric morphometrics of flowers from herbarium specimens

Our study shows that, provided that a selection of the appropriate material is made, GMM are applicable to flowers from herbarium specimens. In the best cases, flowers are dried in such positions to show the taxonomically important features. However, in most cases, flowers are dried in random positions, with many floral parts folded. Selecting appropriate specimens is therefore a crucial prerequisite for generating exploitable photographs in the context of a morphometric study. We had to restrict our sampling in such a way that most of the type specimens we photographed were cast aside. We retained a single type specimen in our final sampling, that of C. gombaultii (see Material and Methods).

Our analyses of floral shapes in the three types based on variation vectors (vectors showing the major direction and intensity of shape variation around the consensus shape, among samples) generated by the tpsRelw software based on the set of samples we studied (result not shown) reveal that the shape of the labium of the dorsal sepal and of the dorsal petal seems much more variable than the shape of the spurs. We interpret this as an artifact due to the compression of the plant material in the process of making herbarium specimens. For this reason, we focus our discussion on the results relative to spur shape.

We here tested the applicability of GMM of herbarium specimens in a non-destructive way.

While it could be interesting to compare our results with those obtained from fresh flowers, a study on herbarium and fresh leaves shows that, although the main conclusions were similar, the results were not strictly identical due to directional changes in shape changes, and that both kinds of material should not be considered together in a study [START_REF] Tomaszewski | Is shape of a fresh and dried leaf the same?[END_REF].

The dorsal module of Delphinieae flowers carries a strong phylogenetic signal

The quantitative description of the flowers of Delphinieae confirms that the shapes of the dorsal sepal and of the dorsal petal(s) are taxonomically discriminant traits. Indeed, the shape of the dorsal module of the flower captured using landmarks and semilandmarks was sufficient to circumscribe the major infratribal clades (genera or subgenera) within Delphinieae.

Our results suggest promising perspectives for the application of GMM in botanical taxonomic studies, as potentially new species can be placed in the space of the dorsal floral module shape. In this way, new species can be preliminarily attributed to a clade within Delphinieae, before confirming the result with morphological and molecular analyses.

In the Delphinium type, three main clusters could be identified, which correspond to the genus Staphisagria, the subgenus Delphinastrum, and the subgenera Delphinium+Anthriscifolium, respectively (Fig. 4A). The shape of the dorsal module of Staphisagria flowers is very different from all the other shapes found in species with the Delphinium type. This is in accordance with the results of molecular phylogeny, which moved the three Staphisagria species out of the genus Delphinium (Jabbour andRenner 2011b, 2012b). Similarly, in the Aconitum and Consolida types, GM applied on the dorsal module of the flower allowed circumscribing the genera and the infragenerical groups (Supp. Lycoctonum flowers. This morphological similarity and evolutionary convergence obviously inspired the botanist who coined the genus name Aconitella [START_REF] Spach | Genre Aconitelle. -Aconitella Spach[END_REF] and was highlighted by the one who subsequently updated the taxonomic treatment of the genus [START_REF] Soják | Aconitella Spach, eine vergessene Gattung der Familie Ranunculaceae[END_REF].

Implications for evo-devo studies

Spur initiation and development has been studied in some model species, such as Linaria (Plantaginaceae; [START_REF] Box | Characterization of Linaria KNOX genes suggests a role in petal-spur development[END_REF]) and Aquilegia (Ranunculaceae; [START_REF] Puzey | Evolution of spur-length diversity in Aquilegia petals is achieved solely through cellshape anisotropy[END_REF][START_REF] Yant | Molecular basis for threedimensional elaboration of the Aquilegia petal spur[END_REF]. In columbines, the three-dimensional elaboration of the spur is initiated by localized, oriented cell divisions surrounding the incipient nectary. This process creates a pouch that is extended by anisotropic cell elongation [START_REF] Puzey | Evolution of spur-length diversity in Aquilegia petals is achieved solely through cellshape anisotropy[END_REF]. [START_REF] Yant | Molecular basis for threedimensional elaboration of the Aquilegia petal spur[END_REF] hypothesized that the development of this evolutionary novelty could be promoted by nonmutually exclusive factors, including 1) prolonged, KNOX-dependent cell fate indeterminacy, 2) localized organ sculpting and/or 3) redeployment of hormone-signaling modules. Whether these three processes and anisotropic cell elongation are also at play during spur development in Delphinieae (belonging in Ranunculaceae as well) has still to be investigated. In addition, it would be interesting to study whether the molecular bases of sepal spur and petal spur development are similar or identical.

We observed that the inner spur(s) [the appendages of the dorsal petal(s)] differentially fit the space provided by the outer spur (the appendage of the dorsal sepal), depending on the species of Delphinieae considered. The development of the inner spur(s) can potentially constrain the development of the outer spur, and vice versa. We hypothesize that the constraints might be of genetic, hormonal, physical or mechanical origins. This landmark-based investigation of Delphinieae flower shape diversity allows classifying the floral shapes into quasi-discrete categories, wherein species of interest could be targeted for evo-devo studies of spur development.

First step towards investigating the evolution of Delphinieae flowers with their pollinators

Floral shape variation often reflects differences in plant pollination systems, particularly for plant species with specialized plant-pollinator interactions. Comparative studies have resulted in qualitative descriptions of these shape differences and multivariate analyses of floral trait measurements have, with varying degrees of success, shown that plant species can cluster in phenotype space according to general pollinator classes [START_REF] Wilson | A multivariate search for pollination syndromes among penstemons[END_REF][START_REF] Bröderbauer | The design of trapping devices in pollination traps of the genus Arum (Araceae) is related to insect type[END_REF]. However, in comparison to these traditional approaches, GMM have the potential to provide additional and more accurate insights into the associations between floral shape variation and pollination biology [START_REF] Van Der Niet | Three-dimensional geometric morphometrics for studying floral shape variation[END_REF][START_REF] Blanco-Pastor | Bees explain floral variation in a recent radiation of Linaria[END_REF]. Geometric morphometric analyses can enrich studies on the phenotypic evolution of flowers and of their pollinators. In the case of the spurred flowers of Delphinieae, for which the pollinators are expected to be specialists (but see [START_REF] Hollens | A spur-ious inference: pollination is not more specialized in long-spurred than in spurless species in Diascia-Rediviva mutualisms[END_REF][START_REF] Vlašánková | The nectar spur is not only a simple specialization for long-proboscid pollinators[END_REF], comparing the form of the dorsal floral module with the front of the pollinator's body would allow testing the mutual fit of both structures. The collections (herbarium, entomological and ornithological) kept at natural history museums would be very useful in this respect.

  .

  The clades characterized by flowers from the type Aconitum are A. subg. Aconitum, A. subg. Lycoctonum, and Gymnaconitum gymnandrum. Until 2013, the genus Aconitum consisted of three subgenera: A. subg. Aconitum; A. subg. Lycoctonum, and A. subg. Gymnaconitum. Wang et al. (2013) stated that Gymnaconitum should be an independent genus based on molecular phylogenetic data. The clades characterized by flowers from the type Delphinium are D. subg. Delphinastrum & Oligophyllon (hereafter called D. subg. Delphinastrum, for simplicity), D. subg. Anthriscifolium, D. subg. Delphinium, and Staphisagria. The genera characterized by flowers from the type Consolida are Aconitella and Consolida.

Fig

  Fig.1).

2

Acknowledgments

We thank the Muséum national d'Histoire Naturelle, Paris, for funding the internship of Yiyun Chen (ATM funding), and the stay of AN as an invited researcher at Paris Herbarium in 2016. We also thank Sophie Nadot (Editor) and one anonymous reviewer for their critical and helpful comments.

Disclosure statement

No potential conflict of interest was reported by the authors.

Notes on contributors

Yiyun Chen was an MSc student at the Muséum national d'Histoire Naturelle, Paris, and did her final year internship in 2017 under the supervision of Sylvain Gerber and Florian Jabbour.

She is currently working in China in a governmental organization for the conservation of endangered species. Jabbour and Renner (2011a, 2011b, 2012a, 2012b), [START_REF] Wang | Gymnaconitum, a new genus of Ranunculaceae endemic to the Qinghai-Tibetan Plateau[END_REF] and [START_REF] Xiang | Recircumscription of Delphinium subg. Delphinium (Ranunculaceae) and implications for its biogeography[END_REF]. The same colour code is used in Fig. 4, Supp. Fig. 1 and Supp. Table 1. The number of species sampled and the total number of species in each infratribal group is given between parentheses. Table 1. The same colour code is used in Fig. 1 and4 and Supp. Table 1. (A) PCA of the variation in the shape of the dorsal module of the flowers from the type Aconitum. (B) PCA of the variation in the shape of the dorsal module of the flowers from the type Consolida. 1. List of species and specimens (barcodes) sampled, and associated flower type. Genera and subgenera ordered as in the phylogenetic tree (bottom-up). Within a subgenus or genus, species are ordered alphabetically. We follow here the same classification of Delphinieae as in Jabbour and Renner (2012b). The same colour code is used in Fig. 1 and 4 and Supp. Fig. 1. "(H)" next to the barcode indicates that the specimen is the holotype of the species name. Digital images of the specimens can be visualized and downloaded using the following link: 1: List of species and specimens (barcodes) sampled, and associated flower type. Genera and subgenera ordered as in the phylogenetic tree (bottom-up). Within a subgenus or genus, species are ordered alphabetically. We follow here the same classification of Delphinieae as in Jabbour and Renner (2012b). "(H)" next to the barcode indicates that the specimen is the holotype of the species name. Digital images of the specimens can be visualized and downloaded using the following link:

Supplementary Table