Francesco Dell'isola

Alessandro Della Corte
email: alessan-dro.dellacorte@uniroma1.it

Leopoldo Greco

Angelo Luongo

Plane Bias Extension Test for a continuum with two inextensible families of fibers: a variational treatment with Lagrange Multipliers and a Perturbation Solution

Keywords: Inextensible fibers, Generalized continua, Variational treatment, Perturbation analysis

In the present paper, a system constituted by two families of nearly inextensible fibers is studied. The interest towards this kind of structures, which can be suitably modeled by 2D continua with inextensible curves, comes form both theoretical and applicative motivations. Indeed, the advantageous weight/strength ratio and the particularly safe behavior in fractures make them attractive in many engineering fields, while the mechanical properties displayed by them are theoretically interesting and still not well understood in full generality. We focused on dead-loading traction boundary-value problems assuming particular forms for the deformation energy associated to every node of the structure. The main result of the paper consists in the determination of the nonlinear integral equations describing fibers directions for a specific class of deformation energies. A perturbative analysis of a particular solution is also performed, and some physical justification for one of the considered forms for the energy density are provided in the Appendix.

Introduction

The conceptual framework of this paper coincides with the one set in the pioneering works by Pipkin and Rivlin, in particular in [START_REF] Ac Pipkin | Plane traction problems for inextensible networks[END_REF][START_REF] Rs Rivlin | Plane strain of a net formed by inextensible cords[END_REF]. We introduce a bidimensional continuum B whose placement is always plane. This continuum models finite and plane deformations of sheets or nets formed by two intersecting families of fibers orthogonal in the reference configuration. The intersection point of two initially orthogonal fibers is always the same material particle, which is identified with a material particle of the introduced 2D continuum. The fibers are assumed to be inextensible and the network resists shearing: this is modeled by a deformation energy which depends (only) on the shear angle between the fibers. The equilibrium configurations are characterized via the condition of stationarity of total energy. The growing interest in such systems stems both from theoretical and technological fields. From a theoretical point of view a deep understanding of the mechanics of such systems can need the further development of an exhaustive and coherent theory of generalized continua. It is well known in fact that by means of suitable homogenization procedure one can recover, from a system similar to the one considered in this paper, higher gradient continuum models (the reader interested is referred to [START_REF] Forest | Mechanics of generalized continua: construction by homogenizaton[END_REF][START_REF] Placidi | Reflection and transmission of plane waves at surfaces carrying material properties and embedded in second-gradient materials[END_REF][START_REF] Neff | Notes on strain gradient plasticity: finite strain covariant modelling and global existence in the infinitesimal rate-independent case[END_REF][START_REF] Alibert | Truss modular beams with deformation energy depending on higher displacement gradients[END_REF][START_REF] Dell'isola | At the origins and in the vanguard of peridynamics, non-local and higher-gradient continuum mechanics: An underestimated and still topical contribution of Gabrio Piola[END_REF][START_REF] Assidi | Equivalent properties of monolayer fabric from mesoscopic modelling strategies[END_REF][START_REF] Reis | Equivalent mechanical properties of auxetic lattices from discrete homogenization[END_REF][START_REF] Dell | A two-dimensional gradient-elasticity theory for woven fabrics[END_REF]). It is important to notice that higher continuum theories can be contextualized in the widest field of micromorphic continua (a general coverage of the field can be found e.g. in [START_REF] Altenbach | On generalized cosserat-type theories of plates and shells: a short review and bibliography[END_REF][START_REF] Garusi | An unsymmetric stress formulation for reissner-mindlin plates: A simple and locking-free rectangular element[END_REF][START_REF] Rinaldi | A microscale second gradient approximation of the damage parameter of quasi-brittle heterogeneous lattices[END_REF][START_REF] Federico | A transversely isotropic composite with a statistical distribution of spheroidal inclusions: a geometrical approach to overall properties[END_REF][START_REF] Neff | On material constants for micromorphic continua[END_REF][START_REF] Boutin | Microstructural effects in elastic composites[END_REF][START_REF] Federico | On the linear elasticity of porous materials[END_REF][START_REF] Neff | A unifying perspective: the relaxed linear micromorphic continuum[END_REF][START_REF] Placidi | Thermodynamics of polycrystalline materials treated by the theory of mixtures with continuous diversity[END_REF][START_REF] Alibert | Second-gradient continua as homogenized limit of pantographic microstructured plates: a rigorous proof[END_REF][START_REF] Neff | A geometrically exact micromorphic model for elastic metallic foams accounting for affine microstructure. modelling, existence of minimizers, identification of moduli and computational results[END_REF][START_REF] Mindlin | Micro-structure in linear elasticity[END_REF][START_REF] Misra | Micromechanical model for cohesive materials based upon pseudogranular structure[END_REF][START_REF] Madeo | Wave propagation in relaxed micromorphic continua: modeling metamaterials with frequency bandgaps[END_REF]).

From a technological point of view such systems are finding very wide applications in different fields (see e.g. [START_REF] Ferretti | Modeling the onset of shear boundary layers in fibrous composite reinforcements by second-gradient theory[END_REF][START_REF] Seppecher | Linear elastic trusses leading to continua with exotic mechanical interactions[END_REF][START_REF] Grillo | Evolution of a fibre-reinforced growing mixture[END_REF][START_REF] Grillo | Remodelling in statistically oriented fibrereinforced materials and biological tissues[END_REF][START_REF] Grillo | Growth, mass transfer, and remodeling in fiberreinforced, multi-constituent materials[END_REF][START_REF] Hamila | Locking in simulation of composite reinforcement deformations. analysis and treatment[END_REF][START_REF] Misra | Micromechanical model for viscoelastic-materials undergoing damage[END_REF][START_REF] Scerrato | A simple non-linear model for internal friction in modified concrete[END_REF]), including biomedical applications (see for instance [START_REF] Laurent | Computer-aided tissue engineering: Application to the case of anterior cruciate ligament repair[END_REF][START_REF] Laurent | Designing a new scaffold for anterior cruciate ligament tissue engineering[END_REF][START_REF] Laurent | Designing a new scaffold for anterior cruciate ligament tissue engineering[END_REF][START_REF] Goda | A micropolar anisotropic constitutive model of cancellous bone from discrete homogenization[END_REF][START_REF] Hollister | Porous scaffold design for tissue engineering[END_REF][START_REF] Grillo | Growth and mass transfer in multi-constituent biological materials[END_REF][START_REF] Andreaus | A 2-D continuum model of a mixture of bone tissue and bio-resorbable material for simulating mass density redistribution under load slowly variable in time[END_REF][START_REF] Giorgio | The influence of different loads on the remodeling process of a bone and bio-resorbable material mixture with voids[END_REF]).

In dead-loading traction boundary-value problems, the vectors representing the directions of the fibers of our system satisfy a pair of coupled integral equations, which are highly nonlinear.

The main result of the paper consists in the determination of such equations for a specific class of deformation energies. The model considered in the present paper is exposed and studied in detail in [START_REF] Dell'isola | Minimization of shear energy in two dimensional continua with two orthogonal families of inextensible fibers: the case of standard bias extension test[END_REF]. In the first section we provide a condensed exposition of the formalism there introduced and of the basic properties of the system, together with a suitable set of boundary conditions that define what is called a standard bias test. In the same section inextensibility and some symmetries of the system are taken into account in order to define a set of constrains and conditions that, in the framework of the principle of stationary action, will bring to the mentioned equations. To the stability of the solution of these equations are devoted sections 3 and 4, in which a perturbative analysis near a homogeneous constant solution is performed. Finally it is shown how the class of energies we consider can be suitable for describing a particular 3D-printed sample.

A variational deduction of equilibrium conditions

For the considered mechanical system the most suitable postulation process is to be framed in the Lagrangian Continuum Analytical Mechanics (see [START_REF] Dell | Variational principles are a powerful tool also for formulating field theories[END_REF][START_REF] Placidi | A variational approach for a nonlinear one-dimensional damage-elasto-plastic secondgradient continuum model[END_REF][START_REF] Placidi | A variational approach for a nonlinear 1-dimensional second gradient continuum damage model[END_REF]). The concept of balance of force itself, in the present context, becomes even more difficult to grasp than usual. 1The main difficulty to be confronted with concerns the need to take into account (in the considered mathematical model) the inextensibility of all material lines belonging to two families of parallel curves. These families are a natural choice for a reference system in the Lagrangian configuration of B. In the fundamental paper by J.Ball ( [START_REF] John | Convexity conditions and existence theorems in nonlinear elasticity[END_REF]) it is clearly stated that < <Note, however, that the constraint of inextensibility [..] is not included [in the presented treatment]. It seems possible, therefore, that solutions do not in general exist for boundary-value problems of inextensible elasticity, and that a higher order theory is required to make such constraints well behaved.> >.

We therefore attack the problem as follows: we characterize the kinematics of considered continuum by finding a suitable set (the space of configurations) which specifies all admissible states. Since it does not seem that we have available (for safely proceeding in our treatment) any general and sound mathematical result, after having characterized the kinematics of the system, we simply choose a set of deformation energies and look for the minimum in the set of placements verifying some precise boundary (kinematical) conditions.

Besides the inextensibility, another feature of the mechanical system we want to investigate can cause difficulties. The modeling of elastic systems with an intrinsically mixed continuousdiscrete nature, like the one herein considered, can indeed involve hard mathematical problems. We will try to outline the key difficulties along the path of achieving a more advanced (generaized) continuum theory in the following sections; however, the classical mathematical tools developed for the mathematical description of elasticity are still of course very useful, and in this regard useful references are [START_REF] Victor | Existence of weak solutions in elasticity[END_REF][START_REF] Federico | On the volumetric-distortional decomposition of deformation in elasticity[END_REF][START_REF] Federico | The linear elasticity tensor of incompressible materials[END_REF].

Geometry of Standard Bias Test. Placements.

Let us consider the set of material particles of the considered continuum B included in the Lagrangian space E L . We introduce an orthonormal system (O, X 1 , X 2 ) in E L to label material particles. In every configuration each material particle is placed, by a map r, into positions belonging to the Eulerian space E:

B r -→ r (B) ⊆ E.
In this paper, we assume that B is a rectangle whose sides have length l and L, described by the following conditions:

B = (X 1 , X 2 ) ∈ E L : X 1 ∈ [0, L] , X 2 ∈ - l 2 , l 2 .
We explicitly observe that the length ratio plays a relevant role. In this paper we will limit ourselves to the case of standard bias test, i.e. L = 3l.

The continuum body B is intended as a suitable homogeneous limit for a mechanical system composed by two (nearly-)inextensible families of fibers which form a uniform orthogonal net intersecting the perimeter of B with an angle of π/4. In order to take into account the material properties related to inextensibility, it is useful to introduce an adapted orthogonal reference system, (O, ξ 1 , ξ 2 ) , which is oriented according to the directions of the inextensible fibers and which uses dimensionless space coordinates. This is done by means of the following definitions:

ξ 1 := 1 l (X 1 -X 2 ) + 1 2 , ξ 2 := 1 l (X 1 + X 2 ) + 1 2 . ( 1
)
This secondly introduced reference system will be called fiber reference and the two directions ξ 1 , ξ 2 fiber directions. Clearly it will be useful to represent all deformation measures in (O, ξ 1 , ξ 2 ), since ξ 1 , ξ 2 represent the orthotropic directions of the material. We can use the coordinates (ξ 1 , ξ 2 )

to label a generic material particle of B and we will indicate with D 1 (ξ 1 , ξ 2 ) and D 2 (ξ 1 , ξ 2 ) the unit vectors tangent to the two families of fibers in the reference configuration at point (ξ 1 , ξ 2 ).

Inextensibility and consequent representation of admissible placements

Let γ ⊂ B be a rectifiable curve whose length is l: we say that γ is inextensible for an admissible placement r of B if for every arc α of γ we have that r(α) has the same length as α. To formulate this condition for every curve with support belonging to the reference configuration C * , it is necessary that the considered placement maps rectifiable curves into rectifiable curves. It seems that the space of all placements verifying such condition is the most suitable one for describing the kinematics of continua with inextensible curves. While it is quite easily seen that continuous functions of Bounded Variation (from D ⊂ R 2 taking values in R 2 ) have the aforementioned property, it is harder to prove that the inverse implication holds, and in the present work this is merely suggested as a reasonable conjecture.

The previous considerations indicate that the regularity requirements for the set of admissible placements must be consistent with inextensibility constraint. However in the present work, having in mind the results by Rivlin [START_REF] Rs Rivlin | Plane strain of a net formed by inextensible cords[END_REF], we will simply assume that the placement field r is piecewise twice continuously differentiable in B. If we indicate with F the gradient of r in the points where it is defined we get obviously that

d 1 = F • D 1 , d 2 = F • D 2 . ( 2
)
where the two vectors d 1 and d 2 represent the directions of the fibers in the current configuration.

Under Rivlin's regularity requirements the inextensibility of the fibers implies that the considered admissible placements have a rather particular form. Indeed the inextensibility constraint implies that a placement field of class C 1 in the neighborhood of a material particle P must verify the following conditions in P :

F • D 1 2 = d 1 2 = F • D 2 2 = d 2 2 = 1. (3)
Moreover, when r ∈ C 2 (∆), where ∆ is an open simply linearly connected subset of B, the inextensibility of the fibers implies that it is possible to decompose the placement field in the following manner (see [START_REF] Rs Rivlin | Plane strain of a net formed by inextensible cords[END_REF]): there exist two vector fields r (∆)

1 (ξ 1 ) and r (∆) 2 (ξ 2 ) respectively defined on the projection of ∆ on the fiber axes ξ 1 and ξ 2 such that

r (∆) (ξ 1 , ξ 2 ) = r (∆) 1 (ξ 1 ) + r (∆) 2 (ξ 2 ) . (4)
We can represent these two vector fields in the basis (D 1 , D 2 ) as follows:

               r (∆) 1 (ξ 1 ) = µ (∆) 1 (ξ 1 ) D 1 + ν (∆) 1 (ξ 1 ) D 2 r (∆) 2 (ξ 2 ) = ν (∆) 2 (ξ 2 ) D 1 + µ (∆) 2 (ξ 2 ) D 2 , (5)
being µ

(∆) 1 , µ (∆) 2 , ν (∆) 1 , ν (∆) 2
suitably regular scalar functions on the projections of ∆ on the fiber axes ξ 1 and ξ 2 .

If we indicate with I 1∆ the projection of ∆ on the axis ξ 1 and with I 2∆ the projection of ∆ on the axis ξ 2 , we have

F • D 1 2 = 1 ⇒ µ (∆) 1,1 2 + ν (∆) 1,1 2 = 1 (6)
⇒ ∃ϑ 1 (ξ 1 ) : (∀ξ

1 ∈ I 1∆ ) F • D 1 = cos (ϑ 1 (ξ 1 )) D 1 + sin (ϑ 1 (ξ 1 )) D 2 F • D 2 2 = 1 ⇒ µ (∆) 2,2 2 + ν (∆) 2,2 2 = 1 (7) ⇒ ∃ϑ 2 (ξ 2 ) : (∀ξ 2 ∈ I 2∆ ) F • D 2 = sin (ϑ 2 (ξ 2 )) D 1 + cos (ϑ 2 (ξ 2 )) D 2
Roughly speaking, the previous equations clearly means that on any subdomain ∆ whose points can be connected with segments parallel to ξ axes, the whole displacement field is known in terms of only two real functions of one real variable. Indeed once, for instance, the functions µ plus some integration constants (to be determined via suitable boundary conditions).

Imposed boundary conditions

The two segments Σ 1 and Σ 2 characterized respectively by the conditions

X 1 = 0, X 2 ∈ - l 2 , l 2 ; X 1 = L, X 2 ∈ - l 2 , l 2
are described, in the fiber reference, by the following corresponding ones:

Σ 1 := {(ξ 1 , ξ 2 ) ∈ B : ξ 1 ∈ [0, 1] , ξ 2 = 1 -ξ 1 } , (8)
Σ 2 := {(ξ 1 , ξ 2 ) ∈ B : ξ 1 ∈ [3, 4] , ξ 2 = 7 -ξ 1 } . (9)
In this paper, we impose the following boundary conditions on the two subsets Σ 1 and Σ 2 of the boundary of B:

1. vanishing displacement of the side Σ 1 , 2. imposed displacement u 0 = u 0 (D 1 +D 2 ) ||(D 1 +D 2 )|| of the side Σ 2 .
As already remarked in [START_REF] Ac Pipkin | Plane traction problems for inextensible networks[END_REF][START_REF] Rs Rivlin | Plane strain of a net formed by inextensible cords[END_REF] because of inextensibility the boundary conditions determine the placement field not only at the boundary but also in some regions inside the continuum at which they are applied. Depending on the geometry of the body and on the nature of applied boundary conditions one can distinguish different regions in considered body each of which is affected in a specific way by the applied boundary conditions. This can be done in our case by defining the following lines:

S 1 := {(ξ 1 , ξ 2 ) ∈ B : ξ 1 ∈ [0, 2] , ξ 2 = 1} , S 2 := {(ξ 1 , ξ 2 ) ∈ B : ξ 1 = 1, ξ 2 ∈ [0, 2]} , S 3 := {(ξ 1 , ξ 2 ) ∈ B : ξ 1 ∈ [1, 3] , ξ 2 = 2} , S 4 := {(ξ 1 , ξ 2 ) ∈ B : ξ 1 = 2, ξ 2 ∈ [1, 3]} (10)
S 5 := {(ξ 1 , ξ 2 ) ∈ B : ξ 1 ∈ [2, 4] , ξ 2 = 3} , S 6 := {(ξ 1 , ξ 2 ) ∈ B : ξ 1 = 3, ξ 2 ∈ [2, 4]} .
It is necessary to consider the projections of these lines on the fiber axes: we indicate with π 1 and π 2 the projection maps respectively on ξ 1 and ξ 2 and introduce the following notations:

I 10 := π 1 (Σ 1 ) = [0, 1] I 13 := π 1 (Σ 2 ) = [3, 4] (11)
I 20 := π 1 (Σ 1 ) = [0, 1] I 23 := π 1 (Σ 2 ) = [3, 4]. (12)
I 1 := π 1 (B), I 2 := π 2 (B). (13)
I 11 := [π 1 (S 2 ) , π 1 (S 4 )] = [1, 2], I 12 := [π 1 (S 4 ) , π 1 (S 6 )] = [2, 3] (14)
I 21 := [π 2 (S 1 ) , π 2 (S 3 )] = [1, 2] I 23 := [π 2 (S 3 ) , π 2 (S 5 )] = [2, 3]. (15) 9
As a consequence, a natural covering of the two intervals I 1 and I 2 , by means of intervals whose intersections reduce to single points, can be introduced (for simplicity we will call it simply a 'partition'):

I 1 = I 10 ∪ I 11 ∪ I 12 ∪ I 13 I 2 = I 20 ∪ I 21 ∪ I 22 ∪ I 23 (16)
Obviously this partition of the intervals I 1 and I 2 naturally induces a covering of B by means of the regions ∆ ij , whose intersections reduce to curves, defined as follows (also in this case we will simply refer to this covering as a 'partition'):

∆ ij := (I 1i × I 2j ) ∩ B, B = 3 i,j=1 ∆ ij . ( 17
)
The partition of B is graphically represented in Fig. 1. We can therefore assume that the placement field is twice continuously differentiable in the interior of any ∆ ij . If we indicate with r (i,j) (ξ 1 , ξ 2 ) the restriction of r to the region ∆ ij , recalling (4) we have

r (i,j) (ξ 1 , ξ 2 ) = r (i) 1 (ξ 1 ) + r (j) 2 (ξ 2 ) . ( 18
)
Because of the the global continuity condition of r, we are allowed to identify the vector fields

r (i) 1 (ξ 1 ) and r (j) 2 (ξ 2 ) along straight lines included in B.
So, if we consider two regions ∆ ij and ∆ hk with i = h or j = k, then we have respectively that

r (i) 1 (ξ 1 ) = r (h) 1 (ξ 1 ), r (j) 2 (ξ 2 ) = r (k) 2 (ξ 2 )
. This property allows us to define the continuous scalar fields µ 1 and µ 2 in the entire intervals I 1 and I 2 respectively. These two scalar functions will belong to

C 1 [1, 3].
2.4 Placement field in ∆ 00 and ∆ 33 and integral conditions implied by boundary conditions

The chosen boundary conditions on Σ 1 and Σ 2 , imply that r is the identity on ∆ 00 and a translation in ∆ 33 (for the proof the reader is referred to [START_REF] Dell'isola | Minimization of shear energy in two dimensional continua with two orthogonal families of inextensible fibers: the case of standard bias extension test[END_REF]).

In ∆ 00 we have that

r (ξ 1 , ξ 2 ) = ξ 1 D 1 + ξ 2 D 2 (19)
while in ∆ 33 we get:

r (ξ 1 , ξ 2 ) = (ξ 1 + u 01 ) D 1 + (ξ 2 + u 02 ) D 2 (20)
The continuity conditions in P i = (1, 1) and P f = (3, 3) of the placement field can be written as (see [START_REF] Dell'isola | Minimization of shear energy in two dimensional continua with two orthogonal families of inextensible fibers: the case of standard bias extension test[END_REF]):

               µ 1 (3) + ´3 1 sin (ϑ 2 (η)) dη = 3 + u 0 µ 2 (3) + ´3 1 sin (ϑ 1 (η)) dη = 3 + u 0 (21)
This condition imposes the continuity of the translation imposed to ∆ 33 with the displacement of the contiguous domains. Once one recalls that (i = 1, 2)

dµ i (ξ i ) dξ i -cos (ϑ i (ξ i )) = 0, µ i (1) = 1, (22)
it is easy to express the found integral conditions exclusively in terms of the fields ϑ i (ξ i ) as follows

               ´3 1 (cos (ϑ 1 (η)) + sin (ϑ 2 (η))) dη = 2 + u 01 ´3 1 (cos (ϑ 2 (η)) + sin (ϑ 1 (η))) dη = 2 + u 02 . (23)
This form will be useful when looking for the equilibrium shapes of considered body.

A remark can be done regarding the interface separating the undeformed triangles ∆ 00 and ∆ 33 from the other regions. In the simple model here considered, indeed, there exists an a priori determinable and geometrically constant interface between regions with different deformation. Since this fact is closely related to the particular set of boundary conditions imposed, one can easily understand that such a priori determination is not possible anymore whether more general and complex conditions are chosen. The theoretical treatment of the problems related to this fact can be rather difficult; a set of useful tools has been developed (also in connection with phase transitions) in [START_REF] Yeremeyev | The stability of the equilibrium of two-phase elastic solids[END_REF][START_REF] Pietraszkiewicz | Extended non-linear relations of elastic shells undergoing phase transitions[END_REF][START_REF] Eremeyev | Phase transitions in thermoelastic and thermoviscoelastic shells[END_REF][START_REF] Eremeyev | Thermomechanics of shells undergoing phase transition[END_REF].

Symmetric placements

In this section we will consider placement fields verifying some specific symmetry conditions.

2.5.1 Placement fields symmetric with respect to the X 1 axis.

This means that, given a point P of coordinates (ξ, η) and his symmetric P s having coordinates (η, ξ) , the following conditions hold:

               d 1 (P ) • D 1 = d 2 (P s ) • D 2 d 1 (P ) • D 2 = d 2 (P s ) • D 1 =⇒                µ 1,1 (ξ) = µ 2,2 (ξ) ν 1,1 (ξ) = ν 2,2 (ξ) (24)
Considering the symmetry of boundary conditions and the first equality in [START_REF] Madeo | Wave propagation in relaxed micromorphic continua: modeling metamaterials with frequency bandgaps[END_REF] we derive directly the following identities:

µ (ξ) := µ 1 (ξ) = µ 2 (ξ) ϑ (ξ) := ϑ 1 (ξ) = ϑ 2 (ξ) ∀ξ ∈ [0, 4] . (25)
Therefore, considering the relation between µ and ν, the kinematics of the problem is completely described by means of an unique field.

It is easily verified that a similar conclusion can be obtained replacing this type of symmetry with the invariance for a reflection with respect to the X 1 axis.

The space of configurations for 2D continua with two families of inextensible fibers in symmetric plane motion

If we assume the hypothesis of symmetry described in 2.5.1, the placement of B is completely determined by one scalar field µ (ξ) only, defined in the real interval [0, 4]. Because of the considered boundary conditions, recalling [START_REF] Placidi | Thermodynamics of polycrystalline materials treated by the theory of mixtures with continuous diversity[END_REF] and [START_REF] Alibert | Second-gradient continua as homogenized limit of pantographic microstructured plates: a rigorous proof[END_REF], µ has to be determined only in the interval I = [START_REF] Ac Pipkin | Plane traction problems for inextensible networks[END_REF][START_REF] Forest | Mechanics of generalized continua: construction by homogenizaton[END_REF].

Thanks to [START_REF] Neff | A geometrically exact micromorphic model for elastic metallic foams accounting for affine microstructure. modelling, existence of minimizers, identification of moduli and computational results[END_REF] we have also the conditions

µ (1) = 1, µ (3) + ˆ3 1 sin (ϑ (η)) dη = 3 + u 0 , (26)
or equivalently the conditions

µ (1) = 1, ˆ3 1 (cos(ϑ (η)) + sin (ϑ (η))) dη = 2 + u 0 . ( 27
)
The space of (admissible) configurations for 2D continua with two families of inextensible fibers (having everywhere orthogonal tangent vectors in the reference configuration) in plane motion is

given by the set of pairs of fields (µ, ϑ) the first of which is piecewise C 2 and the second of which is piecewise C 1 with possible discontinuities in the point ξ = 2 which, together with the conditions [START_REF] Grillo | Evolution of a fibre-reinforced growing mixture[END_REF], verify the local condition

dµ(ξ) dξ -cos (ϑ (ξ)) = 0 ∀ξ ∈ [1, 3] . (28)

Euler-Lagrange equilibrium condition for a class of factorisable energies

The equilibrium configuration will be given by means of two pairs of fields (µ i , ϑ i ) defined in the interval [START_REF] Ac Pipkin | Plane traction problems for inextensible networks[END_REF][START_REF] Forest | Mechanics of generalized continua: construction by homogenizaton[END_REF] which verify the conditions ( 23), [START_REF] Mindlin | Micro-structure in linear elasticity[END_REF], respectively two and one times continuously differentiable in the set [1, 3] -{2} , which minimize the deformation energy functional. Because of the definitions of the fields ϑ i and the definitions (2), (3) and because of the requirement of objectivity, the density of deformation energy W , in the case of first gradient continua, must be a function of the shear angle γ = ϑ 1 + ϑ 2 . As discussed in [START_REF] Dell | A two-dimensional gradient-elasticity theory for woven fabrics[END_REF] when second gradient deformation energies will be needed then the deformation energy must be assumed to depend also on ∇γ.

Some examples of deformation energy density.

We may find several physically reasonable constitutive assumptions for the deformation energy density. In the present paper we will consider the following ones

W 1 (γ) = k 2 (sinγ) 2 = k 2 (cosϑ 1 sinϑ 2 + sinϑ 1 cosϑ 2 ) 2 = k 2 (d 1 • d 2 ) 2 (29)
W 2 (γ) = N i=1 k i (ϑ 1 + ϑ 2 ) 2i (30)
These examples suggest that an interesting class to be considered is the one in which the function W can be represented as the sum of products of two functions each depending only on one of the ϑ i .

The determination of corresponding Euler-Lagrange conditions will be made easier by considering that both fields ϑ i depend only on the variable ξ i .

In the following, we will mainly focus on the energy W 1 , which our experimental work has shown as a particularly relevant case (a related article is in preparation).

Factorizable energy densities and corresponding Euler-Lagrange conditions

In this paper we will consider deformation energy densities having the following structure

W (ϑ 1 + ϑ 2 ) = M α=1 l α 1 (ϑ 1 )l α 2 (ϑ 2 ) (31)
In the previous subsection some examples of factorizable energies which are also objective were given.

The deformation energy functional is therefore given by

E DEF = ˆB M α=1 l α 1 (ϑ 1 )l α 2 (ϑ 2 ) (32)
We start by calculating its first variation relative to the variations (δµ 1 , δϑ 1 ).

To this aim we start by remarking that the region B can be represented as follows

B = ξ 1 ∈[0,4] ξ - 2 (ξ 1 ) , ξ + 2 (ξ 1 )
where the functions ξ ± 2 have as graph two curves whose union forms the perimeter of B and are defined as follows

                       ξ 1 ∈ [0, 1] ξ - 2 (ξ 1 ) := 1 -ξ 1 ξ 1 ∈ [1, 4] ξ - 2 (ξ 1 ) := ξ 1 -1 ξ 1 ∈ [0, 3] ξ + 2 (ξ 1 ) := ξ 1 + 1 ξ 1 ∈ [3, 4] ξ + 2 (ξ 1 ) := 7 -ξ 1 .
By applying Fubini's Theorem we get that the energy functional can be represented in a form which is more suitable for calculating its variation with respect to the variations (δµ 1 , δϑ 1 ). Indeed it can be transformed as follows

E DEF = ˆ4 0 dξ 1 M α=1 l 1 (ϑ 1 ) ˆξ+ 2 (ξ 1 ) ξ - 2 (ξ 1 ) l 2 (ϑ 2 (ξ 2 ))dξ 2 . ( 33
)
Further calculation will become easier if we introduce the integral operator defined as follows:

A 2 [f (ξ)] (η) := ˆξ+ 2 (η) ξ - 2 (η)
f (ξ) dξ.

Among the introduced constraints two are of local nature (see ( 22)), while two others are of integral nature (see [START_REF] Misra | Micromechanical model for cohesive materials based upon pseudogranular structure[END_REF]). Therefore the functional whose stationarity conditions will supply the searched equilibrium conditions is the following one (remark that the integration interval is now [START_REF] Ac Pipkin | Plane traction problems for inextensible networks[END_REF][START_REF] Forest | Mechanics of generalized continua: construction by homogenizaton[END_REF] as in the intervals [0, 1] and [START_REF] Forest | Mechanics of generalized continua: construction by homogenizaton[END_REF][START_REF] Placidi | Reflection and transmission of plane waves at surfaces carrying material properties and embedded in second-gradient materials[END_REF] the fields (µ 1 , ϑ 1 ) are fixed by the imposed boundary conditions)

E T OT = ´3 1 M α=1 A 2 [l α 2 (ϑ 2 )] (ξ 1 ) (l α 1 (ϑ 1 (ξ 1 ))) dξ 1 + ´3 1 λ 1 dµ 1 (ξ 1 ) dξ 1 -cos (ϑ 1 (ξ 1 )) dξ 1 + ´3 1 λ 2 dµ 2 (ξ 2 ) dξ 2 -cosϑ 2 (ξ 2 ) dξ 2 + Λ 1 µ 1 (3) + ´3 1 sin (ϑ 2 (η)) dη -3 -u 01 + Λ 2 µ 2 (3) + ´3 1 sin (ϑ 1 (η)) dη -3 -u 02 +,
where we introduced the global Lagrange multipliers Λ i (i = 1, 2) and the fields of Lagrange multipliers λ i each depending only on the variable ξ i .

By denoting with δ 1 the first variation of functional to be minimized with respect to the variations (δµ 1 , δϑ 1 ) while keeping (µ 2 , ϑ 2 ) fixed, we get

δ 1 E DEF = ´3 1 M α=1 ∂l α 1 ∂ϑ 1 A 2 [l α 2 (ϑ 2 )] + Λ 2 cos (ϑ 1 ) + λ 1 sin (ϑ 1 ) δϑ 1 dξ 1 + ´3 1 -dλ 1 (ξ 1 ) dξ 1 δµ 1 dξ 1 + (λ 1 (3) + Λ 1 ) δµ 1 (3) = 0. (34)
We can obtain similar results interchanging the role of indexes 1 and 22 . Because of the arbitrariness of the variations (δµ 1 , δϑ 1 ), we get the following Euler-Lagrange conditions:

dλ i (ξ 1 ) dξ 1 = 0; Λ i = -λ i (3) (i = 1, 2) (35)
M α=1 ∂l α 1 ∂ϑ 1 A 2 [l α 2 (ϑ 2 )] + Λ 2 cos (ϑ 1 ) -Λ 1 sin (ϑ 1 ) = 0 (36) M α=1 ∂l α 2 ∂ϑ 2 A 1 [l α 1 (ϑ 1 )] + Λ 1 cos (ϑ 2 ) -Λ 2 sin (ϑ 2 ) = 0 ( 37
)
2 This means we can apply Fubini's theorem interchanging the role of variables and express the deformation energy also as

E DEF = ´4 0 dξ2 M α=1 l2(ϑ2) ´ξ+
The Lagrange multipliers Λ i (i = 1, 2) will be determined by means of the integral conditions [START_REF] Grillo | Evolution of a fibre-reinforced growing mixture[END_REF].

Euler-Lagrange condition for the energy W 1

The factorization of the energy W 1 is easily obtained via trigonometric identities

sin 2 (ϑ 1 + ϑ 2 ) = sin 2 θ 1 cos 2 θ 2 + cos 2 θ 1 sin 2 θ 2 + √ 2 2 sin2θ 1 √ 2 2 sin2θ 2
The equation [START_REF] Goda | A micropolar anisotropic constitutive model of cancellous bone from discrete homogenization[END_REF] becomes

sin(2ϑ 1 )A 2 cos 2 (ϑ 2 ) -sin(2ϑ 1 )A 2 sin 2 (ϑ 2 ) + √ 2cos(2ϑ1)A 2 √ 2 2 sin(2ϑ 2 ) + Λ 2 cos(ϑ 1 ) -Λ 1 sin(ϑ 1 ) = 0 (38)
which when recalling the first between equalities ( 22) becomes (when assuming that both ϑ i belong to the interval 0, π 4 ),

2 dµ 1 dξ 1 1 - dµ 1 dξ 1 2 A 2 dµ 2 dξ 2 2 -2 dµ 1 dξ 1 1 - dµ 1 dξ 1 2 A 2 1 - dµ 2 dξ 2 2 + + √ 2 2 dµ 1 dξ 1 2 -1 A 2   √ 2 dµ 2 dξ 2 1 - dµ 2 dξ 2 2   + Λ 2 dµ 1 dξ 1 -Λ 1 1 - dµ 1 dξ 1 2 = 0 (39)
When the symmetry conditions are accepted, we have that Λ 1 = Λ 2 =: Λ, A 1 = A 2 =: A, and the last condition becomes the following integral relation for the field dµ dξ

2 dµ dξ 1 - dµ dξ 2 A dµ dξ 2 -2 dµ dξ 1 - dµ dξ 2 A 1 - dµ dξ 2 + + √ 2 2 dµ dξ 2 -1 A   √ 2 dµ dξ 1 - dµ dξ 2   + Λ dµ dξ -Λ 1 - dµ dξ 2 = 0 (40)
which will be studied later.

Euler-Lagrange condition for the energy

(ϑ 1 + ϑ 2 ) 2
The factorization of considered energy is easily obtained as follows

(ϑ 1 + ϑ 2 ) 2 = (ϑ 1 ) 2 + (ϑ 2 ) 2 + √ 2ϑ 1 √ 2ϑ 2
so that the 36 becomes, in the considered instance,

2ϑ 1 A 2 [1] + √ 2A 2 √ 2ϑ 2 + Λ 2 cos (ϑ 1 ) -Λ 1 sin (ϑ 1 ) = 0
which, when the symmetry conditions are accepted, become

2ϑA [1] + √ 2A √ 2ϑ + Λ (cos (ϑ) -sin (ϑ)) = 0.
The reader will easily note that this Euler-Lagrange condition excludes that the field ϑ can be constant in the interval [START_REF] Ac Pipkin | Plane traction problems for inextensible networks[END_REF][START_REF] Forest | Mechanics of generalized continua: construction by homogenizaton[END_REF].

Problem formulation for symmetric solutions

Once we recall the symmetry conditions given in 2.5.1, the problem is governed by field equations of the following form (the derivation of which is performed in [START_REF] Dell'isola | Minimization of shear energy in two dimensional continua with two orthogonal families of inextensible fibers: the case of standard bias extension test[END_REF]):

1 -2x (s) 2 2 1 -x (s) 2 a 1 (s) -x (s) b 1 (s) + Λ   1 - x (s) 1 -x (s) 2   = 0 s ∈ [1, 2] 1 -2x (s) 2 2 1 -x (s) 2 a 2 (s) -x (s) b 2 (s) + Λ   1 - x (s) 1 -x (s) 2   = 0 s ∈ [2, 3] (41)
sided by the constraint equation:

x (3) + 3 1 1 -x (s) 2 ds = 3 -u (42)
together with the boundary condition:

x (1) = 1 (43)
In Eqs [START_REF] Dell'isola | Minimization of shear energy in two dimensional continua with two orthogonal families of inextensible fibers: the case of standard bias extension test[END_REF], Λ is a Lagrangian multiplier, u the control parameter, and:

a i (s) := β i (s) α i (s) x (t) 1 -x (t) 2 dt, b i (s) := 1 - β i (s) α i (s) 1 -x (t) 2 dt (α 1 (s) , β 1 (s)) = (1, s + 1) , (α 2 , β 2 ) = (s -1, 3) (44)
Note that Eqs [START_REF] Dell'isola | Minimization of shear energy in two dimensional continua with two orthogonal families of inextensible fibers: the case of standard bias extension test[END_REF] are coupled by the integral terms. Namely, when s belong to the left (right) interval, the integrals are extended to this whole interval and a part of the right (left) interval.

Therefore, the integrals can be more conveniently broken as follows:

ˆs+1 1 (•) dt = ˆ2 1 (•) dt + ˆs+1 2 (•) dt s ∈ [1, 2] ˆ3 s-1 (•) dt = ˆ2 s-1 (•) dt + ˆ3 2 (•) dt s ∈ [2, 3] (45)
It can be shown (see the Appendix) that the problem admits solutions y (s) := x (s) which are symmetric with respect to the midpoint s = 2, i.e. y (2 -ξ) = y (2 + ξ). Exploiting the symmetry property, we can write:

I [•] := ˆs+1 1 (•) dt = ˆ2 1 (•) dt - ˆ3-s 2 (•) dt ˆ3 1 (•) ds = 2 ˆ2 1 (•) ds s ∈ [1, 2] (46)
In conclusion, the problem is recast in the form of an integral-algebraic equation in the unknowns y (s) , Λ, namely:

I y (t) 1 -y (t) 2 1 -2y (s) 2 2 1 -y (s) -1 -I 1 -y (t) 2 y (s) + Λ 1 - y (s) 1 -y (s) = 0 s ∈ [1, 2] 1 + 2 2 1 y (s) ds + 2 2 1 1 -y (s) 2 ds = 3 -u (47)
where use has been made of Eq [START_REF] Placidi | A variational approach for a nonlinear one-dimensional damage-elasto-plastic secondgradient continuum model[END_REF]. Once the problem has been solved, the original variable is found by integration:

x (s) = 1 + s 1 y (t) dt (48)
The shear strain ϑ is successively recovered using [START_REF] Mindlin | Micro-structure in linear elasticity[END_REF].

Finally, coming back to the original variable, we have:

x (s) = 1 + √ 3 2 (s -1) (53)

Perturbation solution

In this paper a numerical investigation of the solution in the neighborhood of the critical parameter u c is provided via perturbation method (an explanatory application of which can be found in [START_REF] Luongo | Perturbation methods for nonlinear autonomous discrete-time dynamical systems[END_REF]).

In the considered case the numerical problems are relatively easy. However, if one is interested in performing numerical investigations of 3D continua which model the real-world 3D-printed object graphically represented in Fig. 6, some sophisticated numerical methods would be required. In particular, a refined choice of the mesh is called for, capable to take into account all the peculiarities of the geometry and of kinematics of the object without leading to an explosion in terms of computational costs. Among them, a suitable version of isogeometric FE methods can be of use;

the reader interested can see for example [START_REF] Turco | Elasto-plastic analysis of kirchhoff plates by high simplicity finite elements[END_REF][START_REF] Cazzani | Isogeometric analysis of plane-curved beams[END_REF][START_REF] Greco | A procedure for the static analysis of cable structures following elastic catenary theory[END_REF][START_REF] Carassale | Non-linear discrete models for the stochastic analysis of cables in turbulent wind[END_REF][START_REF] Cazzani | Sensitivity analysis and optimum design of elastic-plastic structural systems[END_REF][START_REF] Descamps | Computational Design of Lightweight Structures: Form Finding and Optimization[END_REF][START_REF] Andreaus | Prediction of micromotion initiation of an implanted femur under physiological loads and constraints using the finite element method[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] Javili | A finite element framework for continua with boundary energies. part i: The two-dimensional case[END_REF][START_REF] Javili | A finite element framework for continua with boundary energies. part ii: The three-dimensional case[END_REF][START_REF] Berezovski | Interfaces in micromorphic materials: Wave transmission and reflection with numerical simulations[END_REF][START_REF] Andreaus | Modeling of Trabecular Architecture as Result of an Optimal Control Procedure[END_REF][START_REF] Andreaus | An optimal control procedure for bone adaptation under mechanical stimulus[END_REF][START_REF] Andreaus | Optimal-tuning pid control of adaptive materials for structural efficiency[END_REF].

The previous analysis revealed that no other homogeneous solutions, in addition to the critical one, are admitted by the system. However, by means of a continuity argument, it should be clear that solutions 'close enough' to the homogeneous state exist in a small neighborhood of the critical parameter u c , which we want to investigate here.

First, the control parameter is written as:

u = u c + εu 1 (54)
where u 1 the increment of u measured from the critical state, and 0 < 1 is a bookeeping (perturbation) parameter (to be removed at the end of the procedure), which remember us that the increment is small with respect to u c . Consistently, the unknowns are expanded in Taylor series of , taking the critical state as the starting point :

y (s) = y c + εy 1 (s) + ε 2 y 2 (s) + O ε 3 Λ = Λ c + εΛ 1 + ε 2 Λ 2 + O ε 3 (55)
where y i (s) , Λ i are unknown. By substituting Eqs ( 54), ( 55) into the problem equations [START_REF] Federico | On the volumetric-distortional decomposition of deformation in elasticity[END_REF] and expanding in series, the order-1 terms mutually cancel; by separately equating to zero the terms with the same power of , the following perturbation equations, at the ε-and ε 2 -order, are derived:

Order ε:

I [y 1 (s)] + 3 √ 3 + 1 √ 3 -1 -2s y 1 (s) + 2 √ 3 -2 √ 3 -1 Λ 1 = 0 2 1 y 1 (s) ds = 1 2 1 - √ 3 u 1 (56)
Order ε 2 :

I [y 2 (s)] + 3 √ 3 + 1 √ 3 -1 -2s y 2 (s) + 2 √ 3 -2 √ 3 -1 Λ 2 = = √ 3 -3 √ 3 -1 I y 2 1 (s) + 2 √ 3 -6 √ 3 -1 y 1 (s) I [y 1 (s)] + 8Λ 1 y 1 (s) - 36 + 6 √ 3 -3 s √ 3 -1 y 2 1 (s) 2 1 y 2 (s) ds = 4 1 - √ 3 2 1 y 2 1 (s) ds (57)
For any given u 1 , Eqs ( 56) and ( 57) must be solved in sequence to find first-and second-order quantities of the series (55).

First-order solution

If we try the solution y 1 = const in equations [START_REF] Greco | A procedure for the static analysis of cable structures following elastic catenary theory[END_REF], we find that this constant can only be zero, according to the fact that the unique constant solution of the nonlinear problem is y c .

Equation (56-a) can be transformed, by performing simple variable substitutions, in a Forward Differential Equation with variable coefficients. Setting:

ζ(s + 1) := ˆs+1 1 y 1 (t)dt, α(s) := - 1 3 √ 3+1 √ 3-1 -2s , β(s) := -2 √ 3 -2 √ 3 -1 Λ 1 α(s) (58)
equation ( 56) indeed becomes:

dζ(s) ds = α(s)ζ(s + 1) + β(s) (59)
For numerical evaluation we will start directly from equations ( 56), and will search a solution in the form of a truncated series of powers of the independent variable s, namely:

y 1 (s) = n k=0 A k s k (60)
where A k are unknown coefficients, and n is the degree of the approximating polynomial. By substituting the power series in Eqs (56-a,b), and separately equating to zero the coefficients of s k , with k = 0, 1, . . . , n (i.e. ignoring the residuals of higher degrees), a linear algebraic system of n + 2 equations in the unknown A 0, A 1 , . . . A n , Λ 1 is derived, having the following form:

                  
0 . . . c n-1,n-2 c n-1,n-1 c n-1,n 0 0 0 . . . 0 c n,n-1 c nn 0 b 0 b 1 . . . b n-2 b n-1 b n 0                                       A 0 A 1 . . . A n-2 A n-1 A n Λ 1                    =                    0 0 0 0 0 0 u 1                    (61)
where the b's and c's are numerical coefficients. The almost-triangular structure of the matrix of coefficients permits, by backward substitution, to get a closed-form solution for any n. By solving Eqs [START_REF] Greco | B-Spline interpolation of Kirchhoff-Love space rods[END_REF] for increasing n's, the solutions reported in Table 1 are found, where Âk := A k /u 1 and similar. There, the residual R n+1 is the value assumed by coefficient of the n + 1-th power of s, when the series truncated at the n-th term is substituted in it. It is seen that the coefficient of the series converge fast, and that the residual goes rapidly to zero. Below we give the exact expressions for the y 1 (s) and for the Lagrangian multipliers Λ 1 by truncating the series respectively at n = 0, 1, 2, 3, for the first order perturbation approximation.

Case n = 0:

y 1 (s) = - 1 + √ 3 4 u 1 (62)
For the sake of completeness it is show the convergence plots of the principal coefficients with the increasing of n We observe that, since all the Â's coefficient are negative, y (s) y c + y 1 (s) < y c ∀s when u -u c > 0(supercritical range) and y (s) y c + y 1 (s) < y c ∀s when u -u c < 0 (subcritical range).

Moreover, y (s) < 0 ∀s in the supercritical range and y (s) > 0 ∀s in the subcritical range, this enatiling a change of concavity of y (s) passing through the critical state. Second-order terms y 2 (s) are always negative, but, being oh higher order, do not change the qualitative behaviour described.

In Fig 5 the solution y(s) is plotted (up to second order) for a set of values of ε.

Appendix: Nonlinear torsion extension of a circular cylinder

In this appendix we present a simple analysis for the assumption of a particular deformation energy associated to every node of the considered model, namely deformation energy density W 2 defined in equation [START_REF] Hamila | Locking in simulation of composite reinforcement deformations. analysis and treatment[END_REF]. We assume that the structure is made of two families of parallel and equally spaced beams overlapping, with the connection between the beams of the two families being realized by means of circular cylinders (see Fig. 6). The shear deformation energy of B is then related to the torsion of the cylinders; since we need to consider large deformations, the torsion, because of Poynting effect, is in general coupled with an elongation. We consider a cylinder subjected to a combined elongation-twisting deformation (with respect to z axis). During the considered deformation the circular cross section remains undeformed. The aim of this section consists in deducing a 1D constitutive law from the 3D case by means of standard Saint Venant assumptions. Although the proposed model is satisfactory for a 3D system such that of Fig. 6 below, a generalization to micro-and nano-scale systems may require taking into account the important effects arising from a large ratio between surface and volume. Effects of this kind are becoming increasingly relevant from a theoretical and practical point of view because of the novel possibilities in micro-and nano-engineering; for useful references, the reader can see e.g. [START_REF] Va Eremeyev | Nonlinear micropolar shells: theory and applications[END_REF][START_REF] Va Eremeyev | The effective stiffness of a nanoporous rod[END_REF][START_REF] Holm Altenbach | Mechanical properties of materials considering surface effects[END_REF].

Kinematic model

We assume that the cross section remains plane during the deformation process, i.e. the warping of the cross section is neglected. The displacement consists in an arbitrary rotation, θ (z), around the centroid z and in a translation w(z) along the z axis independent of x and y coordinates on the cross section: u = x (cos θ (z) -1) -y sin θ (z) v = x sin θ (z) + y (cos θ (z) -1)

w = w(z) (73)
in which u, v are the orthogonal components in the cross section plane. The strain gradient deformation tensor associated to ( 73) is given by

F =      
cos(ϑ) -sin(ϑ) ϑ (-x sin(ϑ) -y cos(ϑ)) sin(ϑ) cos(ϑ) ϑ (x cos(ϑ) -y sin(ϑ))

0 0 w (z) + 1       (74)
The Green-Lagrange tensor, 2 E = F T F -I, associated to the displacement [START_REF] Piccardo | A complete dynamic approach to the generalized beam theory cross-section analysis including extension and shear modes[END_REF], is:

{E} ij =       0 0 -y 2 ϑ 0 0 x 2 ϑ -y 2 ϑ x 2 ϑ w + 1 2 (w ) 2 + r 2 2 (ϑ ) 2       (75)
where r 2 = x 2 + y 2 . The only non vanishing terms are:

ε z := w + 1 2 (w ) 2 + (r θ ) 2 2 , γ := γ 2 xz + γ 2 yz = rθ . ( 76
)
form W 2 is interesting and suggests that probably a more refined reduction procedure has to be considered in order to better capture the peculiarities of the mechanics of the system.

Conclusions

The considered mechanical system entails rich prospects of theoretical and application-oriented developments, as it can be generalized in many possible directions:

• More general geometry, e.g. non orthogonal or non rectilinear fibers.

• More general kinematics, e.g. non symmetrical/non uniform boundary conditions, non planar placements, non planar reference (unstressed) configuration.

• More general beam model, e.g. non exactly inextensible and shear deformable beams (on generalized beam theory the reader can see e.g. [START_REF] Bîrsan | Deformation analysis of functionally graded beams by the direct approach[END_REF][START_REF] Piccardo | A complete dynamic approach to the generalized beam theory cross-section analysis including extension and shear modes[END_REF][START_REF] Piccardo | Dynamic response of euler-bernoulli beams to resonant harmonic moving loads[END_REF][START_REF] Boutin | Generalized beams and continua. dynamics of reticulated structures[END_REF])

• More general forms for the deformation energy.

• Multiphysics systems, e.g. beams displaying piezo-flexo electric effects (for general references see [START_REF] Giorgio | Multimode vibration control using several piezoelectric transducers shunted with a multiterminal network[END_REF][START_REF] Vescovo | Dynamic problems for metamaterials: review of existing models and ideas for further research[END_REF]).

• Investigation of dynamical effects due to inertia, and in particular wave propagation. In this connection, one can observe that peculiar geometry and kinematics like the ones considered can entail the onset of various kinds of instabilities (see e.g. [START_REF] Luongo | A linear curved-beam model for the analysis of galloping in suspended cables[END_REF][START_REF] Piccardo | Some research perspectives in galloping phenomena: critical conditions and post-critical behavior[END_REF]).

Besides these possible generalizations, some delicate theoretical questions, touched in the present paper (in particular concerning the most suitable function space in which one has to settle), have to be more deeply addressed and rigorous results have to be provided. Moreover, numerical investigations, as already observed, have to be widened and results related to 2D nd 3D modeling have to be systematically compared. Finally, the development of new experimental results is also crucial

Figure 1 :

 1 Figure 1: Partition of B in ∆ ij regions.

c 1

 1 00 c 01 . . . c 0,n-2 c 0,n-1 c 0n c 0Λ c 10 c 11 . . . c 1,n-2 c

  562494 -0.567583 -0.058844 -0.009073 -0.001597 0.002796

Figure 3 :Figure 4 :

 34 Figure 3: Convergence behavior of the coefficients B0 , B1 , B2 and Λ2 for the second order perturbation approximation.

Figure 5 :

 5 Figure 5: Solution y(s)f ors ∈ [1, 3] for different values of ε in a neighborhood of the critical solution.

Figure 6 :

 6 Figure 6: 3D view of a possible physical realization of B, with a zoom showing the cylindric pivots.

  Figure 2: Convergence behavior of the coefficients Â0 , Â1 , Â2 and Λ1 for the first order perturbation approximation.

		B 0 u 2		B 1 u 2		B 2 u 2
	2.3800		0.11200
					0.00200
	2.4000		0.11400
	2.4200	A 0 u 1		A 1 u 1	0.00150	A 2 u 1
			0.11600
	2.4400				0.0090
	0.5800		0.0590 0.11800		0.00100
	2.4600
	0.6000 2.4800		0.0595 0.12000		0.0091 0.00050
	2.5000		0.0600 0.12200
	2.5200 0.6200		0.0605 0.12400		0.0092 0.00000
	0.6400 2.5400		0.0610 0.12600		0.00050
	0 2 4 6 8 10	0 2 4 6 8 10	0.0093	0 2 4 6 8 10
	0.6600	n	0.0615	n		n
	0.6800		0.0620	2 u 2	0.0094
		0 2 4 6 8 10	0 2 4 6 8 10	0 2 4 6 8 10
			48.000
		n		n		n
				1 u 1
			50.000
			6.6000
			6.8000 52.000
			7.0000 54.000
			7.2000
			56.000
			7.4000
			0 2 4 6 8 10
			7.6000	n
			7.8000
			0 2 4 6 8 10
				n

As a brief historical remark, we simply want to recall that in the Traité de Dynamique by D'Alembert one finds the following (very impressive) statement:I have proscribed completely the forces relative to the bodies in motion, entities obscure and metaphysical, which are capable only to throw darkness on a Science which is clear by itself.

(ξ

) ξ - 1 (ξ 2 )l1(ϑ1 (ξ1))dξ1 .

Critical solution

A simple inspection of Eq (47-a) reveals that y = y c = const are, in principle, solutions, representing homogeneous states of strain for the system. However, since the problem is nonlinear, the number, values and character (real or complex) of the the constant solutions are unknown. We denote by a subscript c all the quantities associated with the constant solution, namely (y c , Λ c , u c ), and refer to it as the 'critical solution', since it possibly occurs at critical values of the control parameter.

By substituting the ansatz in Eq (47-a), and performing the integrations, we find:

By requiring that it is satisfied for any s, two algebraic equations are derived:

They admit the unique non-trivial solution in the real field:

after having disregarded the opposite one for physical reasons.

With y = y c now known, Eq (47-b)) provides the critical value of the control parameter:

consistently with the numerical results.

Case n = 1:

Case n = 2:

Case n = 3:

Second-order solution

In the equations ( 57) the first-order quantities y 1 (s) and Λ 1 are known. They can be transormed in FDE in a similar way as done before, and here are solved again in power series, by taking:

where B k are unknown coefficients. By proceeding as before, a linear algebraic system in the unknown is obtained, in which the matrix is the same than Eq [START_REF] Greco | B-Spline interpolation of Kirchhoff-Love space rods[END_REF], but the known terms are much more involved. By summarizing, the second-order solution (with ε reabsorbed) reads:

The original variable, consequently is:

Elastic energy and constitutive equations

By assuming a linear elastic material, the potential energy for the considered model reduces to

where A is the cross section area. By substituting [START_REF] Giorgio | Multimode vibration control using several piezoelectric transducers shunted with a multiterminal network[END_REF], and integrating with respect to A, we obtain:

where:

Let w = ε and θ = κ. The density of the energy for the 1D model becomes then:

so that the constitutive equation are given by:

in which N and M are the the axial stress and the twisting moment respectively. The non linear coupled constitutive equations for the reduced 1D theory are given by:

(82)

Condensation of the constitutive equations

If there is no axial stress, we can assume N (z) = 0 to define the ε = ε(κ) and finally one has M = M (κ), i.e. exists a local axial deformation on the generic fiber that globally produces N = 0.

But the condition, N = 0, in not solvable in an easy analytical way, so we assume for the ε the subsequent even power series expansion:

in this way the elongation id independent on the sign of the torsion. Substituting in the equation of N = 0 by means of the polynomial identity principle we get a linear set of equations defining the coefficients c i . In this way, considering only the first three terms of the series (83), we have:

so that, the twisting moment is M (ε, κ) = M (κ):

in which only odd powers of κ are present, and then the sign of M (κ) changes according to the sign of the κ deformation component. It is easily verified that the related deformation energy is of the polynomial form W 2 above considered. As already mentioned, experimental results are indicating the alternative form W 1 as a reasonable assumption. The fact that basic assumptions lead to the in this line of investigation. The joint effort of researchers from many areas is therefore called for to achieve substantial progresses on this topic, as is always the case with rich and potentially fruitful subjects.