S Lanteri

C Scheid

N Schmitt

J Viquerat

Discontinuous Galerkin Time Domain Methods for Nonlocal Dispersion Models and Electron Beam Modeling in the Context of Nanoplasmonics

We present recent advancements of the development of our DGTD solvers for computational nanophotonics, particularly for metallic nanostructures irradiated by laser pulses or electron beams. By this means, we firstly discuss the numerical treatment of a nonlocal dispersion model for the electron gas which is necessary for structures in the regime of 2 nm to 25 nm. Subsequently, we deal with the modeling of electron beams in 3D simulations for e.g. electron energy loss spectroscopy or cathodoluminescence.

INTRODUCTION

The relatively new field of nanophotonics has gained significant importance in the last twenty years. Mainly due to technological advancements on the fabrication side that do now allow the controlled manufacturing of detailed nanostructures for photonic experiments and applications.

Computational nanophotonics is mainly concerned with the accurate modeling of the appearing physical effects in the frequency range of interest and the subsequent solution of the modeling equations. This allows access to detailed understanding of microscopic effects which are fairly often extremely hard to obtain via measurements. Hence, appropriate physical and mathematical models allow systematic optimization, uncertainty quantification, and many other kinds of studies that can be performed at relatively low costs compared to experiments.

Our here presented contribution is the continuation of the presented work in [START_REF] Lanteri | Numerical Modeling of Light/Matter Interaction at the Nanoscale with a High Order Finite Element Type Time-domain Solver[END_REF]. For the considered geometry and frequency regime, we model electromagnetics by macroscopic Maxwell's equations with additional material equations for the electron gas of the metal structure. Depending on the geometrical size different models have to be considered. We present recent results of time domain simulations of a nonlocal dispersion (also known as nonlocal Drude) model for structures in the regime of approx. 2 nm to 25 nm. The second part of this work deals with the modeling of electron beams for structures in the regime of approx. 20 nm to 500 nm where local dispersion (e.g. Drude, Drude-Lorentz) models are appropriate. In order to efficiently and precisely solve the above-mentioned equations we use a special time domain Finite Element (FE) method. This method is referred to as Discontinous Galerkin (DG) method. The DG method is easily parallelizable, allows local mesh refinement, conforming meshes, curvilinear elements, etc. and can thus precisely handle large structures in reasonable computational times on modern cluster computers.

DISCONTINOUS GALERKIN TIME DOMAIN METHOD

The DG Time Domain (DGTD) method is based on a local FE formulation with particular basis functions [START_REF] Fezoui | Convergence and stability of a discontinuous Galerkin time-domain method for the 3D heterogeneous Maxwell equations on unstructured meshes[END_REF]. This leads to discontinuous solutions in contrast to classical FE methods. Due to the local formulation of DGTD methods only local mass matrices are inverted which yields a significant advantage by means of computational efficiency in comparison to continuous FE methods. Hence, DGTD methods can combine high order basis functions with explicit time stepping and unstructured meshes where the latter is particularly important for complex geometries.

NONLOCAL DISPERSION MODEL

Nonlocal dispersion models model the nonlocal nature, i.e. mutual electron interaction in the electron gas for metallic nanostructures. In contrast to local models (Drude, Drude-Lorentz,...), nonlocal models take the microscopic electron interaction directly into account and allow additional solutions by means of electron density waves that can travel inside the bulk. [START_REF] Jauho | Unusual resonances in nanoplasmonic structures due to nonlocal response[END_REF][START_REF] Schmitt | A DGTD method for the numerical modeling of the interaction of light with nanometer scale metallic structures taking into account non-local dispersion effects[END_REF] The underlying equa- tions are:

∇ × H = ∂ t D, ∇ × E = -µ 0 ∂ t H, (1)
∂ t D = ε 0 ε ∞ ∂ t E + J, ω 2 p ε 0 ∂ t E = ∂ tt J + γ∂ t J -β 2 ∇(∇ • J), (2)
where ( 2) is a linearized nonlocal Drude model. Here, H, D, E, J, µ 0 , ε 0 , ε ∞ , ω p , γ, and β respectively are the magnetic field, electric displacement, electric field, polarization current, vacuum permeability, vacuum permittivity, relative permittivity at infinite frequencies, plasma frequency, damping constant, and the quantum parameter.

We have successfully implemented equations ( 1) and ( 2) in our DG time domain framework. Figure 1 shows the electric field distribution for a single and coupled infinitely long nanowire.

ELECTRON BEAM MODELING

The second part of our contribution addresses the modeling of electron beams in the context of nanoplasmonics. Electron beams traveling in the vicinity or even penetrating plasmonic structure excite plasmons. Microscopy techniques like Electron Energy Loss Spectroscopy (EELS) and Cathodoluminescence (CL) exploit this effect in order to measure/map plasmonic mode patterns [START_REF] Garcia De Abajo | Optical excitations in electron microscopy[END_REF]. We have recently implemented electron beams and the modeling of EELS as well as CL into our DG time domain framework. First simulation results will be presented.

Figure 1 :

 1 Figure 1: Fourier transformed field plot of |E| at ω = 1.1963 ω p . The left plot shows the bulk plasmon mode for an single infinitely long nanowire and the right plot two coupled wires. The coupling of the two wires in terms of the mode pattern is obvious.

ACKNOWLEDGMENT

The authors gratefully acknowledge support from the University of Nice-Antipolis which funds N. Schmitt.


M Cuomo
email: mcuomo@unict.it

F Dell'isola

L Greco

N L Rizzi

First versus second gradient energies for planar sheets with two families of inextensible fibres: investigation on deformation boundary layers, discontinuities and geometrical instabilities

Keywords: Second gradient continua, inextensible networks, Plane fabrics, Bias test

We present a homogenized model for the analysis of a 2D continuum with two straight families of inextensible fibres embedded in it. The kinematics of the continuum is analyzed and, motivated by phenomenological observations, it is assumed that the strain energy depends on the shear deformation of the fibres and on their bending curvature. It is shown that in order to account for the latter deformation it is necessary to introduce second gradient strains. The problem is formulated as a nonlinear constrained minimization, after introducing a suitable discretization of the domain. Some deformation processes are simulated using different constitutive hypotheses, comparing the predictions obtained assuming the presence of only first gradient or second gradient deformations, or a combination of both. It is found that the first gradient model leads to the presence of discontinuities in the rotation of the fibres, while the second gradient model regularizes these discontinuities by means of boundary layers. In particular in some deformation processes an instability of geometrical nature is observed when the second gradient model is used, that can be suppressed by the first gradient contribution.

Introduction

The characterization of the mechanical behaviour of woven fabrics is of paramount importance in the design of mechanical elements and of the material itself, in order to take advantage of all the peculiar properties of the composite. With respect to traditional composites, woven fabrics present dimensional stability, balanced properties in the reinforcement plane, high impact resistance, good drapability, so that they are suitable to be formed as doubly curved components. In this work we refer to 2D fabrics, obtained interlocking two families of fibres. This is the most common type of fabric used in the engineering practice, and owes its mechanical properties to the interaction between the fibres. While extension in the fibre direction is almost negligible, due to the large elastic modulus of the fibres (often carbon, glass or other high performance materials are used), shear deformation, i.e. the ability of the fibres to undergo relative rotation, can usually be large [START_REF] Lee | Non-orthogonal constitutive model for woven composites incorporating tensile effect on shear behavior[END_REF]. The latter mode of deformation is mostly responsible of the easiness of the fabric to modify its shape and to adapt to the final form required in the design. Friction and slippage of fibres is also a major issue in plane fabrics, when the deformation becomes large. In this context, some significant studies are presented in [START_REF] Nadler | A model for frictional slip in woven fabrics[END_REF][START_REF] Giorgio | Multi-scale concrete model with rate-dependent internal friction[END_REF].

Disregarding the latter effects, that introduce energy dissipation, a fabric characterized by shear deformation alone, under the action of external loads, would develop sharp variations of the rotation along lines, giving rise to kinks in the deformation pattern. Indeed this feature is exploited in simplified models for interpreting the results of experimental tests [START_REF] Dell'isola | A two-dimensional gradient-elasticity theory for woven fabrics[END_REF][START_REF] Dell'isola | Plane bias extension test for a continuum with two inextensible families of fibers: a variational treatment with Lagrange multipliers and a perturbation solution[END_REF][START_REF] Cuomo | Simplified analysis of a generalized bias test for fabrics with two families of inextensible fibres[END_REF]. A careful examination of the deformation pattern, as done in [START_REF] Dell'isola | Pantographic 2D sheets: Discussion of some numerical investigations and potential applications[END_REF] reveals, however, that such sharp kinks do not occur, since the fibres bear some finite flexural rigidity, so that bending is observed.

Apparently, shear and bending are the main deformation modes of woven fabrics, at least for plane deformation and for moderately large displacements. The overall response of the fabric depends on a competition between these two modes. If the shear stiffness prevails, sharp transition zones appear in the element, and eventually the deformation locks. On the contrary, if the bending stiffness prevails, a smoother deformation is obtained, but the element, or part of it, can present floppy modes.

In the work we analyze the relative influence of shear and bending stiffness on the plane deformation of a bidirectional tissue, disregarding all other mechanisms of deformation. Particularly, in addition to the dissipation mechanisms already mentioned, we disregard the out of plane displacement due to wrinkling, that occurs when some part of the element undergoes compression stresses. This aspect will be the object of future investigations, where the methodologies proposed among others in [START_REF] Rizzi | Initial postbuckling behavior of thin-walled frames under mode interaction[END_REF][START_REF] Rizzi | The effects of warping on the postbuckling behaviour of thin-walled structures[END_REF][START_REF] Aminpour | A one-dimensional continuum with microstructure for single-wall carbon nanotubes bifurcation analysis[END_REF] will be employed.

The source of instability in the fabric is an important issue, and its origin needs to be investigated. In models proposed in the literature for the onset of instability the fabric is considered as a thin plate, and classical von Karman's theory is applied [START_REF] Boisse | Simulation of wrinkling during textile composite reinforcement forming. Influence of tensile, in-plane shear and bending stiffnesses[END_REF][START_REF] Harrison | Characterising the sheartension coupling and wrinkling behaviour of woven engineering fabrics[END_REF]. However, those studies do not consider the presence of direction of quasi inextensibility, nor the in plane bending mechanism of deformation. In the paper it will be analyzed the effect of the two deformation mechanisms (shear and bending) on the occurrence of in plane unstable deformation modes.

Testing methods for the mechanical characterization of fabric concentrate on the measure of the shear energy associated to the relative rotation between the fibres. Specific test have been designed to this end, like the direct shear test, the Picture Frame test, and especially the Bias Extension Test that is largely used in the engineering practice, see [START_REF] Harrison | Shear characterisation of viscous woven textile composites: a comparison between picture frame and bias extension experiments[END_REF][START_REF] Cao | Characterization of mechanical behavior of woven fabrics: experimental methods and benchmark results[END_REF][START_REF] Harrison | Normalization of Shear Test Data for Rateindependent Compressible Fabrics[END_REF][START_REF] Harrison | Normalisation of biaxial bias extension test results considering shear tension coupling[END_REF][START_REF] Härtel | Evaluation of normalisation methods for uniaxial bias extension tests on engineering fabrics[END_REF][START_REF] Lee | Bias-extension of woven composite fabrics[END_REF][START_REF] Launay | Experimental analysis of the influence of tensions on in plane shear behaviour of woven composite reinforcements[END_REF][START_REF] Hivet | A contribution to the analysis of the intrinsic shear behaviour of fabrics[END_REF]. In [START_REF] Cuomo | Simplified analysis of a generalized bias test for fabrics with two families of inextensible fibres[END_REF][START_REF] Greco | In plane shear and bending for first gradient inextesible pantographic sheets: numerical study of deformed shapes and global constraint reactions[END_REF] the latter test, and some modified versions of it, were examined, and it was shown that, in general, both shear and bending deformations influence the response of the specimen. According to the relative stiffness of the two modes of deformation, different response were obtained. However, the results in [START_REF] Cuomo | Simplified analysis of a generalized bias test for fabrics with two families of inextensible fibres[END_REF] were found on the basis of simplified kinematic assumptions, that, as will be shown in the present work, only hold when the deformation tends to a limit value. In general, a non uniform deformation state arises in the specimen. In this paper a more accurate analysis of plane deformation states of rectangular samples will be carried out.

The material will be considered as an homogenized continuum, at each point of which there are two inextensible directions, coinciding with the fibre directions. Two families of straight parallel fibres will be considered, according to the most common types of woven fabrics.

The material will be modeled as a second gradient continuum, identifying the first gradient deformation with the shear strain (extensional strains, if needed, could be added), while the bending deformation of the fibres will be identified with the second gradient deformation. In this way a reduced couple stress model is obtained. Analogous approaches can be founded in [START_REF] Dell'isola | A two-dimensional gradient-elasticity theory for woven fabrics[END_REF][START_REF] Dell'isola | Pantographic 2D sheets: Discussion of some numerical investigations and potential applications[END_REF][START_REF] Dell'isola | The postulations á la D'Alembert and á la Cauchy for higher gradient continuum theories are equivalent: a review of existing results[END_REF]. In this context the symmetry analysis presented in [START_REF] Eremeyev | Local Symmetry Group in the General Theory of Elastic Shells[END_REF][START_REF] Eremeyev | Material symmetry group and constitutive equations of micropolar anisotropic elastic solids[END_REF] can be very useful for characterizing the anisotropic properties of the model. Several homogenization approaches for macro model of woven fabrics can be found in [START_REF] Steigmann | Mechanical response of fabric sheets to threedimensional bending, twisting, and stretching[END_REF][START_REF] Dell'isola | Synthesis of Fibrous Complex Structures: designing Microstructure to Deliver Targeted Macroscale Response[END_REF][START_REF] Rahali | Homogenization à la Piola produces second gradient continuum models for linear pantographic lattices[END_REF][START_REF] Goda | Equivalent mechanical properties of textile monolayers from discrete asymptotic homogenization[END_REF][START_REF] Soubestre | Non-local dynamic behavior of linear fiber reinforced materials[END_REF] or in [START_REF] Pideri | A second gradient material resulting from the homogenization of an heterogeneous linear elastic medium[END_REF][START_REF] Alibert | Second-gradient continua as homogenized limit of pantographic microstructured plates: a rigorous proof[END_REF]. Alternatively micro-models as those proposed [START_REF] Dell'isola | Large deformations of planar extensible beams and pantographic lattices: heuristic homogenization, experimental and numerical examples of equilibrium[END_REF][START_REF] Turco | Hencky-type discrete model for pantographic structures: numerical comparison with second gradient continuum models[END_REF][START_REF] Turco | Large deformations induced in planar pantographic sheets by loads applied on fibers: Experimental validation of a discrete Lagrangian model[END_REF] can be used for the identification of the mechanical parameters.

Summarizing aim of the work is to investigate the relative influence of the shear and bending deformation on the response of the woven fabric. A simple geometry will be investigated, that is a rectangular strip with the fibres symmetrically oriented with respect to the specimen axis, with one end built in.

Kinematics

In this section the representation of the placement of a bi-dimensional sheet with inextensible fibres will be introduced. A theorem due to Rivlin and Pipkin [START_REF] Pipkin | Some developments in the theory of inextensible networks[END_REF][START_REF] Pipkin | Plane traction problems for inextensible networks[END_REF][START_REF] Rivlin | Plane strain of a net formed by inextensible chords[END_REF] will be exploited in order to get an implicit account of the inextensibility constraint.

Constrained placement

Let B 0 be the reference domain occupied by the sheet, having an uniform thickness t (for a woven fabric the concept of thickness has to be intended in an average sense). We consider the case of two families of straight parallel inextensible fibers embedded in the continuum, with directions given by the vectors D 1 , D 2 . It is convenient to assume a reference system with axes parallel to the fiber directions, so that, calling ξ 1 , ξ 2 two coordinates having values in a subset Ξ ⊂ R 2 , the differential of the reference placement is:

dp 0 = dξ 1 D 1 + dξ 2 D 2 . ( 1
)
Here D i = ∂p 0 ∂ξ i are the covariant bases of the manifold described by [START_REF] Lee | Non-orthogonal constitutive model for woven composites incorporating tensile effect on shear behavior[END_REF]. The vectors D i can be taken, without loss of generality, unitary. Calling p the position at a generic configuration, we have

dp = dξ 1 d 1 + dξ 2 d 2 , ( 2
)
with

d 1 = ∂p ∂ξ 1 = D 1 + ∂u ∂ξ 1 , d 2 = ∂p ∂ξ 2 = D 2 + ∂u ∂ξ 2 , ( 3
)
having indicated by u the displacement vector.

Since the fibers are inextensible, the directors cannot stretch but can only rotate. Let's indicate with R 1 (θ 1 ), R 2 (θ 2 ) two rotation matrices, the first anticlockwise, the second clockwise, and by D 1 , D 2 the contravariant basis vectors (unitary if D i are unitary). The directors (3) can then be represented as (see figure 1):

d 1 =R 1 (θ 1 ) D 1 = d 1 • D 1 D 1 + d 1 • D 2 D 2 = cos ( θ 1 + π 2 -α ) D 1 + sin(θ 1 ) D 2 , d 2 =R 2 (θ 2 ) D 2 = d 2 • D 1 D 1 + d 2 • D 2 D 2 = sin(θ 2 ) D 1 + cos ( θ 2 + π 2 -α ) D 2 .
(4) For the rotations to represent a compatible deformation it must be:

∂d 1 ∂ξ 2 = ∂d 2 ∂ξ 1 (5)
that is -sin

( θ 1 + π 2 -α ) ∂θ 1 ∂ξ 2 D 1 + cos(θ 1 ) ∂θ 1 ∂ξ 2 D 2 = -sin ( θ 2 + π 2 -α ) ∂θ 2 ∂ξ 1 D 2 + cos(θ 2 ) ∂θ 2 ∂ξ 1 D 1 (6) Equating the components in the two directions D 1 , D 2 one gets two homogeneous equations in ∂θ 1 ∂ξ 2 , ∂θ 2 ∂ξ 1 : ( cos(α -θ 1 ) cos(θ 2 ) cos(θ 1 ) cos(α -θ 2 ) ) ( ∂θ 1 ∂ξ 2 ∂θ 2 ∂ξ 1 ) = ( 0 0 ) (7)
which admits as solutions only ∂θ 1 ∂ξ 2 = 0, ∂θ 2 ∂ξ 1 = 0, unless the determinant of the matrix on the left hand side of (7) becomes zero, which happens only for θ 1 + θ 2 = α, that is when the fibres superpose to each other, which is a degenerate condition that will be excluded from the analysis.

Therefore the rotation angle θ 1 is a function of ξ 1 only and θ 2 is a function of ξ 2 only. The results generalizes Rivlin's theorem that for a net constituted by two families of orthogonal inextensible fibres the deformation along each of the two families of lines is independent form the other direction.

In this work it will be examined the case of two orthogonal families of fibres, so that (4) specialize to:

d 1 =R 1 (θ 1 ) D 1 = cos(θ 1 ) D 1 + sin(θ 1 ) D 2 , d 2 =R 2 (θ 2 ) D 2 = sin(θ 2 ) D 1 + cos(θ 2 ) D 2 .
(8)

Strain tensors

From equation (2) it follows that the Green's strain tensor is:

E = 1 2 [( 1
d 1 • d 2 d 1 • d 2 1
) -

( 1 D 1 • D 2 D 1 • D 2 1 )] . ( 9
)
The component

E 12 = 1 2 (d 1 • d 2 -D 1 • D 2
) is denoted by γ 2 and represents the shear deformation. With the notation of figure 1

it is γ = cos (α -(θ 1 + θ 2 )) -cos(α) = cos(α) (cos(θ 1 + θ 2 ) -1) + sin(α) sin(θ 1 + θ 2 ) (10)
that, for small relative rotations of the directors, can be linearized as lim

θ 1 +θ 2 →0 γ = sin(α)(θ 1 + θ 2 ) (11)
In the case of orthogonal fibres it is

γ = d 1 • d 2 = sin (θ 1 + θ 2 ) . ( 12
)
whose linearized approximation is lim

θ 1 +θ 2 →0 γ = θ 1 + θ 2 . ( 13
)
Second gradient deformation is also considered in this work. Denoting by ∇E the strain gradient tensor, its only non vanishing components are:

(∇E) 121 = (∇E) 211 = 1 2 sin (α -(θ 1 + θ 2 )) ∂θ 1 ∂ξ 1 (∇E) 122 = (∇E) 212 = 1 2 sin (α -(θ 1 + θ 2 )) ∂θ 2 ∂ξ 2 (14)
In the limit for vanishing relative rotation we get the linearized forms:

lim θ 1 +θ 2 →0 (∇E) 121 = 1 2 sin(α) ∂θ 1 ∂ξ 1 lim θ 1 +θ 2 →0 (∇E) 122 = 1 2 sin(α) ∂θ 2 ∂ξ 2 (15)
The latter expressions for α = π/2 specializes to:

(∇E) 121 = (∇E) 211 = 1 2 cos(θ 1 + θ 2 ) ∂θ 1 ∂ξ 1 (∇E) 122 = (∇E) 212 = 1 2 cos(θ 1 + θ 2 ) ∂θ 2 ∂ξ 2 (16)
whose linearized forms are (∇E)

lin 121 = 1 2 ∂θ 1 ∂ξ 1 , (∇E) lin 122 = 1 2 ∂θ 2 ∂ξ 2 .
The strain gradient components account then for the bending deformation of the fibres, since the line derivatives of the rotations represent the bending curvature of the fibres. However in the exact expression of the strain gradient (16) the measure of the curvature is influenced by the amplitude of the relative rotation of the fibres. The more the fibres tend to overlap the smaller is the bending strain.

Constitutive models

Only elastic deformation is considered in the paper, so that dissipation mechanisms (friction, damage etc.) are ruled out. The hypothesis appears restrictive, but it is sufficient for simulating a monotonic loading process, at least up to fairly large deformation. The strain energy is assumed as the sum of two independent contributions, the first depending on the second invariant of the Green's strain tensor (first energy model), and the second on the strain gradient components, which represents the part of energy absorbed by bending (second gradient energy):

Φ = ∫ B 0 g 1 (γ) dB + ∫ B 0 g 2 (∇γ) dB ( 17
)
For both functions a quadratic form is assumed. Particularly, g 1 is taken as a function of the second invariant of the Green's strain tensor (the only non null invariant), I 2 = γ 2 , that is

g 1 (E) = 1 2 k 1 γ 2 (18)
with k 1 a positive constant having the dimension of a stress. Also g 2 is assumed as a quadratic form of the components of ∇E,

g 2 (∇E) = 1 2 k 2 ( (∇E) 2 121 + (∇E) 2 122 ) . ( 19
)
The stiffness coefficients k 1 , k 2 are independent properties, however, it is possible to give an interpretation for their ratio, attributing k 1 to the shear energy per unit volume needed to deform the fabric, directly proportional to the width l of the fibre, and k 2 to the bending stiffness of a fibre, proportional to the cube of its width, so that one has

k 2 k 1 ∼ l 2 (20)
In [START_REF] Dell'isola | Large deformations of planar extensible beams and pantographic lattices: heuristic homogenization, experimental and numerical examples of equilibrium[END_REF] a sound physical argument has been given for the ratio [START_REF] Hivet | A contribution to the analysis of the intrinsic shear behaviour of fabrics[END_REF] on the basis of an homogenization procedure, together with a more rigorous evaluation. Denoting by l a scale parameter representative of the microscale of the network, not necessarily equal to the width of the fibres, we can take k 2 = k 1 l 2 , so that equation [START_REF] Härtel | Evaluation of normalisation methods for uniaxial bias extension tests on engineering fabrics[END_REF] becomes

Φ = 1 2 k 1 (∫ B 0 γ 2 dB + l 2 ∫ B 0 ∥∇E∥ 2 dB ) (21)
The parameter l will then rule the relative weight of the two contributions of the energy to the overall deformation. For a vanishing l a first gradient material is obtained, while a large l represents a material for which the second gradient (bending) deformation is the most relevant energy storage phenomenon.

Geometrical model

Geometry of the sample

The rectangular strip of figure 2(a) is considered in the paper. Its left edge Σ 1 is built in, while its right edge Σ 2 is free to translate and rotate but is restrained to remain undeformed. The width and the length of the specimen are W and H. A double family of orthogonal inextensible fibres at 45 • degrees with the specimen axis is embedded in the sheet. In figure 2(a) is represented the reference frame used in the numerical simulations.

A length to width ratio 3:1 is used unless otherwise specified. The specimen can be divided in regions, as highlighted in figure 2(b), in which different pattern of deformation can be recognized. Regions ∆ 00 and ∆ 33 , whose sides cannot elongate, remain undeformed, region ∆ 33 undergoes, then, a rigid body motion.

The geometry was already investigated in [START_REF] Cuomo | Simplified analysis of a generalized bias test for fabrics with two families of inextensible fibres[END_REF][START_REF] Dell'isola | Bias extension test for pantographic sheets: numerical simulations based on second gradient shear energies[END_REF] in connection with the Bias Extension test and a modified form of it.

Boundary conditions

Displacements or forces are applied at the free edge Σ 2 of the specimen, and, since the region ∆ 33 remains rigid, those displacements can be expressed as a function of the three degrees of freedom,

u 1 O , u 2 O , ϕ o , of the point O, see figure 3(a)
. In particular it is convenient to express the boundary conditions referring to point Q as:

u Q (θ 1 , θ 2 ) = ūQ ( 22
)
where u Q (θ 1 , θ 2 ) is the displacement of point Q evaluated from the kinematic fields, and ūQ is the displacement of point Q obtained from the imposed displacements at point O as:

ūQ = u O + (R(ϕ O ) -I) (Q -O) (23)
with u O is the displacement vector of point O, R(•) is the anti clockwise unit rotation operator around point O and I is the identity matrix.

Since P is a fixed point, the term u Q (θ 1 , θ 2 ) can be evaluated as

(u Q -u P ) = ∫ Γ (dp -dp 0 )dΓ (24)
being Γ a line joining points Q and P. Substituting the expressions (8) one gets:

∫ Γ (dp -dp 0 )dΓ = ∫ ξ 1 Q ξ 1 P (cos(θ 1 )D 1 + sin(θ 1 )D 2 ) dξ 1 + + ∫ ξ 2 Q ξ 1 P (cos(θ 2 )D 2 + sin(θ 2 )D 1 ) dξ 2 + -(ξ 1 P -ξ 1 Q )D 1 -(ξ 2 P -ξ 2 Q )D 2 (25)

Problem Formulation

The deformation field of the 2D sheet described in section 3 is obtained solving a constrained minimization problem in the rotation angles θ 1 , θ 2 . The objective functional is given by the energy potential [START_REF] Härtel | Evaluation of normalisation methods for uniaxial bias extension tests on engineering fabrics[END_REF], depending only on θ 1 (ξ 1 ) and θ 2 (ξ 2 ), under the inequality constraints on the rotation angle fields, that enforce the condition of non-overlapping rotation of the fibres:

min (θ 1 ,θ 2 ) Φ ( θ 1 , θ 2 , dθ 1 dξ 1 , dθ 2 dξ 2 ) -λ • (u Q (θ 1 , θ 2 ) -ūQ ) , (26)
subjected to : -

π 2 ≤ θ 1 + θ 2 ≤ π 2 (27)
with Φ defined by [START_REF] Greco | In plane shear and bending for first gradient inextesible pantographic sheets: numerical study of deformed shapes and global constraint reactions[END_REF]. The second term in the objective functional enforces the boundary conditions expressed by the condition [START_REF] Dell'isola | The postulations á la D'Alembert and á la Cauchy for higher gradient continuum theories are equivalent: a review of existing results[END_REF], by means of the Lagrangian multiplier vector λ, that represents the reacting force applied at point Q. Notice that in the case the displacement is not assigned, so that ūQ vanishes, the latter term in [START_REF] Dell'isola | Synthesis of Fibrous Complex Structures: designing Microstructure to Deliver Targeted Macroscale Response[END_REF] represents the external work done by a load applied at point Q.

In the case displacement boundary conditions are applied, the reaction T dual of the applied displacement is obtained by means of Catigliano's theorem, see [START_REF] Cuomo | Simplified analysis of a generalized bias test for fabrics with two families of inextensible fibres[END_REF] for details.

The problem is solved discretizing the domain in quadrilateral elements, and considering a piece-wise constant Q 0 -interpolation for the two rotation fields θ 1 (ξ 1 ) and θ 2 (ξ 2 ) on the uniform quadrilateral mesh. The primal kinematical variables are the two rotation fields while the displacement is obtained by means of an integration procedure. The discretized values of the rotations are denoted by θ 1,ij , θ 2,ij , with i = 1, ..., n e and j = 1, ..., 3n e 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 figures 3(b) and 3(c), where are also shown the elements whose rotation are imposed by the boundary conditions. In figure 3(d) is synthesized the distribution of the sum of the two rotations in the regions of the specimen.

Since a piece-wise constant interpolation has been used, the gradients of the rotations cannot be directly evaluated. They have therefore been collocated in the vertices of the mesh, evaluating 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,ij √ 2h e . ( 28
)
where h e = W/n e is the element side length. The second gradient deformation has been calculated in a similar way.

Numerical Investigations

For each of the cases investigated in the paper we will compare the predictions of the model obtained assuming that only the first gradient term (resp. the second gradient term) appears in the energy functional [START_REF] Greco | In plane shear and bending for first gradient inextesible pantographic sheets: numerical study of deformed shapes and global constraint reactions[END_REF]. Furthermore, for each simulation two cases will be considered, that differs for the evaluation of the strain. In the first case, the exact non linear expressions of the strains ( 12), ( 16) will be considered. In the second case the linearized expressions for the strains will be used. Notice that in both cases the strain energy, being quadratic, is convex, so that material instability are ruled out from the model. In the first case, however, the strain displacement relation is non linear, so that a geometric effect may arise. The question is of great importance in the occurrence of an unstable response in the plane.

The results obtained with the first gradient model will be labeled by 1gS, 1gQ, where S, Q indicate that the non linear or the linear expression has been used for the shear. Likewise, the results obtained with the second gradient model will be labeled 2gS, 2gQ.

Bias Extension Test

In this subsection we consider the case that an axial translation is applied to the specimen, in order to simulate a standard Bias Extension Test. In [START_REF] Cuomo | Simplified analysis of a generalized bias test for fabrics with two families of inextensible fibres[END_REF] it has been shown that in this case a limit displacement is reached when in the central zone of the specimen the fibres get aligned (θ 1 + θ 2 = π/2), so no further extension is allowed. The value of the limit displacement is

u lim = (√ 2 -1 ) (H -W ).
The results of the simulation are presented in dimensionless form. The dimensionless displacement is u Q /u lim , ranging in (0, 1). The reaction T is divided by the initial slope S 0 of the reaction-displacement curve, T adim = T /(S 0 u lim ). In this way all the forcedisplacement curves have the same initial tangent, so that they can be directly compared. The dimensionless form of the strain energy used is Φ/(tW 2 k 1 ), where t is the thickness of the sample. Moreover the dimensionless scale parameter α = l/W is also used. Since α ranges in (0, ∞), it has been replaced by α 2 = ζ 1-ζ , so that ζ = 0 for α = 0 (first gradient model), and ζ = 1 for α = ∞ (second gradient model).

The dimensionless plots of the reaction vs. the imposed displacement are presented in figure 4(a) for the two energy models considered and for the two forms of the strain. In all cases the reaction tends to diverge when the displacement approaches the limit value. However the first gradient model shows a sharper neck in the force-displacement curve when the limit displacement is approached, and lower values of the dimensionless reaction for an assigned displacement.

In figure 4(b) is plotted the dimensionless strain energy vs. the applied displacement. The first gradient energy approaches a limit value, when the deformation in the sample reaches its limit configuration, shown in figure 5(b). In it the relative rotation θ 1 + θ 2 in the central zone becomes π/2, while θ 1 + θ 2 = π 4 in regions ∆ 01 , ∆ 10 , ∆ 23 , ∆ 32 . The dimensionless strain energy [START_REF] Greco | In plane shear and bending for first gradient inextesible pantographic sheets: numerical study of deformed shapes and global constraint reactions[END_REF] becomes in this condition equal to 1 if the non linear expression of the strain is used, and equal to 2.159 if the linear expression of the strain is used. In the second gradient energy model, on the contrary, it appears that no limit value of the energy is attained, due to the presence of a boundary layer where the deformation concentrates. As has been observed in previous works ( [START_REF] Dell'isola | Plane bias extension test for a continuum with two inextensible families of fibers: a variational treatment with Lagrange multipliers and a perturbation solution[END_REF][START_REF] D'agostino | Continuum and discrete models for sructures including quasi inextensible elasticae with a view to the design and modeling of composite reinforcements[END_REF][START_REF] Greco | In plane shear and bending for first gradient inextesible pantographic sheets: numerical study of deformed shapes and global constraint reactions[END_REF]), the deformation of the sample in the first gradient model presents sharp discontinuities along the boundaries of the regions, as is highlighted in figure 5(a), 5(b), whatever is the level of the imposed displacement. The presence of jumps in the strain field implies that along the discontinuity lines arise concentrated line forces that equilibrate the difference in the shear stresses. This peculiar phenomenon will be investigated in a future work, where an interpretation of the concentrated line forces will be proposed in analogy with what occurs in cable nets [START_REF] Greco | A procedure for the static analysis of cable structures following elastic catenary theory[END_REF][START_REF] Greco | On the force density method for slack cable nets[END_REF].

In the deformation obtained with the second gradient energy model, instead, this sharp discontinuities are substituted by layers where a smooth transition occurs (figures 6(a), 6(b)). Figure 7 shows that for the first gradient model the deformation, in addition to present sharp discontinuities, is not uniform in the central area of the specimen, but reaches the theoretical piecewise constant deformation for the limit displacement. The plot of the distribution of the energy shows a similar trend. The situation when the second gradient model is used is very different (figure 8). The deformation does not present discontinuities, and even for values of the imposed displacement very close to the limit one a transition layer is present. The distribution of the energy is completely different form the first gradient results. For the latter case the plot of the strain energy density is similar to the distribution of the shear deformation. In the former case, instead, the main contributions to the strain energy tend to concentrate in the transition layers. For values of the imposed displacement close to the limit one almost all the strain energy is stored in these layers, corresponding to the fact that the deformation in the specimen gets close to the piecewise constant distribution found in the first gradient energy, except for the presence of the transition layers. A comparison of the relative rotations of the fibres along the central axis of the specimen is presented in figure 9 for an imposed displacement u = u lim /2 and u = 0.99 u lim . Both cases of the exact and of the linearized kinematics are presented. The presence of a boundary layer on the second gradient deformation as opposed to the jump experienced by the first gradient solution is apparent. The use of the linearized kinematics, for which the second gradient strain is only function of the curvature and not of the relative rotation also, does not appear to modify the shape of the boundary layer. It is interesting to estimate the width of the transition layer, that should be related to the ratio α 2 = (l/W ) 2 of the stiffness constants k 1 and k 2 . However for the pure second gradient model (ζ = 1) this parameter diverges, so that the interpretation of the internal scale is not trivial. To this end in figure 10 has been plotted the rotation θ 1 along the line directed in the ξ 1 direction shown in the figure (along this line the rotation θ 2 is zero), for several values of the parameter ζ ∈ (0, 1), using the exact kinematics and for 4 values of the imposed displacement: u = u lim /4, u = u lim /2, u = 3 u lim /4, u = 0.99 u lim . On the horizontal axis are indicated the number of elements along the selected line. Two observations may be done. First, examining the plot relative to a single step, it can be seen that for ζ = 0 there exists a jump in the value of the rotations, which disappears for any other value of the parameter ζ > 0. Also, apparently, for values of ζ larger than 0.1 it seems that the boundary layer width remains almost constant. Notice that for ζ = 0.3 one gets l ≃ W/3. However, examining the last figure 10(d) relative to a deformation close to the limit one, the boundary layer is still present but gets thinner and thinner. Moreover, as already observed, the deformation in the specimen outside of the layer tends to become constant, more rapidly than it happens with the first gradient solution (ζ = 0).

The deformation of the specimen for u = u lim thus tends in the limit to develop also for the second gradient model sharp kinks, where a diverging bending deformation occurs. It should be observed that numerically it is very difficult to find this limit state, since convergence becomes slower and slower. The plots of figure 11 presents the convergence for the strain energy and for the reaction at several levels of the imposed displacement as function of the number of element along the width of the sample, for the second gradient model. While when the first gradient model is used convergence is always reached with a reasonable discretization, for the second gradient solution it can be observed that convergence slows down dramatically for large values of the imposed displacement.

Shear Test

In this simulation a monotonically growing tangential displacement is imposed to the point O , and hence to point Q, as illustrated in figure 12. The region ∆ 33 translates rigidly, while the central zone of the specimen deforms. The test is similar to the double bending of a cantilever, however, due to the inextensibility constraints, the deformation is quite different. Figures 13 and14 show the evolution of the deformation for the two limit cases that the strain energy depends on only the first (resp. second) gradient deformation. When the first gradient energy model is used it can be observed that, similarly to what happens in the Bias Extension Test, jumps in the rotations of the fibres occur along the boundaries of the triangular regions in the specimen. Consequently, a limit state is reached in which the rotation is piecewise constant, equal to zero in regions ∆ 00 , ∆ 01 , ∆ 33 , ∆ 32 , and to π/2 in the rest of the specimen, so that all the fibres get aligned and no further deformation is possible. The limit displacement is v lim = H -W . At this stage the first gradient strain energy stored in the sample is equal to Φ = 1 2 k 1 γ 2 tH(H -W ), so that its dimensionless value, for the 1:3 sample examined, becomes Φ ad = Φ t k 1 W 2 = 1, if for γ the exact kinematics is used, and 2.467 if the linear kinematics is used.

The deformation obtained with the second gradient model does not show kinks, which are replaced by transition layers. Since the shear stiffness in this model is zero, the sample undergoes larger rotations with respect to the first gradient model, and the fibres tend to overlap along a diagonal line, as can be seem comparing for instance figures 13(c) and 14(c). The plot of the reaction vs. the applied displacement is presented in dimensionless form in figure 15(a), for both energy models and for the exact and linearized forms of the strains. In figure 15(b) is plotted the strain energy as function of the dimensionless applied displacement v/v lim . Differently than in the case of the Bias Extension Test, the second gradient model presents a softening behaviour. For an applied displacement v/v lim ≃ 0.6 the reaction reaches a maximum and then decreases. Subsequently, approaching the limit displacement, the reaction grows again, tending to diverge, as with the other models. The presence of this instability is confirmed in the plot of the strain energy for the second gradient model, that shows a discontinuity on the slope. Notice that the results obtained with the second gradient model but using the linearized strain do not present instability, thus suggesting that the origin of the instability is of a geometrical nature related to the strain-displacement relation.

2 Ξ 1 0 1 2 3 4 0 1 2 3 4 5 (a) v = 1/4 v lim . Ξ 2 Ξ 1 0 1 2 3 4 0 1 2 3 4 5 (b) v = 1/2 v lim . Ξ 2 Ξ 1 0 1 2 3 4 0 1 2 3 4 5 6 (c) v = 3/4 v lim . Ξ 2 Ξ 1 0 1 2 3 4 0 1 2 3 4 5 6 (d) v = v lim .
2 Ξ 1 0 1 2 3 4 0 1 2 3 4 5 (a) v = 1/4 v lim . Ξ 2 Ξ 1 0 1 2 3 4 0 1 2 3 4 5 (b) v = 1/2 v lim . Ξ 2 Ξ 1 0 1 2 3 4 0 1 2 3 4 5 6 (c) v = 3/4 v lim .
In order to further investigate the onset of the instability, a series of simulations have been carried out on samples with increasing slenderness H/W . The results are summarized in figure 16, where the end reaction has been plotted vs. the displacement normalized with respect to the limit value v lim = H -W . The limit loads decreases with increasing slenderness of the sample, strengthening the attribution of the instability to a geometrical phenomenon. It is interesting to notice that the limit load occurs for all the slenderness ratios investigated at the same value of the normalized displacement. Graphs of the relative rotations along the central line of the specimen for four values of the applied displacement are presented in figure 17 for several values of the parameter ζ. The contour plots also presented in the figure refer to the relative rotation evaluated with the model 2gS, ζ = 1. In this example the results obtained changing the parameter ζ differ substantially. Notice that the thickness of the boundary layers can be appreciated from the contour plots reported in the figure for the second gradient model.

The distribution of the strain energy shows a similar trend to what found in the Bias Extension Test. While for the first gradient energy model the strain energy tends to become uniform in the central zone of the specimen, it tends instead to concentrate along the transition regions in the second gradient model, while vanishing in the remaining part of the sample (figures [START_REF] Lee | Bias-extension of woven composite fabrics[END_REF][START_REF] Launay | Experimental analysis of the influence of tensions on in plane shear behaviour of woven composite reinforcements[END_REF].

In order to better investigate the behaviour of the sample it has been evaluated its response in the case both the first gradient and the second gradient strains contribute to the strain energy, that is using a finite scale parameter α = l/W . The results are presented in figures 20(a) and 20(b) for the reaction and the strain energy. In these plots the reaction has been adimensionalized dividing by the initial slope of the force displacement curve, while the dimensionless strain energy as been defined as (1 -ζ)Φ/(k 1 tW 2 ), in order to avoid the divergency of the scale parameter. For values of ζ > 0.2, (α > 0.5) the load-displacement curve presents a limit point. The first gradient strain appears then to regularize the instability that occurs when the second gradient strain only is considered. The relative influence of the first and second gradient strain on the plane deformation of a pantographic model was already examined in [START_REF] Scerrato | On the effect of shear stiffness on the plane deformation of linear second gradient pantographic sheets[END_REF].

End Rotation Test

A clockwise rotation of the free edge around point O is now applied. Figure 21 shows a sketch of the simulation. Also in this case a limit rotation exists, when the fibres in region ∆ 23 get aligned but, since there is no symmetry in the deformation, the limit rotation cannot be predicted, as was already indicated in [START_REF] Cuomo | Simplified analysis of a generalized bias test for fabrics with two families of inextensible fibres[END_REF]. In the case examined here, the limit rotation has been found to be π/3. The deformation of the sample, presented in figure 22, 23 for the first and second gradient models with the exact kinematics, shows that indeed the specimen axis bends. Also in this case the force displacement curve presents an instability for the 2gS model at a rotation about equal to π/6. In this case also the 1gS model presents a softening in the stiffness before the limit rotation is reached, At this stage, as can be observed from figure 23(c), a concentration tends to develop along a diagonal line in the specimen. Figures 25,[START_REF] Dell'isola | Synthesis of Fibrous Complex Structures: designing Microstructure to Deliver Targeted Macroscale Response[END_REF] highlight this phenomenon. Comparing them it can again be observed the presence of boundary layers that characterize the second gradient solution, where the energy tends to concentrate, as opposite to the solution of the first gradient model, that presents also in this case jumps in the values of the rotations. A phenomenological explanation of this behaviour has to be sought in the competing actions of the shear deformation, whose rigidity decreases with increasing deformation, and of the inextensibility constraint, whose effect becomes dominant when the fibres tend to align. The use of a combined model eliminates the limit point for a value of ζ roughly smaller than 0.5 (α = 1), see figure 27. Figure 28 shows the evolution of the relative rotation of the fibres along the central line of the sample as the parameter ζ grows from 0 to 1.

Ξ 2 Ξ 1 0 1 2 3 4 0 1 2 3 4 (a) ϕ o = π 12 . Ξ 2 Ξ 1 0 1 2 3 4 0 1 2 3 4 (b) ϕ o = π 6 . Ξ 2 Ξ 1 0 1 2 3 4 0 1 2 3 4 (c) ϕ o = π 4 . Ξ 2 Ξ 1 0 1 2 3 4 0 1 2 3 4 (d) ϕ o = π 3 .
Figure 22: End Rotation Test -first gradient model (1gS-model): reference and deformed configurations for the four values of the applied rotation.

Ξ 2 Ξ 1 0 1 2 3 4 0 1 2 3 4 (a) ϕ o = π 12 . Ξ 2 Ξ 1 0 1 2 3 4 0 1 2 3 4 (b) ϕ o = π 6 . Ξ 2 Ξ 1 0 1 2 3 4 0 1 2 3 4 (c) ϕ o = π 4 . Ξ 2 Ξ 1 0 1 2 3 4 0 1 2 3 4 (d) ϕ o = 0.95 × π 3 .

Cantilever

In the final test all degrees of freedom of the free edge are together activated. A monotonically increasing tangential conservative load is applied at point Q, as indicated in figure 29(a) The deformation of the specimen is represented in figures 30 for the first gradient model, and 31 for the second gradient model. The steps at which the deformation is represented are indicated by points A,...D on the load displacement plot of figure 29(b). In the figure is plotted the value of the applied force vs. the component of the displacement of point Q in the direction of the force. In the first gradient model the deformation initially concentrates near the built in end of the specimen, until the applied force aligns with the deformed axis of the specimen. From this instant on the specimen deforms in a way similar to an Extension Test, with the formation of lines of discontinuity. In the case of the second gradient model, being the deformation ruled by bending of the fibres, a boundary layer quickly appears on the top side of the sample, which is reflected in the softening of the load displacement curve that occurs between points A and B.

Conclusions

In the paper we have examined the plane deformation of a bidimensional tissue constituted by two orthogonal families of fibres that have been considered inextensible. The kinematics of the model has been obtained on the basis of a theorem of Pipkin, that allows to take as independent variables describing the deformation of the continuum the rotation fields of the fibres. It has been recognized that the mechanical beahviour of such a material depends on the shear deformation (relative rotation of the fibres) and on the bending deformation of the fibres. The first deformation model has been associated with the first gradient shear strain measure, while the second deformation is related to the gradient of the strain. A numerical model has been used to investigate the response of a planar sheet to some simple tests. A rectangular specimen with one end built in has been analysed. Different constitutive hypotheses have been used, that is it has been examined the response of the specimen in the case the fabric behaves like a first gradient or like a second gradient material, or in the case the material response is a combination of the two.

It has been found that in the case a first gradient material model is employed discontinuities in the deformation field are always found, for any value of the applied displacement, since the fibres bear no resistance to bending. The discontinuities disappear if a second gradient material model is employed, substituted by boundary layers whose thickness, however, tends to zero when a limit deformation is approached where the inextensibility effect prevails. The thickness of the boundary layer in the examples analysed remains approximately the same also if a linearized measure of the second gradient deformation is used.

It has been found that for some non symmetric deformation field the second gradient material model presents an in plane instability if a non linear measure of the bending (second gradient) strain is used. This instability is a different phenomenon from wrinkling, that is an out of plane deformation, but can interact with the latter. The introduction of a contribution of the first gradient shear deformation in the model tends to stabilize the response. This can be an useful design suggestion, in order to produce fabrics that have a stable in plane behaviour.

A structural application has been presented with reference to a cantilever, whose deformation field is very different from the one exhibited by a standard elastic cantilever, either a first gradient or a second gradient model is used. This is mainly due to the inextensibility constraint introduced by the fibres. Future studies will examine the links between the macroscopical mechanism of deformation and the microstructure of the fabric. Estimates of the constitutive parameters of the material will be sought using either homogenization procedures or defining a set of fabric parameters defining the microstructure of the composite as done in [START_REF] Misra | Identification of Higher-Order Elastic Constants for Grain Assemblies based upon Granular Micromechanics[END_REF].

Figure 1 :

 1 Figure 1: Definitions of the rotation angles.

Figure 2 :

 2 Figure 2: Reference configuration: (a) Initial configuration and referential inextensible directors, (b) Partitioning and labeling of subregions.

Figure 3 :

 3 Figure 3: Boundary conditions evaluation and distribution of the rotations in the sample (a) evaluation of the boundary conditions, (b) distribution of the angle θ 1 , (c) distribution of the angle θ 2 , (d) regions where the θ 1 + θ 2 get different forms.

Figure 4 :

 4 Figure 4: Bias Extension Test: (a) Dimensionless force versus u Q /u lim and (b) Strain energy versus u Q /u lim for the first and second gradient models, (ζ = 0 and ζ = 1).

Figure 5 :

 5 Figure 5: Bias Extension Test -model 1gS, first gradient energy model: reference and deformed configurations for two values of the displacements.

Figure 6 :

 6 Figure 6: Bias Extension Test -model 2gS, second gradient energy model: reference and deformed configurations for two values of the displacements.

  (a) (ϑ 1 + ϑ 2 ) max = 0.259946. (b) h 2 g 1,max = 0.0000809379. (c) (ϑ 1 + ϑ 2 ) max = 0.539297. (d) h 2 g 1,max = 0.000323304. (e) (ϑ 1 + ϑ 2 ) max = 0.849624. (f) h 2 g 1,max = 0.000692719. (g) (ϑ 1 + ϑ 2 ) max ∼ = π/2. (h) h 2 g 1,max = 0.00124388.

Figure 7 :

 7 Figure 7: Bias Extension test with 1gS-model (ζ = 0) -relative rotation and first gradient shear energy distribution for several imposed displacements: (a-b) u = u lim /4, (c-d) u = u lim /2, (e-f) u = 2u lim /3, (g-h) u = u lim . h = mesh size

Figure 8 :

 8 Figure 8: Bias Extension test with 2gS-model (ζ = 1) -relative rotation and second gradient shear energy distributions for several imposed displacements: (a-b) u = u lim /2, (c-d) u = 2u lim /3, (e-f) u = 0.95u lim , (g-h) u = 0.99u lim . 80 × 240 elements mesh. h = mesh size

Figure 9 :

 9 Figure 9: Bias Extension Test. Top: Relative rotation of the fibres θ 1 + θ 2 along the central axis of the specimen, first and second gradient energy models. Bottom: contour plots of θ 1 + θ 2 for the model 2gS (left) and 2gQ (right): (a) 1gS-and 2gS-models for u 1 = u 2 = (0.99 u lim )/2, (b) 1gQ-and 2gQ-models for u 1 = u 2 = (0.99 u lim )/2, (c) 1gS-and 2gS-models for u 1 = u 2 = 0.99 u lim , (d) 1gQand 2gQ-models for u 1 = u 2 = 0.99 u lim .

Figure 10 :Figure 11 :

 1011 Figure 10: Bias Extension Test. θ 1 along the line indicated in the figure, for increasing values of ζ and for several values of the displacement: (a) u = (0.99 u lim )/4, (b) u = (0.99 u lim )/2, (c) u = 3(0.99 u lim )/4, (d) u = 0.99 u lim . The contour plots refer to the model 2gS (ζ = 1).

Figure 12 :

 12 Figure 12: Schematics of the Shear test.

Ξ

Figure 13 :

 13 Figure 13: Shear Test -First gradient energy model (1gS-model): reference and deformed configurations for several values of the applied displacement.

Ξ

  v = 0.95 v lim .

Figure 14 :

 14 Figure 14: Shear Test -Second gradient energy model (2gS-model): reference and deformed configurations for several values of the applied displacement.

Figure 15 : 2 Figure 16 :

 15216 Figure 15: Shear Test: (a) shear reaction vs. applied displacement, (b) strain energy vs. applied displacement.

Figure 17 :

 17 Figure 17: Shear Test. Top: distribution of θ 1 + θ 2 along the central line of the specimen, for several values of ζ. Bottom: θ 1 + θ 2 for model 2gS, (ζ = 1).: (a) v = (0.99 v lim )/4, (b) v = (0.99 v lim )/2, (c) v = 3(0.99 v lim )/4, (d) v = 0.99 v lim .

Figure 18 :

 18 Figure 18: Shear Test with 1gS-model (ζ = 0) -Relative rotation of the fibres (left) and normalized first gradient shear energy distribution on the specimen (right) for several displacements: (a-b) v = 1/4 v lim , (c-d) v = 1/2 v lim , (e-f) v = 3/4 v lim , (g-h) v = v lim .

Figure 19 :Figure 20 :

 1920 Figure 19: Shear Test with 2gS-model (ζ = 1) -Relative rotation of the fibres (left) and normalized second gradient shear energy distribution on the specimen (right) for several displacements: (a-b) v = 1/4 v lim , (c-d) v = 1/2 v lim , (e-f) v = 3/4 v lim , (g-h) v = 0.95 v lim .32

Figure 21 :

 21 Figure 21: Schematics of the End Rotation Test.

Figure 23 :Figure 24 :

 2324 Figure 23: End Rotation Test -second gradient energy model (2gS-model): reference and deformed configurations for the four values of the applied rotation.

Figure 25 :

 25 Figure 25: End Rotation Test with 1gS-model (ζ = 0) -Relative rotation of the fibres (left) and normalized first gradient shear energy distribution (right) on the specimen for several rotations: (a-b) ϕ o = π/12, (c-d) ϕ o = π/6, (e-f) ϕ o = π/4, (g-h) ϕ o = π/3.

(a) (θ 1 + θ 2 )

 12 max = 0.521131. (b) h 2 g 2,max = 0.00116466. (c) (θ 1 + θ 2 ) max = 1.08234. (d) h 2 g 2,max = 0.00402006. (e) (θ 1 + θ 2 ) max = π/2. (f) h 2 g 2,max = 0.00971574. (g) (θ 1 + θ 2 ) max = π/2.(h) h 2 g 2,max = 0.0995264.

Figure 26 :Figure 27 :

 2627 Figure 26: End Rotation Test with 2gS-model (ζ = 1) -Relative rotation of the fibres (left) and normalized second gradient shear energy distribution (right) on the specimen for several rotations: (a-b) ϕ o = π/12, (c-d) ϕ o = π/6, (e-f) ϕ o = π/4, (g-h) ϕ o = 0.99 × π/3.39

Figure 28 :

 28 Figure 28: End Rotation Test. Top: plots of θ 1 + θ 2 along the central line of the specimen, for various ζ. Bottom: θ 1 +θ 2 for the model 2gS; (a) ϕ o = (0.99 ϕ lim )/4, (b) ϕ o = (0.99 ϕ lim )/2, (c) ϕ o = 3(0.99 ϕ lim )/4, (d) ϕ o = 0.99 ϕ lim .

Figure 29 :

 29 Figure 29: Cantilever test: (a) Set-up of the problem, (b) Force versus displacement, for the 1gS-and 2gS-models.

  u F /W = 1.6814, point A.

  u F /W = 2.1346, point B.

  u F /W = 2.5958, point D.

Figure 30 :

 30 Figure 30: Cantilever test, 1gS-model: reference and deformed configuration.

  u F /W = 0.7035, point A.

  u F /W = 1.3557, point B.

  u F /W = 2.3553, point D.

Figure 31 :

 31 Figure 31: Cantilever test, 2gS-model: reference and deformed configuration.