Francesco Dell'isola
email: francesco.dellisola@uniroma1.it

Massimo Cuomo
email: mcuomo@dica.unict.it

Leopoldo Greco
email: leopoldo.greco@virgilio.it

Alessandro Della Corte
email: alessandro.dellacorte@uniroma1.it

Bias extension test for pantographic sheets: numerical simulations based on second gradient shear energies

Keywords: Bias extension test, Constrained optimization, Inextensible fibre network, Pantographic sheet, Second gradient theory

We consider a bi-dimensional sheet consisting of two orthogonal families of inextensible fibres. Using the representation due to Rivlin and Pipkin for admissible placements, i.e. placements preserving the lengths of the inextensible fibres, we numerically simulate a standard bias extension test on the sheet, solving a non-linear constrained optimization problem. Several first and second gradient deformation energy models are considered, depending on the shear angle between the fibres and on its gradient, and the results obtained are compared. The proposed numerical simulations will be helpful in designing a systematic experimental campaign aimed at characterizing the internal energy for physical realizations of the ideal pantographic structure presented in this paper.

Introduction

The objective of this work is the analysis of a rectangular sheet composed of two orthogonal families of inextensible fibres with low bending stiffness. The two families intersect the sides of the rectangle at an angle of ±45 • in the initial configuration; we refer to this model as a 2D-continuum pantographic model. Mechanical systems of this kind are of great practical relevance in engineering applications. Indeed, the forming procedure of woven fabric represents the initial step in the production of many composite materials in both the aerospace and automotive industries. Hence, simulations dealing with the formation of woven fabrics (impregnated with polymer or dry) are very useful for pattern designing in the aforementioned context (e.g. [START_REF] Nikopour | Concentrated loading of a fibre-reinforced composite plate: experimental and computational modeling of boundary fixity[END_REF][START_REF] Selvadurai | Transverse elasticity of a unidirectionally reinforced composite with an irregular fibre arrangement: Experiments, theory and computations[END_REF][START_REF] Hamila | Locking in simulation of composite reinforcement deformations. Analysis and treatment[END_REF][START_REF] Boisse | Composite reinforcements for optimum performance[END_REF][START_REF] Nikopour | Torsion of a layered composite strip[END_REF]).

A bias extension test, in which a displacement at 45 • with the fibres is imposed, is widely used to characterize the mechanical behaviour of a woven fabric under large shear deformations, the curve force-displacement (along the bias direction) presenting the typical hardening effect due to geometric effects. Most commonly used fibres have a low (but non-zero) bending stiffness. In this work we analyse the bias extension test considering both first and second gradient deformation energy densities, continuing in the line of investigation introduced 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][START_REF] Dell'isola | A two-dimensional gradient-elasticity theory for woven fabrics[END_REF][START_REF] Scerrato | Three-dimensional instabilities of pantographic sheets with parabolic lattices: numerical investigations[END_REF][START_REF] Scerrato | On the effect of shear stiffness on the plane deformation of linear second gradient pantographic sheets[END_REF] or, with the same topology but considering extensible beams instead of inextensible fibres, in [START_REF] Dell'isola | Elastic pantographic 2D lattices: a numerical analysis on the static response and wave propagation[END_REF][START_REF] Dell'isola | Large deformations of planar extensible beams and pantographic lattices: Heuristic homogenization, experimental and numerical examples of equilibrium[END_REF].

In this work, a new numerical scheme is implemented. Two main novelties characterize the present work. The first consists of a change in variables which reduces the order of differentiation of the problem formulation, making it a first gradient (rather than a second gradient) problem. Specifically, we adopt the rotations of the two fibres as degrees of freedom and assume for each rotation field a particular structure conforming with the Pipkin-Rivlin representation. The inextensibility constraint is therefore automatically satisfied. The displacement is successively evaluated by means of numerical integration of the two rotation fields. This new formulation of the problem leads clearly to an improved convergence as the computational cost is reduced with respect to a displacement-based formulation with Lagrangian multipliers accounting for the inextensibility constraint. The second novelty is the simulation of new energy models involving the shear angle and its gradient in the form g 2,1 (see subsequent section on numerical results).

The shape ratio of the specimen is 3:1, i.e. the short clamped side is one-third of the long side. The standard boundary conditions of the bias extension test are considered, i.e. the long sides of the rectangle are free, one of the short sides is fixed and on the other a prescribed displacement, parallel to the long sides, is imposed.

The paper is structured as follows. In the second section we give a brief description of the model, which is basically the one set 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][START_REF] Dell'isola | Plane Bias Extension Test for a continuum with two inextensible families of fibers: a variational treatment with Lagrange Multipliers and a Perturbation Solution[END_REF], showing that by means of geometry and symmetry considerations, the kinematics of the system can be described in terms of one scalar field (representing the value of the shear deformation angle). In the third section a numerical formulation of the problem is developed, and a certain set of possible energies is given, with some heuristic considerations on our choices. In the fourth section the numerical results are presented and discussed.

Theoretical and technological contextualization of proposed line of investigation

The growing importance of the type of systems considered here arises both from their high strength-to-weight ratio and from their peculiarly safe behaviour under fracture. In fact, some theoretical and experimental results show that fracture can be forecast rather accurately, and it is also worth noting that from the start of the fracture to the complete split of the sample, the system is still able to store a considerable amount of deformation energy. Moreover, the development of 3D printing and other computer-aided manufacturing techniques allows for a physical realization of the most exotic geometries (see [13] for recent applications of advanced manufacturing). Therefore, a thorough understanding of the mechanical behaviour of the simple geometry proposed is a necessary first step towards the development of increasingly complex mechanical systems. The considered structure, because of its characteristics, can be indeed considered as a basic case of metamaterial, which can be enriched in different ways by generalizing its geometry, its mechanical properties or its kinematics (for a recent review on metamaterials, see [START_REF] Vescovo | Dynamic problems for metamaterials: review of existing models and ideas for further research[END_REF]).

The theoretical framework in which it is natural, in the opinion of the authors, to study the behaviour of the considered system (and in general that of systems displaying complex microstructures) is the variational approach (see [START_REF] Dell'isola | Variational principles are a powerful tool also for formulating field theories[END_REF] for a general discussion), which of course is also the most suitable for the numerical formulation of the problem. The theoretical tools used for analysing the test are those set out in the works of Rivlin and Pipkin (see [START_REF] Rivlin | Plane strain of a net formed by inextensible cords[END_REF][START_REF] Rivlin | Plane strain of a net formed by inextensible cords[END_REF] and [START_REF] Pipkin | Some developments in the theory of inextensible networks[END_REF][START_REF] Pipkin | Plane traction problems for inextensible networks[END_REF][START_REF] Wang | Inextensible networks with bending stiffness[END_REF][START_REF] Wang | Plane deformations of nets with bending stiffness[END_REF]). On the basis of the principle of stationary action, some suitable forms of energy, depending only on shear deformations, have been postulated, considering both first and second gradient energy models. The problem lies in the field of non-linear constrained optimization in which the objective function is represented by the energy functional subjected to impenetrability and kinematical compatibility constraints. We formulate the problem directly considering Pipkin's representation.

It is clear nowadays that a suitable homogenization (see [START_REF] Reis | Equivalent mechanical properties of auxetic lattices from discrete homogenization[END_REF][START_REF] Goda | A micropolar anisotropic constitutive model of cancellous bone from discrete homogenization[END_REF][START_REF] Alibert | Second-gradient continua as homogenized limit of pantographic microstructured plates: a rigorous proof[END_REF][START_REF] Forest | Mechanics of generalized continua: construction by homogenizaton[END_REF][START_REF] Reis | Construction of micropolar models from lattice homogenization[END_REF][START_REF] Assidi | Equivalent properties of monolayer fabric from mesoscopic modelling strategies[END_REF][START_REF] Goda | Construction of micropolar models from lattice homogenization[END_REF][START_REF] Reis | Homogenized elastoplastic response of repetitive 2D lattice truss materials[END_REF][START_REF] Chaouachi | A micromechanical model of woven structures accounting for yarn-yarn contact based on Hertz theory and energy minimization[END_REF] for some relevant examples) of such systems leads to models that can be properly contextualized in the field of micromorphic continua (see [START_REF] Misra | Based micromorphic model predicts frequency band gaps[END_REF][START_REF] Misra | Identification of higher-order elastic constants for grain assemblies upon granular micromechanics[END_REF]). A reader interested in the theory of micromorphic continua can see the classical references [START_REF] Mindlin | Micro-structure in linear elasticity[END_REF][START_REF] Toupin | Theories of elasticity with couple-stress[END_REF][START_REF] Eringen | Microcontinuum field theories: I. Foundations and solids[END_REF][START_REF] Eringen | Theory of micropolar fluids[END_REF] and more recent ones [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] Neff | A unifying perspective: the relaxed linear micromorphic continuum[END_REF][START_REF] Altenbach | Acceleration waves and ellipticity in thermoelastic micropolar media[END_REF][START_REF] Altenbach | On generalized Cosserat-type theories of plates and shells: a short review and bibliography[END_REF][START_REF] Eremeyev | Nonlinear micropolar shells: theory and applications[END_REF][START_REF] Boutin | Microstructural effects in elastic composites[END_REF], while for applications the following references are relevant: [START_REF] Berezovski | Interfaces in micromorphic materials: wave transmission and reflection with numerical simulations[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][START_REF] Madeo | Wave propagation in relaxed micromorphic continua: modeling metamaterials with frequency band-gaps[END_REF][START_REF] Federico | On the linear elasticity of porous materials[END_REF][START_REF] Misra | Micromechanical model for cohesive materials based upon pseudo-granular structure[END_REF][START_REF] Misra | Micromechanical model for viscoelastic-materials undergoing damage[END_REF][START_REF] Andreaus | Modeling of the interaction between bone tissue and resorbable biomaterial as linear elastic materials with voids[END_REF][START_REF] Scerrato | A simple non-linear model for internal friction in modified concrete[END_REF][START_REF] Scerrato | A micro-structural model for dissipation phenomena in the concrete[END_REF].

We propose in this paper some types of energy depending on the second gradient of displacement, embedding our model in the mathematical framework of generalized continua, and in particular of higher gradient continua theory, which by now relies on a solid theoretical and experimental background. For more on this topic, the interested reader may consult [START_REF] Germain | The method of virtual power in continuum mechanics. Part 2: Microstructure[END_REF][START_REF] Mindlin | Second gradient of strain and surface-tension in linear elasticity[END_REF] as standard references and [START_REF] Alibert | Truss modular beams with deformation energy depending on higher displacement gradients[END_REF][START_REF] Yang | Higher-order continuum theory applied to fracture simulation of nanoscale intergranular glassy film[END_REF][START_REF] Yang | Higher-order stress-strain theory for damage modeling implemented in an element-free Galerkin formulation[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] Rinaldi | A microscale second gradient approximation of the damage parameter of quasi-brittle heterogeneous lattices[END_REF][START_REF] Placidi | A variational approach for a nonlinear 1-dimensional second gradient continuum damage model[END_REF][START_REF] Placidi | A variational approach for a nonlinear one-dimensional damage-elasto-plastic second-gradient continuum model[END_REF][START_REF] Placidi | Reflection and transmission of plane waves at surfaces carrying material properties and embedded in second-gradient materials[END_REF] for more recent results in the field.

Theoretical model and set-up for bias extension test

In this section we consider a continuous model at every point of which two fibres (belonging to two orthogonal families) intersect. The fibres are assumed to be inextensible. It should be remarked that, in general, difficult theoretical and numerical problems can arise from constraints of this kind (for related problems, see for instance [START_REF] Andreaus | Soft-impact dynamics of deformable bodies[END_REF][START_REF] Selvadurai | On the surface displacement of an isotropic elastic halfspace containing an inextensible membrane reinforcement[END_REF][START_REF] Federico | The linear elasticity tensor of incompressible materials[END_REF]). As we will see, indeed, the inextensibility condition is naturally expressed by means of integral terms, which are in general more difficult to handle than differential ones. The motion is described in a set of coordinates ξ 1 , ξ 2 which have the same directions as fibres.

Rivlin-Pipkin decomposition theorem

Let D 1 and D 2 be the unit tangent vectors along the fibres in the reference configuration (Fig. 2a), and let

d 1 = R 1 (θ 1 )D 1 , d 2 = R 2 (θ 2 )D 2 (1)
be the in-plane rotated vectors by means of two independent rotation operators indicated by R 1 (θ 1 ) and R 2 (θ 2 ). Since we consider only plane motion, the axial vector of the rotation operators is simply characterized by the rotation angles θ 1 , θ 2 . Notice that the rotations θ 1 , θ 2 are considered positive if the angle between the directions D 1 , D 2 decreases. Because of the inextensibility constraints, the deformation tensor is then given by

F = (R 1 D 1 ) ⊗ D 1 + (R 2 D 2 ) ⊗ D 2 . ( 2
)
Rivlin (in [START_REF] Rivlin | Plane strain of a net formed by inextensible cords[END_REF]) proved that, in our assumptions, the placement splits into the sum of two vector functions as

r (ξ 1 , ξ 2 ) = r 1 (ξ 1 ) + r 2 (ξ 2 ). ( 3
)
Therefore, we have that

θ 1 = θ 1 (ξ 1 ) and θ 2 = θ 2 (ξ 2 ).
The displacement field can be expressed in integral form observing that

d i = D i + ∂u ∂ξ i (4)
and, using Eqs. ( 1) and ( 3), we obtain by integration the relative displacement between two points P, Q of the specimen:

u Q x = ξ 1 Q ξ 1 P (cos(θ 1 ) -1) dξ 1 + ξ 2 Q ξ 2 P sin(θ 2 )dξ 2 D 1 • e x + ξ 1 Q ξ 1 P sin(θ 1 )dξ 1 + ξ 2 Q ξ 2 P (cos(θ 2 ) -1) dξ 2 D 2 • e x + u P x , u Q y = ξ 1 Q ξ 1 P (cos(θ 1 ) -1) dξ 1 + ξ 2 Q ξ 2 P sin(θ 2 )dξ 2 D 1 • e y + ξ 1 Q ξ 1 P sin(θ 1 )dξ 1 + ξ 2 Q ξ 2 P ) cos(θ 2 ) -1)dξ 2 D 2 • e y + u P y , (5)
where e x , e y are the unit vectors of the external reference frame.

Energy functional

We proceed now to the identification of a suitable energy functional for the considered structure.

Let γ = d 1 • d 2 = sin(ϑ 1 + ϑ 2 )
be the shear between the fibre directions. In order for the transformation to be monodrome (single-valued), the shear must belong to the interval

θ 1 + θ 2 ∈ - π 2 , π 2 ⇒ γ ∈ (-1, 1). ( 6
)
The Cauchy deformation tensor is then

C = D 1 ⊗ D 1 + γ (D 1 ⊗ D 2 + D 2 ⊗ D 1 ) + D 2 ⊗ D 2 . ( 7
)
The Green's strain tensor is

E = γ 2 (D 1 ⊗ D 2 + D 2 ⊗ D 1 ) . ( 8
)
Here we consider also the contribution of the second gradient of the displacements to the deformation. Since the basis vectors D i are constant, the only non-zero terms of the gradient of the strain tensor E are

∇ E 121 = ∇ E 211 = 1 2 cos(θ 1 + θ 2 ) ∂θ 1 ∂ξ 1 , ∇ E 122 = ∇ E 212 = 1 2 cos(θ 1 + θ 2 ) ∂θ 2 ∂ξ 2 . ( 9
)
The elements of the gradient of the strain tensor can be given a very useful interpretation. Indeed, it can be easily shown that ∂θ i ∂ξ i is the in-plane curvature of i-family fibres. The elements of ∇E are the projection of these curvatures on the deformed directors and can be interpreted as

∇ E 121 = d 1 × d 2 • e z ∂θ 1 ∂ξ 1 , ( 10
)
that is, the curvatures multiplied by the area of the deformed element, so that the second gradient deformation vanishes as the shear tends to 1, i.e. when the fibres tend to overlap. This interpretation will be useful for the choice of the second gradient energy. Any isotropic first gradient energy functional must then be a function of γ 2 :

F 1 = B g 1 γ 2 dB. ( 11
)
We also include in the energy functional a contribution of the strain gradient, given by the norm of ∇E,

∇E 2 = 1 2 ∇γ 2 , ( 12
)
so that, also accounting for the impenetrability constraints (6), the energy functional F becomes

F = B αg 1 γ 2 + βg 2 ∇γ 2 dB + U K (γ ), ( 13
)
in which we have introduced two constitutive parameters α and β, and U K is the indicator function of the set

K = (θ 1 , θ 2 ) : γ ξ 1 , ξ 2 ∈ (-1, 1), ∀(ξ 1 , ξ 2 ) , U K = 0, (θ 1 , θ 2 ) ∈ K , +∞, otherwise. (14)
The latter condition implies a constraint on the optimization problem that can be included in various ways (e.g. penalty, Lagrangian multipliers).

First and second gradient energy models

For a first gradient deformation problem, a sufficient condition for the existence of a solution of the BVP (Boundary value problem) is that the energy functional is a convex function of the strain tensor. Therefore, only convex functionals F 1 are considered in this work. The following four forms of g 1 (γ ) are used in the numerical analyses:

-First gradient shear energy model:

g 1,1 (γ ) = 1 2 γ 2 = 1 2 sin 2 (θ 1 + θ 2 ); (15)
-First gradient quadratic energy model:

g 1,2 (γ ) = 1 2 arcsin 2 (γ ) = 1 2 (θ 1 + θ 2 ) 2 ; (16)
-First gradient quartic energy model:

g 1,3 (γ ) = 1 4 arcsin 4 (γ ) = 1 4 (θ 1 + θ 2 ) 4 ; (17)
-First gradient quadratic and quartic energy model, [α = 1 in (13)]:

g 1,4 (γ ) = α 1 arcsin 2 (γ ) + α 2 arcsin 4 (γ ) = α 1 (θ 1 + θ 2 ) 2 + α 2 (θ 1 + θ 2 ) 4 . ( 18
)
Figure 1a shows a sketch of the four energy models considered over the existence interval for the variable γ . Note that the functions g 1,2 and g 1,4 can be considered as series expansions of g 1,1 , so that for small shears they coincide. However, the function g 1,3 is quite different, as is evidenced from the plot of the conjugated internal stress reproduced in Fig. 1b, showing that the stress remains very small for small values of the shear, then rapidly grows. This function can simulate the effect of the densification of the fibres as the shear grows. Two energy functionals accounting for the second gradient deformation are considered; the first, as stated, is the L 2 norm of the second gradient strain tensor, and the other is a truncated expansion of it:

g 2,1 (∇γ ) = ∇γ 2 = 1 2 cos 2 (θ 1 + θ 2 ) ∂θ 1 ∂ξ 1 2 + ∂θ 2 ∂ξ 2 2 ; (19)
-Second gradient quadratic energy model:

g 2,2 (γ ) = 1 2 ∂θ 1 ∂ξ 1 2 + ∂θ 2 ∂ξ 2 2 . ( 20
)
The existence and uniqueness of the minimum for the given energy functional is in general a non-trivial issue and is linked to the properties of the convexity, coercivity and lower semi-continuity of the functional. We conjecture that in our case there are no local minima in symmetric extensions, but this point, as well as a rigorous mathematical characterization of the proposed models, deserve future investigations in specific works. The dual stresses, for the functional (13), are given by

σ 1 j = α ∂ θ 1 (g 1, j ), σ 2 j = α ∂ θ 2 (g 1, j ), j = 1, . . . , 4, τ 1 i = β ∂ θ 1 (g 2,i ), τ 2 i = β ∂ θ 2 (g 2,i ), i = 1, 2, (21)
where the definitions of the first gradient energies, g 1, j , are given in ( 15), ( 16), ( 17) and ( 18), while the second gradient energies, g 2,i , are given in [START_REF] Pipkin | Plane traction problems for inextensible networks[END_REF] and [START_REF] Wang | Inextensible networks with bending stiffness[END_REF]. In what follows, we limit ourselves to boundary conditions involving free normal derivatives of the kinematical descriptors and, therefore, null double forces.

Geometry of the problem

A quadrilateral specimen with a ratio of 1:3 of the sides subjected to a bias extension test is considered, as indicated in the introduction. Figure 2 shows the geometry of the sample and the intrinsic coordinate system considered. The x-axis coincides with the longitudinal direction of the specimen. Under our assumptions the displacement field of the specimen can be divided into seven subregions, Δ 00 , Δ 01 , Δ 10 , Δ 23 , Δ 32 , Δ 33 and Δ = Δ 11 ∪ Δ 12 ∪ Δ 21 ∪ Δ 22 (Figs. 2a and3c), in each of which the displacement field obtained by integration of the strain field attains a specific functional form. In regions Δ 00 and Δ 33 the deformation vanishes. This is experimentally verified by the simple extension test on a 3D-printed polymeric specimen (Fig. 3).

A limit condition exists corresponding to the situation where, in the central region, the relative rotation between the fibres becomes π/2. No further deformation is allowed beyond this point since the fibres overlap and cannot stretch. The displacement that corresponds to this limit condition can be evaluated through the first equation, Eq. ( 5), setting θ 1 = θ 2 = π/4, and is found to be

u lim = 2 √ 2( √ 2 -1) = 0.82 √ 2.

Boundary conditions

Indicating by u the imposed displacement along the free short side of the specimen, and by u x = u • e x and u y = u • e y its components, the bias extension test corresponds to the case u x = u, u y = 0. Using the first of expressions (5) the boundary condition becomes (Fig. 2a)

3 1 cos(θ 1 (ξ 1 ))dξ 1 + 3 1 sin(θ 2 (ξ 2 ))dξ 2 -2 = u √ 2 , 3 1 cos(θ 2 (ξ 2 ))dξ 2 + 3 1 sin(θ 1 (ξ 1 ))dξ 1 -2 = u √ 2 , ( 22
)
which represents two line integrals along the line defined by the two points {1, 1} and {3, 3} (the internal vertexes of the two rigid triangles Δ 00 and Δ 33 ). Concerning higher-order boundary conditions, on the short side on which the displacement is imposed, no condition is imposed on the normal derivative of the field θ 1 and θ 2 , and therefore, according to the weak formulation of the problem, the corresponding double forces (i.e. the dual quantities of the normal derivative in the energy) are zero. For the same reason, on the free sides of the specimen, both the force and double force are assumed to be zero.

Numerical formulation

In this work we approach the problem as a constrained minimization in which the objective functional is expressed in terms of two discretized rotation fields, ϑ 1 (ξ 1 ) and ϑ 2 (ξ 2 ), characterizing the unit tangent vectors to the fibres. A quadrilateral mesh consistent with the ratio 1:3 of the geometry of the specimen is introduced. We consider a piecewise constant (P 0 -)interpolation for the two rotation fields ϑ 1 (ξ 1 ) and ϑ 2 (ξ 2 ) on the uniform quadrilateral mesh. The discretized values of the rotations are denoted by ϑ 1,i j and ϑ 2,i j , with i = 1, . . . , n e -1 and j = 1, . . . , 3n e -1 and n e the number of elements along the short side of the specimen. According to Rivlin's representation, all the elements having the same value ξ 1 (resp. ξ 2 ) of the centroid coordinate have the same value for ϑ 2 (resp. ϑ 1 ), as represented in Fig. 3, where the element with zero rotation according to the boundary conditions is also shown.

a θ 1,i j , b θ 2,i j , c Θ i j = θ 1,i j + θ 2,i j

Gradients of rotations

Since a piece-wise constant interpolation (or P 0 -interpolation) has been used, the gradients of the rotations cannot be directly evaluated. They have therefore been collocated in the vertices of the mesh, and we evaluate for each node i, j the gradients of the two scalar rotation fields as

∂ S 1 ϑ 1 i, j = ϑ 1,i j+1 -ϑ 1,i+1 j √ 2h e , ∂ S 2 ϑ 2 i, j = ϑ 2,i+1 j+1 -ϑ 2,i j √ 2h e , ( 23
)
with h e = L/n e the element side length, while the second gradient deformation was calculated as

∇ E 121 = cos ϑ m 1,i j + ϑ m 2,i j ∂ S 1 ϑ 1 i, j , (24)
in which ϑ m 1,i j and ϑ m 2,i j are mean values evaluated by averaging the rotations in the four elements connected at the selected node, and S 1 is the arclength corresponding to the parametric abscissa ξ 1 .

Constrained minimization problem

The discrete problem reduces to a constrained optimization in which the energy functional, which depends only on the angle of rotation fields, represents the objective function, with inequality constraints on the angle fields (the incompressibility conditions) and equality constraints for the imposition of the displacement at the free edge in terms of the rotation:

min (ϑ 1 ,ϑ 2 ) F ϑ 1 , ϑ 2 , ∂ϑ 1 ∂ξ 1 , ∂ϑ 2 ∂ξ 2 subject to : ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ -1 < ϑ 1,i j + ϑ 2,i j < 1, ( Inequality Constraint) l e d e h=1 (cos(ϑ 1h ) + sin(ϑ 2h )) -2 = u/ √ 2, (Equality Constraint) l e d e h=1 (cos(ϑ 2h ) + sin(ϑ 1h )) -2 = u/ √ 2, (Equality Constraint) (25)
A Newton method with an interior point constraint was used to solve the constrained minimization problem. The objective function is given by one of the expressions ( 15), ( 16), ( 17), [START_REF] Pipkin | Some developments in the theory of inextensible networks[END_REF] for the first gradient energy models and [START_REF] Pipkin | Plane traction problems for inextensible networks[END_REF] or [START_REF] Wang | Inextensible networks with bending stiffness[END_REF] for the second gradient energy models previously considered.

Numerical results

In this section we present and compare the results obtained with the energy models considered. The distribution of the rotation angles as well as the deformation of the specimen will be analysed. The deformation of the specimen undergoing a bias extension test reproduces the regions indicated in Fig. 2a. In the central part of the specimen, the rotations show a continuous variation, either increasing or decreasing, going from the centre of the specimen towards its fixed part. It will be shown that this variation characterizes the influence of the second gradient deformation on the solution; therefore, an indicator parameter is introduced (Fig. 4). With reference to Fig. 5, we measure the variation in the angle between the fibre directions at points p and c located at the vertex of the fixed region Δ 00 and at the centre of the specimen,

Δϑ p = π 2 -(ϑ 1 (p) + ϑ 2 (p)) and Δϑ c = π 2 -(ϑ 1 (c) + ϑ 2 (c)
), and we define as the indicator parameter the difference Δϑ pc :=Δϑ c -Δϑ p . It will be shown that in the case of a first gradient energy model the angle Δϑ pc turns out positive, that is the total rotation at the centre of the specimen is smaller than at the vertex, and the opposite occurs for second gradient dominated solutions.

First gradient energy models

To verify the ability of the numerical model to correctly represent the solution of the problem, a convergence analysis for the relative energy error was performed (considering an overkill solution as an exact solution), and its results are presented in Fig. 6, which refers to a value of the displacement u 1 = u 2 = 0.6 for three different energy cases specified in Sect. 2.3 (shear, quadratic and quartic energy). The square of the energy error is reported against the mesh size h in a log-log plot. In any case, the rate of convergence reaches the expected rate of 2, with the same accuracy for all considered cases.

Shear energy model

First we analyse the results obtained using the quadratic energy in γ [START_REF] Dell'isola | Variational principles are a powerful tool also for formulating field theories[END_REF]. Plots of the relative rotation between the fibres ϑ 1 +ϑ 2 for three different values of the end displacement are shown in Fig. 7. For any value of the displacement the total rotation exhibits sharp jumps corresponding to the boundaries of the regions represented in Fig. 2. In the two triangles next to the ends of the specimens, the rotations are zero; then two sets of discontinuity lines are found. The two fields ϑ 1 and ϑ 2 and their sum are also represented as contour plots in Fig. 8, in which the same scale is used for the three values considered for the displacement, i.e. u 1 = u 2 = 0.6, √ 3 -1 and 0.8 respectively. As appears from the plots, the rotation is not constant in the central region, so that the parameter Δϑ pc indicating the difference in the rotation between points p and c (Fig. 5) is non-zero. It is interesting to examine more closely the results for the relative rotation of the central region of the specimen (given by Δ 11 ∪ Δ 12 ∪ Δ 21 ∪ Δ 22 ; see Fig. 2a) shown in Fig. 7. From Fig. 7a shows that for small values of the imposed displacement (u y < u c = √ 3 -1), the relative rotation at the centre of the specimen (point C; see Fig. 5) is lower than the relative rotation at the boundary (point P; see Fig. 5). The contour plot of the relative angles is shown in Fig. 8a.

However, there exists a critical value of the end displacement component,

u 1 = u 2 = u c = √ 3 -1,
for which both rotation fields become piece-wise constant. This value corresponds to a change in the sign of Δϑ pc (Fig. 7). For this value of the end displacement we find Δϑ pc = 0 (Fig. 7b).

For large values of the imposed displacement the relative rotation at the centre becomes greater than the rotation at the boundary of the central region of the specimen. The change in the sign of the difference of these relative rotations is a characteristic trend only for the first gradient shear energy model.

The relative rotations Δϑ p at point p, the angle Δϑ c at point c and their difference Δϑ pc are plotted in Fig. 9 as a function of the imposed end displacement. It can be observed that for values of the displacement u < u c one has Δϑ pc > 0, while for u > u c the indicator Δϑ pc becomes negative. At the critical displacement the angle at the centre is equal to the angle at the vertex, and then Δϑ pc = 0. As observed in [START_REF] Dell'isola | Plane Bias Extension Test for a continuum with two inextensible families of fibers: a variational treatment with Lagrange Multipliers and a Perturbation Solution[END_REF], the configuration corresponding to the critical displacement value remains stable, but note that, in general, stability issues concerning systems with a complex microstructure can be difficult. Useful tools for studying such problems can be found, for example, in [START_REF] Luongo | A unified perturbation approach to static/dynamic coupled instabilities of nonlinear structures[END_REF].

Quadratic energy model

Plots of the relative rotation ϑ 1 + ϑ 2 between the fibres for the values of the end displacement 0.65, √ 3 -1 and 0.8 are presented in Fig. 10. The plots must be compared with those of Fig. 7. The constitutive parameter α is set at 1. Also, in this case, the rotation is discontinuous along the diagonals of the specimen. Contrary to the first gradient shear energy model, in this case there is no change in the sign of the difference between the angles of the fibres at the centre and at the vertex of the rigid triangle of the specimen.

Indeed, the indicator angle Δϑ pc (Fig. 11a) is always positive for every value of the imposed displacement u ∈ [0, 0.8]. By differentiating the computed energy with respect to the imposed displacement, we obtain the resulting reaction on the loaded side (Fig. 11b). This curve is characterized by a hardening effect due to the alignment of the inextensible fibres. As the displacement approaches the limit value 2( √ 2 -1), the end reaction tends to diverge.

Finally, in Fig. 12, the deformed configurations for the imposed displacements 0.65, 0.73 and 0.8 are shown. The central part of the specimen gets thinner and thinner, while its boundaries remain straight.

Quartic energy

The quartic energy density g 1,3 is considered next. As was observed at the end of Sect. 2.3, while the quadratic energy can be viewed as a truncated power series expansion of the energy quadratic on the shear, the quartic energy simulates a different network behaviour. Despite this, the results of the simulation are very similar to the those in the previous case, as shown by the plot of the relative rotations presented in Fig. 13. The indicator angle Δϑ pc also remains in this case always positive (Fig. 14a). Moreover, Fig. 15 proves that also for any combination of the quadratic and quartic energy modes, model g 1,3 , the indicator angle Δϑ pc does not change signs. The black line in the plot corresponds to the values α 1 = 1 and α 2 = 1.

Figure 14b summarizes the force-displacement curve obtained for the first gradient models examined. Again, the quartic energy model gives a different response than the other two.

In conclusion, for all the first gradient models considered, the distribution of the rotation is discontinuous. Only with the shear energy model is a change in the sign of Δϑ pc in the deformation for a critical displacement u c observed.

Pure second gradient energy models

In this section we examine the results for the bias extension test in the case where the energy of the material is completely governed by the second gradient deformation. The model relates to the case where no energy is associated to the shear deformation between fibres, and the internal energy is completely due to the bending deformation that

= u 2 = √ 3 -1, c u 1 = u 2 = 0.8
arises in them. The main effect would be to regularize the discontinuities found in the deformation of the sheet with the first order energy models. The rate of convergence of the square of the energy error for the second gradient model is lower than for the first gradient model (Fig. 16), since a P 0 -interpolation for the rotations was used, and the second gradient deformations were obtained using a collocation procedure. As a consequence, the accuracy level obtained in the present case is also smaller than in the previous first gradient case. As for the two different second gradient energies g 2,1 and g 2,2 , the accuracy level is the same (Fig. 17).

It is well known that in second gradient models a characteristic length determines the thickness of the boundary layer. In the presented numerical simulation, we selected a mesh size h fine enough to accurately describe this phenomenon, i.e. much smaller than the observed boundary layer. As is clear from a comparison between Fig. 18b andc, a coarse mesh is not able to reproduce the onset of boundary layers without significant jumps. For a detailed discussion of the theoretical prevision of the thickness of the boundary layer in this type of structure, see [START_REF] Placidi | Identification of two-dimensional pantographic structure via a linear D4 orthotropic second gradient elastic model[END_REF].

Second gradient shear energy (S-model)

The expression g 2,1 of the deformation energy is considered for this model. The plots of the relative rotations Δϑ i j for free end displacements of 0.4, 0.6 and 0.8, presented in Fig. 18a-c respectively, show that in this case the solution does not present discontinuities in the distribution of the deformation (Fig. 19). The deformation of the long side is not straight, and smooth transition zones appear between the different regions. A neck clearly appears in the central part of the specimen. For values of the imposed displacement between 0.6 and 0.8, the relative rotation of the fibres in the neck reaches a limit value of π/2, so that the width of the specimen cannot be reduced any further. The phenomenon recalls the familiar necking in ductile plastic materials whose deformation is driven exclusively by shear. The latter phenomenon can also be observed from the plot of Fig. 17a, where the results of the model g 2,1 are indicated by solid thick lines. It can be observed that the relative rotation at the centre of the specimen reaches its limit value for an end displacement around 0.7.

Second gradient quadratic energy (Q-model)

The energy functional is given in this case by the expression g 2,2 . Also, in this case, there are no discontinuities in the rotation field, but the transition is less smooth than in the previous case, as can be observed from Fig. 20, where the scalar field Δϑ i j for three values of the free end displacement is shown. With respect to the model g 2,1 , smaller deformations are obtained, as can be observed in Fig. 17a, where the relative angles at points c and p for the presented model are indicated by solid thin lines. However, even in this case, in the central part of the specimen, the limit relative rotation Δθ = π/2 is reached (Fig. 20c), although at a value of the end displacement greater than in the case of the energy based on the gradient of the shear. The plot of the reaction vs. the end displacement is presented in Fig. 17b for both second gradient models. For comparison, the results obtained for the first order models are also reported. The second gradient solution is characterized by the same initial stiffness as the first gradient solution, but the stiffness of the specimen strongly increases for the highest values of the displacement, while in the first gradient solution there is a sharper transition to the asymptotic value.

First and second gradient energy models

In this subsection, combinations of the first and second gradient energy models are considered, with the aim of analysing the influence of the second gradient regularization on the solution of the bias extension test. The following four cases are examined:

-Model (SS): First gradient shear energy with second gradient shear energy:

F SS = B (α g 1,1 + β g 2,1
)dB;

-Model (QQ): First gradient quadratic energy with second gradient quadratic energy:

F Q Q = B (α g 1,2 + β g 2,2
)dB;

-Model (SQ): First gradient shear energy with second gradient quadratic energy:

F S Q = B (α g 1,1 + β g 2,2
)dB;

-Model (QS): First gradient quadratic energy with second gradient shear energy:

F QS = B (α g 1,2 + β g 2,1 )dB.
The parameters α and β are used as weights. Notice that changing the weights but leaving their ratio constant does not modify the deformation of the sheet for a given end displacement, only the reaction changes. Therefore, in what follows, only the deformation obtained with the four models is investigated, while the force-displacement relation is not reported.

The deformation of the specimen relative to the free end displacement Fig. 18 Pure second gradient shear energy model (with β = 1): relative rotations ϑ 1,i j + ϑ 2,i j for several imposed displacements: a

u 1 = u 2 = 0.4, b u 1 = u 2 = 0.6, c u 1 = u 2 = 0.8
The transition from the second gradien-dominated deformation, with a smooth continuous variation of the rotations between the parts of the specimen characterized by different regimes, to the first gradient-dominated deformation, with sharp discontinuities along the borders of the different regions of the specimen, is apparent in both cases. It can be observed that the second gradient has a strong regularizing effect, the first order solution with sharp transitions being recovered only when the weight α is at least one order of magnitude greater than β. Fig. [START_REF] Pipkin | Plane traction problems for inextensible networks[END_REF] Pure second gradient shear energy model (with β = 1) -reference and deformed configurations for several imposed displacements:

a u 1 = u 2 = 0.4, b u 1 = u 2 = 0.6, c u 1 = u 2 = 0.8
In the following paragraphs, for the four cases considered, the rotation fields ϑ 1,i j + ϑ 2,i j are reported. They are calculated for the end displacement u 1 = u 2 = 0.6 and for three values of the weights. The differences among the energy models are better understood by examining the indicator angles Δϑ c , Δϑ p and their differences, which relate to the deformation of the central section of the specimen, i.e. the one that presents greater differences in the numerical experiments performed.

First and second gradient quadratic energy: QQ-model

Figure 22 shows the relative angle of rotation for the cases (α = 5, β = 1) and (α = 50, β = 1), enlightening the existence of boundary layers among the regions of the specimen, whose width reduces as the weight of the first gradient energy increases.

It is interesting to examine the plots of the indicator angles, reported in Fig. 23 grouped in three ranges: the second gradient-dominated range, the transition range and the first gradient-dominated range. In the second gradientdominated range, characterized by the weights β α, the second gradient energy predominates, so that smooth boundary layers are observed, and the angle Δϑ p at the corner of the fixed zone of the specimen is always grater than the angle Δϑ c at the centre of the specimen, that is, the centre of the specimen deforms more than the extremities, so that the difference Δϑ pc is always negative (Fig. 23a). In the transition range, the weights have comparable values (actually α is one order of magnitude greater than β), Δθ pc ≈ 0, the two angles become comparable, and there is an inversion in the sign of Δϑ pc (Fig. 23b). In the third first gradient-dominated range, α β, and the indicator angle Δθ pc becomes positive (Fig. 23c). In no case does the latter angle show a change of sign during the elongation.

First and second gradient shear energy: SS-model

Analogously to the previous case, Fig. 24 shows the rotation field Δϑ i j for the same values of α and β and for the same displacement previously considered. The solution in the proximity of the borders of the regions presents a more regular character with respect to the previous case. That is, the shear energy appears to have a greater regularizing effect with respect to the quadratic energy.

The indicator angles are examined in Fig. 25 for the same combination of weights used in the case of model Q Q. In this case for β > α not only is Δϑ pc always negative (Fig. 25a), but in addition the angle at the centre of the specimen Δϑ c becomes zero when the end displacement is between 0.6 and 0.7, indicating that the fibres at the centre are getting aligned. In the transition range Δϑ c is still negative, even though the angle at the centre Δϑ c does not become zero (Fig. 25b). In the first gradient-dominated range Δϑ pc first changes its sign, and then the first order solution for which Δϑ pc goes through the zero is recovered (Fig. 25c).

Shear first gradient energy with quadratic second gradient energy: SQ model

The rotation fields are reported in Fig. 26. The indicator angles Δϑ p and Δϑ c and their difference are represented in Fig. 27a-c, analogously to what was done earlier. In this case, either in the second gradient-dominated range and in the transition range the angle at the centre Δϑ c is always smaller than Δϑ p , but the difference is smaller than in the previous case (model SS). In the first gradient-dominated range, the two angles overlap, and the difference Δϑ pc changes its sign during the evolution of the test, as was observed in model S with shear first gradient energy.

In comparison with model SS it appears, then, that the quadratic second gradient energy has less of an ability to regularize the solution than the second gradient shear energy. at the centre of the specimen reduces to zero for a value of u 1 = u 2 between 0.6 and 0.7. In the transition range, Δϑ pc is always negative and changes its sign only for the highest values of α but, contrary to the previous case, never becomes zero during the elongation.

Conclusions

The deformation of a sheet of woven fabrics composed of two sets of orthogonal inextensible fibres was numerically analysed. To intrinsically enforce the inextensibiity constraint, the rotations of the fibres were used as degrees of freedom of the model. An intrinsic reference frame was introduced, with the axes running in the directions of the two families of inextensible fibres. Using Pipkin's result stating that the deformation field can be decomposed into the sum of two functions, each one depending only on one of the intrinsic coordinates, an expression for the boundary conditions in terms of the rotation fields was furnished. The problem was then formulated as a constrained minimization, with equality constraints enforcing the considered boundary conditions. A suitable piecewise constant interpolation of the two independent rotation fields was used. Three classes of cases were analysed, namely models with energy depending only on the first gradient deformation (which actually reduces to the shear between the fibres), models whose energy depends only on the second gradient deformation, and mixed models. It was shown that the results obtained with first gradient models reproduce the theoretical deformation field, which presents lines of discontinuity on the rotations. The energy based on a second observations can be made. The transition between the different zones of the specimen becomes increasingly sharper as the parameter α grows, that is, as the first gardient energy becomes dominant. For comparison, Fig. 26 relative to the first gradient shear energy should be examined.

The second observation is that the various models trested mainly differ in the deformation pattern obtained in the central zone of the specimen. As observed earlier, only in models where the first gradient shear energy dominates is there a change in the sign of Δθ pc in the plot of the rotations in the central region. This can be observed from Fig. 31, which collects the results for the indicator angles for most of the models examined.

These observations appear to be of great relevance for the experimental characterization of the internal energy of woven composites with a simple bias extension test. More complex cases, involving more general external actions, a richer geometry of the reference configuration or the generalization to 3D homogenized continua formed by inextensible fibres, will probably require suitably refined numerical tools, such as, for instance, isogeometric analysis or enriched finite-element analysis, such as those used in [START_REF] Cuomo | A varational model based on isogeometric interpolation for the analysis of cracked bodies[END_REF][START_REF] Turco | Elasto-plastic analysis of Kirchhoff plates by high simplicity finite elements[END_REF][START_REF] Cazzani | Constitutive models for strongly curved beams in the frame of isogeometric analysis[END_REF][START_REF] Cazzani | Isogeometric analysis: a powerful numerical tool for the elastic analysis of historical masonry arches[END_REF][START_REF] Cazzani | Isogeometric analysis of plane-curved beams[END_REF] and [START_REF] Cuomo | A varational model based on isogeometric interpolation for the analysis of cracked bodies[END_REF][START_REF] Ciancio | On inter-element forces in the FEM-displacement formulation, and implication for stress recovery[END_REF][START_REF] Ciancio | Crack opening at corner nodes in FE analysis with cracking along mesh lines[END_REF].

Finally, the systematic experimental study of the considered structures is of course crucial for further progress. Experimental studies on the subject are carried out in [START_REF] Dell'isola | Designing a light fabric metamaterial being highly macroscopically tough under directional extension: first experimental evidence[END_REF], and new comparisons of the obtained numerical results with experimental ones will be analysed in future works.

Fig. 1

 1 Fig. 1 First gradient energy models: a first gradient energy functions, b conjugated shear stress

Fig. 2

 2 Fig. 2 Set-up of bias extension test: a geometry and subregions, b mesh (for n e = 6)

Fig. 3

 3 Fig. 3 Assumed piece-wise constant interpolations of rotation fields according to Rivlin-Pipkin constraint: a θ 1,i j , b θ 2,i j , c Θ i j = θ 1,i j + θ 2,i j

Fig. 4 Fig. 5 Fig. 6

 456 Fig.[START_REF] Boisse | Composite reinforcements for optimum performance[END_REF] Photograph of 3D-printed pantographic specimen under bias extension test. The regions Δ 00 , Δ 01 , Δ 10 , Δ 23 , Δ 32 and Δ 33 are highlighted (see also Fig.2a). It is clear that the regions Δ 00 and Δ 33 are undeformed

2 Fig. 7 Fig. 8 Fig. 9

 2789 Fig.[START_REF] Dell'isola | A two-dimensional gradient-elasticity theory for woven fabrics[END_REF] First gradient shear energy model -relative rotation ϑ 1,i j + ϑ 2,i j for the specimen 1:3 and α = 1 for several imposed displacements:a u 1 = u 2 = 0.65 < u c , b u 1 = u 2 = u c = √ 3 -1, c u 1 = u 2 = 0.8 > u c

2 Fig. 10 8 Fig. 11 Fig. 12

 21081112 Fig.[START_REF] Dell'isola | Elastic pantographic 2D lattices: a numerical analysis on the static response and wave propagation[END_REF] First gradient quadratic energy model -relative rotation ϑ 1,i j + ϑ 2,i j for specimen 1:3 with α = 1 for imposed displacements:a u 1 = u 2 = 0.65, b u 1 = u 2 = √ 3 -1, c u 1 = u 2 = 0.8

2 Fig. 13 8 Fig. 14

 213814 Fig. 13First gradient quartic energy model -relative rotation, ϑ 1 i j + ϑ 2,i j , for free end displacement u 1 = u 2 = 0.8

Fig. 15 Fig. 16

 1516 Fig. 15 Quartic and quadratic first gradient energy models -indicator angle Δϑ pc for several {α 1 , α 2 } constitutive parameters

Fig. 17

 17 Fig. 17 Pure second gradient energy models (with β = 1): a indicator angles Δϑ p , Δϑ c and the difference Δϑ pc . (Thick line: second gradient shear model, thin line: second gradient quadratic model); b force-displacement curves

u 1 =

 1 u 2 = 0.6 for model SS are shown in Fig. 21. In each figure, six cases are considered, namely the first gradient energy model, the combinations (α = 50, β = 1), (α = 25, β = 1), (α = 5, β = 1), (α = 1, β = 1), and the pure second gradient energy model. In the case of model Q Q the deformed configurations are qualitatively the same but for the same values of the considered weights present a smaller necking effect with respect to model SS.

2 Fig. 20 Fig. 21

 22021 Fig. 20 Pure second gradient quadratic energy model (with β = 1) -relative rotation ϑ 1,i j + ϑ 2,i j for several imposed displacements: a u 1 = u 2 = 0.4, b u 1 = u 2 = 0.6, c u 1 = u 2 = 0.8

4. 4 2 Fig. 22 Fig. 23 2 Fig. 24 Fig. 25 2 Fig. 26 Fig. 27 2 Fig. 28

 4222232242522627228 Fig. 22 First gradient quadratic energy with second gradient quadratic energy (QQ model) -relative rotation ϑ 1,i j + ϑ 2,i j for fixed displacement u 1 = u 2 = 0.6 and β = 1 for the two values of α parameter: a α = 5, b α = 50

Fig. 29 Fig. 30 Fig. 31

 293031 Fig. 29 Indicator angles for QS model: a second gradient-dominated range α = 1, β = 10 2 ; b transition range α = 50, β = 1; c first gradient-dominated range α = 10 3 , β = 1

Acknowledgments

The authors gratefully acknowledge the financial support of the International Research Center on the Mathematics and Mechanics of Complex Systems (M&MoCS). The interesting discussions with the participants of the workshop Computational Mechanics of Generalized Continua and Applications to Materials with Microstructure [Catania, 29-31 October 2015, Scuola Superiore di Catania, and the François Cosserat-Tullio Levi Civita International Laboratory associated to CNRS (LIA Coss&Vita)] were very fruitful and greatly influenced this work.