J Yvonnet
email: julien.yvonnet@univ-paris-est.fr

N Auffray

V Monchiet

Computational second-order homogenization of materials with effective anisotropic strain-gradient behavior

Keywords: Computational homogenization, Strain-gradient, Second-order homogenization, Finite Element

Computational second-order homogenization of materials with effective anisotropic strain-gradient behavior

Introduction

Architectured materials are multiscaled materials for which the inner structure can be various, ranging from randomly distributed phases to perfectly organized architectures. In the last case, we are in the realm of periodic materials and the mesoscopic scale is composed of a basic cell repeated according to a regular lattice. Intermediate configurations are possible, for instance in biological architectured materials, in which almost periodic structures or perturbed lattices are often encountered [START_REF] Barthelat | Architectured materials in engineering and biology: fabrication, structure, mechanics and performance[END_REF]. The great influence of the mesostructure on the overall properties of materials is known in condensed matter physics, at least from the end of the 19th century, where it was found to drive optical and electric properties of crystals [START_REF] Curie | Sur la symétrie dans les phénomènes physiques, symétrie d'un champ électrique et d'un champ magnétique[END_REF]. At the scale of mechanical engineering, interest in architectured materials lies in the possibility to obtain overall behaviors which are not found in nature. Non-standard elastodynamics properties such as, for instance, the cloaking effect, are at the very beginning of the studies of what is presently known as metamaterials [START_REF] Milton | On cloaking for elasticity and physical equations with a transformation invariant form[END_REF][START_REF] Norris | Elastic cloaking theory[END_REF][START_REF] Liu | Chiral effect in plane isotropic micropolar elasticity and its application to chiral lattices[END_REF]. But interest is not confined to dynamics, since non-standard static effects such as large auxetic behavior have been investigated [START_REF] Reis | Equivalent mechanical properties of auxetic lattices from discrete homogenization[END_REF][START_REF] Dirrenberger | Effective elastic properties of auxetic microstructures: anisotropy and structural applications[END_REF][START_REF] Saxena | Three decades of auxetics researchmaterials with negative poisson's ratio: a review[END_REF], or pure strain-gradient designed structures [START_REF] Seppecher | Linear elastic trusses leading to continua with exotic mechanical interactions[END_REF][START_REF] Dell'isola | Pantographic metamaterials: an example of mathematically driven design and of its technological challenges[END_REF][START_REF] Abdoul-Anziz | Strain gradient and generalized continua obtained by homogenizing frame lattices[END_REF]. The theoretical formulation of continuum mechanics with microstructure has been investigated since the mid-60's with the seminal works of Mindlin, Toupin, and others [START_REF] Toupin | Elastic materials with couple-stresses[END_REF][START_REF] Eringen | Nonlinear theory of simple micro-elastic solids-i[END_REF][START_REF] Mindlin | Micro-structure in linear elasticity[END_REF]. All these works demonstrated that the classical formulation of elasticity (Cauchy elasticity) is not rich enough to model architecture effects in a continuum fashion, and that generalized models have to be considered instead. However, despite their intrinsic interests and until recently, those results were confined to essentially theoretical speculations due to, at least, two major problems:

(1) How to practically design architectured materials? (2) How to determine the very large amount of parameters required by these models?

With regard to the first point, recent developments in additive manufacturing processes have led to a paradigm shift and, in recent years, new direct manufacturing methods have been developed. Today, 3D printing techniques have allowed fabricating lattices [START_REF] Bückmann | Tailored 3D mechanical metamaterials made by dip-in direct-laser-writing optical lithography[END_REF] or bi-materials [START_REF] Qu | Micro-structured two-component 3D metamaterials with negative thermal-expansion coefficient from positive constituents[END_REF] with controlled microstructures over a wide range of scales and materials, from nano to micro lattices [START_REF] Vaezi | A review on 3D micro-additive manufacturing technologies[END_REF] and made of metallic, polymer or ceramics materials, among many others.

However, due to the complex internal geometry of architectured materials, Direct Numerical Simulations (DNS, i.e. when the geometry of the inner structure is explicitly meshed in simulations) are often limited to few unit cells, or few Representative Volume Elements (RVEs). In other words, the full-field computation of a structure made out of an architectured material can be prohibitively expensive, especially in the design phase in which the structure is optimized. To solve these technical difficulties, a possible way is to substitute the actual material by a homogeneous one having equivalent effective properties. However, due to the meso-structure coarseness, the usual scale separation assumption may no more hold and standard homogenization techniques break down. Further, to correctly catch the behavior of the original heterogeneous material, the effective homogeneous medium cannot be Cauchy elastic and should incorporate higher-order effects in its formulation. At this level the questions are: (a) which generalized continuum to consider among the multiple possibilities (Cosserat, Koiter, Strain-gradient, Micromorphic,...) (b) how to compute the coefficients of the retained model from the knowledge of the microstructure? It is where we meet the second problem previously evoked.

In the present paper, the model retained for the overall description of the heterogeneous continuum is a strain gradient elastic continuum as described by Mindlin [START_REF] Mindlin | Micro-structure in linear elasticity[END_REF][START_REF] Mindlin | On first strain-gradient theories in linear elasticity[END_REF]. The reasons for such as choice are three-fold: (a) SGE model can be naturally constructed by asymptotic analysis [START_REF] Boutin | Microstructural effects in elastic composites[END_REF][START_REF] Bacigalupo | Second-order computational homogenization of heterogeneous materials with periodic microstructure[END_REF][START_REF] Tran | A micromechanics-based approach for the derivation of constitutive elastic coefficients of strain-gradient media[END_REF][START_REF] Bacigalupo | Second-gradient homogenized model for wave propagation in heterogeneous periodic media[END_REF]; (b) in statics, SGE elasticity is rich enough to describe emergent properties without involving as many parameters as the micromorphic elasticity [START_REF] Auffray | A complete description of bidimensional anisotropic strain-gradient elasticity[END_REF]; (c) we have a good understanding of the anisotropic features of this model [START_REF] Auffray | Derivation of anisotropic matrix for bidimensional strain-gradient elasticity behavior[END_REF][START_REF] Auffray | A complete description of bidimensional anisotropic strain-gradient elasticity[END_REF]. The question now lies in the determination of the parameters of the model.

With regard to this modelling issue, there is a large interest in the literature for strain-gradient homogenization of heterogeneous structures [START_REF] Ben-Amoz | A dynamic theory for composite materials[END_REF][START_REF] Forest | Cosserat overall modelling of heterogeneous materials[END_REF][START_REF] Forest | Mechanics of generalized continua: construction by homogenizaton[END_REF][START_REF] Forest | Asymptotic analysis of heterogeneous cosserat media[END_REF][START_REF] Trinh | Evaluation of generalized continuum substitution models for heterogeneous materials[END_REF][START_REF] Tran | A micromechanics-based approach for the derivation of constitutive elastic coefficients of strain-gradient media[END_REF][START_REF] Rahali | Homogenization à la piola produces second gradient continuum models for linear pantographic lattices[END_REF]. For instance, in [START_REF] Bacigalupo | Second-order computational homogenization of heterogeneous materials with periodic microstructure[END_REF], Bacigalupo and Gambarotta provided a computational homogenization procedure for strain-gradient media with numerical solutions. In [START_REF] Bacigalupo | Homogenization of periodic hexa-and tetrachiral cellular solids[END_REF], the same authors applied micropolar and strain-gradient homogenization schemes to auxetic materials (hexachiral and tetrachiral honeycombs) and obtained analytical expressions for effective higher-order coefficients, as well as numerical values based on FEM calculations. An anisotropic effective straingradient behavior has been constructed based on RVE of fibrous materials in [START_REF] Goda | Construction of first and second order grade anisotropic continuum media for 3D porous and textile composite structures[END_REF] using heuristic tests as boundary conditions. In [START_REF] Yuan | A micromechanical approach of nonlocal modeling for media with periodic microstructures[END_REF], a micromechanical model for macro behavior when prescribing a polynomial field based on asymptotic analysis was introduced and a critical review of the different approaches proposed in the literature was provided. Following [START_REF] Boutin | Microstructural effects in elastic composites[END_REF][START_REF] Smyshlyaev | On rigorous derivation of strain gradient effects in the overall behaviour of periodic heterogeneous media[END_REF], Tran et al. [START_REF] Tran | A micromechanics-based approach for the derivation of constitutive elastic coefficients of strain-gradient media[END_REF] derived closed-form estimates for stratified composites using asymptotic expansion methods and a Fast Fourier Transform approach to compute the effective coefficients of the strain-gradient elastic model. In [START_REF] Goda | Construction of first and second order grade anisotropic continuum media for 3D porous and textile composite structures[END_REF][START_REF] Berkache | Construction of second gradient continuum models for random fibrous networks and analysis of size effects[END_REF], second-order coefficients have been obtained in a linear context using heuristic boundary conditions while, in [START_REF] Kaczmarczyk | Scale transition and enforcement of rve boundary conditions in second-order computational homogenization[END_REF], the issue of enforcing boundary conditions in secondorder computational homogenization has been addressed using constrained minimization. In [START_REF] Li | A numerical approach for the establishment of strain gradient constitutive relations in periodic heterogeneous materials[END_REF], Li and Zhang combined asymptotic expansion and FFT and discussed gradient effects. The approach is not restricted linear behaviors, and second-order computational homogenization has been extended to the nonlinear context through the FE 2 method in [START_REF] Kouznetsova | Multi-scale constitutive modeling of heterogeneous materials with gradient enhanced computational homogenization scheme[END_REF][START_REF] Kouznetsova | Multi-scale second order computational homogenization of multi-phase materials: a nested finite element solution strategy[END_REF][START_REF] Feyel | A multilevel finite element method (FE 2 ) to describe the response of highly non-linear structures using generalized continua[END_REF]. In [START_REF] Lesičar | A second-order two-scale homogenization procedure using c 1 macrolevel discretization[END_REF], a FE 2 extension of second-order homogenization with C 1 macro discretization has been introduced.

Amongst the approaches used in the above examples, a very common strategy to compute higher-order moduli is the use of Quadratic Boundary Conditions (QBC) over the RVE [START_REF] Gologanu | Recent extensions of gurson's model for porous ductile metals[END_REF][START_REF] Forest | Mechanics of generalized continua: construction by homogenizaton[END_REF][START_REF] Forest | Asymptotic analysis of heterogeneous cosserat media[END_REF]. This method is popular due to its inherent simplicity, and to the fact that it is a straightforward extension of the classical Kinematic Uniform Boundary Conditions (KUBC) used in classical homogenization [START_REF] Huet | Application of variational concepts to size effects in elastic heterogeneous bodies[END_REF][START_REF] Kanit | Determination of the size of the representative volume element for random composites: statistical and numerical approach[END_REF][START_REF] Yvonnet | Computational Homogenization of Heterogeneous Materials with Finite Elements[END_REF]. This approach has been used in numerous situations to determine strain-gradient effective parameters as in [START_REF] Geers | Gradient-enhanced computational homogenization for the micro-macro scale transition[END_REF][START_REF] Kouznetsova | Multi-scale second order computational homogenization of multi-phase materials: a nested finite element solution strategy[END_REF] and further extended to micromorphic effective media [START_REF] Forest | Cosserat overall modelling of heterogeneous materials[END_REF][START_REF] Forest | Asymptotic analysis of heterogeneous cosserat media[END_REF][START_REF] Bouyge | A micromechanically based couple-stress model of an elastic orthotropic two-phase composite[END_REF][START_REF] Jänicke | Two-scale modelling of micromorphic continua[END_REF][START_REF] Forest | Generalized continua and non-homogeneous boundary conditions in homogenisation methods[END_REF]. Despite its appealing simplicity, as shown in [START_REF] Forest | Generalized continua and non-homogeneous boundary conditions in homogenisation methods[END_REF][START_REF] Yuan | A micromechanical approach of nonlocal modeling for media with periodic microstructures[END_REF][START_REF] Barboura | Establishment of strain gradient constitutive relations by using asymptotic analysis and the finite element method for complex periodic microstructures[END_REF], the method has two major flaws: (a) the effective strain-gradient properties remain non zero either when the scale are separated or when the material is homogeneous [START_REF] Forest | Generalized continua and non-homogeneous boundary conditions in homogenisation methods[END_REF] which induces a contradiction with the local Cauchy behavior of the phases; (b) the effective straingradient properties do not converge with respect to the RVE size [START_REF] Li | A numerical approach for the establishment of strain gradient constitutive relations in periodic heterogeneous materials[END_REF][START_REF] Barboura | Establishment of strain gradient constitutive relations by using asymptotic analysis and the finite element method for complex periodic microstructures[END_REF], which is also questionable. In [START_REF] Li | Establishment of strain gradient constitutive relations by homogenization[END_REF][START_REF] Barboura | Establishment of strain gradient constitutive relations by using asymptotic analysis and the finite element method for complex periodic microstructures[END_REF], Barboura and Li proposed a correction to the asymptotic analysis combined with the finite element method to remove these spurious effects.

In the present paper, we provide a computational homogenization method for evaluating higher-order tensors of a general effective anisotropic straingradient (Mindlin) model using Finite Element Calculations on RVEs based on the superposition principle. As compared to previous methods, body forces are added to the Quadratic Boundary Conditions and a correction of the localization tensor in a straightforward application of the superposition principle is introduced to remove the above-mentioned issues. We show that the resulting method is consistent with results obtained by asymptotic expansion techniques [START_REF] Tran | A micromechanics-based approach for the derivation of constitutive elastic coefficients of strain-gradient media[END_REF]. Finally, a C 1 -FEM for the effective anisotropic Mindlin model is provided, which can serve to validate the homogenized model.

The outline of the paper is as follows. After notational preliminaries in section 2, the localization problem with quadratic boundary conditions and additional body forces is introduced in section 3. In section 4 the effective higher-order tensors are derived. In section 5, a full C 1 -FEM to calculate the homogeneous structure with the computed anisotropic strain-gradient behavior is provided. Finally, numerical examples are presented in section 5.3 to validate and illustrate the approach, as well as presenting some current limitations in the case of soft inclusions.

Preliminary notations

Throughout this paper, the physical space is modelled on the Euclidean space E d with R d its associated vector space. Once an arbitrary reference point chosen, those spaces can be associated and P = {e k } 1≤k≤d will denote an orthonormal basis of R d .

Vectors and second order tensors, as well as matrices, are denoted by bold letters A. Third-order tensors are denoted by calligraphic uppercase letters G, fourth-order, fifth-order and sixth-order tensors are denoted by double case letters A. The simple, double and third-order contractions are written •, :, . . . respectively. In components, with respect to P, these notations correspond to

a • b = a i b i , A : B = A ij B ij , A . . .B = A ijk B ijk , (1)
where the Einstein summation on repeated indices is used. The order on which components are summed has been chosen in such a way that the full contraction between same order tensors is a scalar product.

The outer product between two vectors is defined by (a ⊗ b) ij = a i b j . The gradient operator is denoted by ∇(.) and the divergence operator by ∇•(.). The classical second order identity tensor is denoted by I. When needed secondorder and fourth-order tensors will be represented as vectors and matrices using the classical Voigt's convention.

To characterize the symmetry of the architectured materials considered in this paper, some definitions from group theory will be used. The following matrix groups are considered:

• Due to its practical importance, let us detail the case d = 2. As a matrix group, O(2) can be generated by:

R(θ) = ⎛ ⎜ ⎝ cos θ -sin θ sin θ cos θ ⎞ ⎟ ⎠ , 0 ≤ θ < 2π, and P(e 2 ) = ⎛ ⎜ ⎝ 1 0 0 -1 ⎞ ⎟ ⎠ ,
in which R(θ) is a rotation by an angle θ and P(n) is the reflection across the line normal to n. SO(2) corresponds to the group of rotations generated by R(θ). The following finite subgroups of O(2) will be used:

• 1, the identity group;

• Z k , the cyclic group with k elements generated by R(2π/k). It is the symmetry group of chiral figures; • D k , the dihedral group with 2k elements generated by R(2π/k) and P(e 2 ).

It is the symmetry group of regular polygons. The group D 1 will also be denoted by Z π 2 , and is the symmetry group of a figure solely invariant by P(e 2 ).

Let u be the displacement vector and x a material coordinate, since the present paper is devoted to linear strain-gradient elasticity we introduce the linear strain tensor ε and ∇ε its gradient, defined by, respectively:

ε ij = 1 2 ∂u i ∂x j + ∂u j ∂x i , ∇ε ijk = 1 2 ∂ 2 u i ∂x j ∂x k + ∂ 2 u j ∂x i ∂x k . (2)
The second gradient of displacement is defined as

G ijk = ∂ 2 u i ∂x j ∂x k . (3)
The two tensors ∇ε and G are related by (see Appendix A):

G ijk = ∇ε ijk + ∇ε ikj -∇ε jki . (4)
In this work, we adopt the formalism of Mindlin strain-gradient elasticity as a function of the strain-gradient ∇ε, this is the type II formulation as described in [START_REF] Mindlin | Micro-structure in linear elasticity[END_REF][START_REF] Mindlin | On first strain-gradient theories in linear elasticity[END_REF]. The theory can also be formulated in terms of the second derivative of displacement G (type I formulation), the correspondence between the two formalisms can be found e.g. in [START_REF] Mindlin | On first strain-gradient theories in linear elasticity[END_REF][START_REF] Auffray | A complete description of bidimensional anisotropic strain-gradient elasticity[END_REF].

• (•) indicates a macroscopic quantity;

• ⟨(•)⟩ V denotes the volume averaging of the quantity (•) over the domain V :

⟨(•)⟩ V = 1 |V | V (•)dV.
3 Localization problem and boundary conditions

Definition of the localization problem

We consider a structure defined in a domain Ω ⊂ R 2 , which characteristic size is L, constituted of a finite number of unit cells of elementary domain Y.

The material filling Ω is hence considered as a periodic medium characterized by a RVE defined in a domain ω ⊂ R 2 whose boundary is denoted by ∂ω (see Fig. 1 (c)). It has to be noted that ω can be larger than Y, hence the notions of RVE and of unit cell have to be distinguished. Note that in this work we distinguish ω from Y because non periodic boundary conditions are used on the boundary of the RVE. The RVE, whose characteristic size is denoted by ℓ, is supposed constituted of N different Cauchy linear elastic phases characterized by their elastic tensors C i , i = 1, 2, ..., N. This RVE ω is assumed to be subjected to (ε, ∇ε) a couple of homogeneous strain and strain gradient tensors, which are prescribed through boundary conditions and body forces specified in the following. The scale separation ratio is defined as ϵ = ℓ/L and is not vanishingly small. In that situation, the assumption of scale separation is no longer valid and the effective medium may include strain-gradient effects [START_REF] Rosi | On the validity range of strain-gradient elasticity: a mixed static-dynamic identification procedure[END_REF][START_REF] Poncelet | An experimental evidence of the failure of cauchy elasticity for the overall modeling of a noncentro-symmetric lattice under static loading[END_REF] (see section 4).

To extend the classical homogenization procedure based on Kinematic Uniform Boundary Conditions (KUBC) [START_REF] Huet | Application of variational concepts to size effects in elastic heterogeneous bodies[END_REF][START_REF] Kanit | Determination of the size of the representative volume element for random composites: statistical and numerical approach[END_REF], Quadratic Boundary Conditions (QBC) have been introduced in several works (see e.g. [START_REF] Gologanu | Recent extensions of gurson's model for porous ductile metals[END_REF][START_REF] Forest | Mechanics of generalized continua: construction by homogenizaton[END_REF][START_REF] Forest | Asymptotic analysis of heterogeneous cosserat media[END_REF][START_REF] Auffray | Strain gradient elastic homogenization of bidimensional cellular media[END_REF]) to prescribe an effective strain-gradient over the RVE (see Appendix B):

u(x) = ε • x + 1 2 G : x ⊗ x on ∂ω, ( 5
)
where G is the macroscopic second gradient of displacement, which depends on ∇ε according to (4). However, one obvious issue with such condition arises when considering a homogeneous RVE characterized by a unique elastic tensor C 1 .

In that case, and for ε = 0, it is expected that the local strain solution within the RVE should be equal to:

ε(x) = ∇ε • x ∀x ∈ ω. (6)
However generally ( 6) is not a statically admissible solution for boundary conditions (5) since:

∇ • C 1 : ∇ε x = (C 1 : ∇ε) : I ̸ = 0. ( 7
)
The inequality [START_REF] Bacigalupo | Second-order computational homogenization of heterogeneous materials with periodic microstructure[END_REF] holds because in the present work ∇ε can be chosen arbitrarily. Therefore, as observed in [START_REF] Forest | Generalized continua and non-homogeneous boundary conditions in homogenisation methods[END_REF][START_REF] Li | A numerical approach for the establishment of strain gradient constitutive relations in periodic heterogeneous materials[END_REF], fluctuations remain even when the local continuum is homogeneous, leading to persistent non-physical gradient effects. Indeed, when the local medium is Cauchy homogeneous, there is no dependence on an internal length and the overall medium cannot be of generalized type. To cure this problem, and following the analysis conducted in [START_REF] Monchiet | Strain-gradient homogenization: a bridge between asymptotic expansion and quadratic boundary condition methods[END_REF], we propose to prescribe body forces in addition to QBC [START_REF] Auffray | Complete symmetry classification and compact matrix representations for 3D strain gradient elasticity[END_REF] to enforce a constant strain-gradient within the RVE when the material is homogeneous.

The new localization problem involves solving the equilibrium equation:

∇ • σ (u(x)) = f (∇ε) ∀x ∈ ω, (8)
where

f (∇ε) = ∇ • C 0 (x) : (∇ε • x) . (9)
In the definition of f , C 0 (x) is an arbitrary elastic tensor field that has to be specified. At this point, and without loss of generality, we assume a twophase composite whose elastic properties are described by C 1 and C 2 , and in which the phase 1 has the highest volume fraction. The RVE is elastically homogeneous if either (a) the volume fraction of phase 2 goes to zero, i.e. f 1 → 1, or (b) if the contrast between phase properties goes to one, i.e. ∥C 2 ∥ → ∥C 1 ∥. For each of these two conditions, the tensor C 0 (x) should satisfy:

C 0 (x) → C 1 if ⎧ ⎪ ⎨ ⎪ ⎩ f 1 → 1, or ∥C 2 ∥ → ∥C 1 ∥ (10)
to satisfy [START_REF] Bacigalupo | Homogenization of periodic hexa-and tetrachiral cellular solids[END_REF].

Several choices are possible to respect condition [START_REF] Barboura | Establishment of strain gradient constitutive relations by using asymptotic analysis and the finite element method for complex periodic microstructures[END_REF], such as (among others):

Pointwise body force correction:

C 0 (x) = C(x), Effective body force correction: C 0 (x) = C,
while the standard one is

Null body force (Standard QBC): C 0 (x) = O,
where C is the effective classical elastic tensor of the RVE. In the present work, we consider the Effective Body Force Correction, but discuss the cases C 0 (x) = C(x) (Pointwise Body Force Correction) as well as the case C 0 (x) = 0 (Null body force (Standard QBC)) for comparison. For this last choice, the condition (10) is obviously not satisfied.

To summarize, the localization problem to be solved on the RVE is the following one:

Given ε and ∇ε, find u(x) such that:

∇ • σ (u(x)) = f (∇ε), ∀x ∈ ω, ( 11
) with σ (x) = C(x) : ε (u(x)) , (12)
and subjected to the boundary conditions:

u(x) = εx + 1 2 G : x ⊗ x ∀x ∈ ∂ω, (13)
and the source term:

f (∇ε) = ∇ • C 0 (x) : (∇ε • x) . (14)
The problem defined by the equations ( 11)-( 13) can be solved by classical finite elements as described below.

FEM Formulation

Denoting by [ε(x)] the vector form associated with ε(x) and [σ(x)] the vector form associated with σ(x), the weak form associated with ( 11) is given by:

ω [σ] • [ε(δu)]dΩ = ω C(x)[η] • [ε(δu)]dΩ, ( 15
)
where C is the matrix form associated to C and by setting:

[η] = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ∇ε 111 x 1 + ∇ε 112 x 2 ∇ε 221 x 1 + ∇ε 222 x 2 ∇ε 121 x 1 + ∇ε 122 x 2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ . ( 16
)
Using classical FEM discretization, we obtain a linear system in the form:

Ku = F, (17)
with

K = ω B(x) T C(x)B(x)dΩ, and F = ω B T (x)C(x)[η]dΩ, (18)
where B(x) are the shape function derivatives in the elements of the RVE.

4 Calculation of effective higher-order tensors

The Strain-Gradient constitutive law

In its Type II formulation, the energy of a macroscopic linear strain-gradient elastic medium has the following expression [START_REF] Mindlin | Micro-structure in linear elasticity[END_REF][START_REF] Mindlin | On first strain-gradient theories in linear elasticity[END_REF],:

E = 1 2 ε : C : ε + ε : M . . . ∇ε + 1 2 ∇ε . . . G . . . ∇ε, (19)
where C denotes the effective elastic tensor, M is a fifth-order tensor coupling first and second order strains, and G is the sixth-order strain-gradient elasticity tensor. These tensors satisfy the following index permutation symmetry:

C (ij) (kl) ; M (ij)(kl)m ; G (ij)k (lm)n (20)
where the notation (..) stands for the minor symmetries, whereas .. stands for the major one. In the case where the microstructure of a material is centrosymmetric, the fifth-order elastic stiffness tensor M vanishes. This assumption is often made in order to uncoupled the equations and hence to simplify the problem. However, as demonstrated for example by [START_REF] Boutin | Microstructural effects in elastic composites[END_REF] from an asymptotic analysis, the moduli of M are of order ϵ while those of G are of order ϵ 2 and therefore, when they are present, the influence of M can be dominant over G.

In statics, some recent experiments on a non-centrosymmetric architectured beam [START_REF] Poncelet | An experimental evidence of the failure of cauchy elasticity for the overall modeling of a noncentro-symmetric lattice under static loading[END_REF] have evidenced the necessity of involving M in modelling its overall behavior. In dynamics, the coupling described by M induces changes in wave polarization and is associated to the acoustic activity of crystals and to the so-called gyrotropic effects [START_REF] Toupin | Elastic materials with couple-stresses[END_REF][START_REF] Portigal | Acoustical activity and other first-order spatial dispersion effects in crystals[END_REF][START_REF] Maranganti | A novel atomistic approach to determine straingradient elasticity constants: Tabulation and comparison for various metals, semiconductors, silica, polymers and the (ir) relevance for nanotechnologies[END_REF]. In this work, we consider arbitrarily anisotropic RVEs, with M possibly different from zero. Then the second-order stress tensor σ and the third-order hyperstress tensor S are defined, respec-tively, by:

σ = ∂E ∂ε = C : ε + M . . . ∇ε, (21)
S = ∂E ∂∇ε = M T : ε + G . . . ∇ε. ( 22
)
with

M T ijklm = M lmijk .
According to the symmetries of the RVE, this constitutive law can belong, in 2D, to 14 non equivalent symmetry classes. Their characteristics are reported together with their number of independent components in Table 1:

Name Oblique Rectangular Digonal Orthotropic Trichiral Trigonal Tetrachiral [G L ] [1] [Z π 2 ] [Z 2 ] [D 2 ] [Z 3 ] [D 3 ] [Z 4 ] # indep (L) 45 (
[G L ] [D 4 ] [Z 5 ] [D 5 ] [Z 6 ] [D 6 ] [SO(2)] [O(2)] # indep (L) 9 9 (8) 7 9 (8) 7 7
6 Table 1 Names, Sets of subgroups [G L ] and numbers of independent components # indep (L) for the 14 symmetry classes of L. The in-parenthesis number indicates the minimal number of components of the law in an appropriate basis.

in which L = (C, M, G). Details concerning this classification can be found in the following references [START_REF] Auffray | Derivation of anisotropic matrix for bidimensional strain-gradient elasticity behavior[END_REF][START_REF] Auffray | A complete description of bidimensional anisotropic strain-gradient elasticity[END_REF]. The same analysis has also been conducting in 3D [START_REF] Auffray | Matrix representations for 3D straingradient elasticity[END_REF][START_REF] Auffray | Complete symmetry classification and compact matrix representations for 3D strain gradient elasticity[END_REF], but results are not reported here for the sake of conciseness.

Definition of the localization tensors

The local problem defined by the equations ( 11)-( 13) is linear, hence using the superposition principle, the local strain field ε(x) can be expressed as:

ε(x) = A 0 (x) : ε + A 1 (x) . . . ∇ε, (23)
where A 0 ijkl (x) is the strain solution ε ij (x) obtained by solving the problem (11)-( 13) for:

ε = 1 2 (e k ⊗ e l + e l ⊗ e k ) , ∇ε = 0 ( 24
)
and A1 ijklm (x) is the strain solution ε ij (x) obtained by solving the problem (11)-( 13) for:

ε = 0, ∇ε = 1 2 (e k ⊗ e l + e l ⊗ e k ) ⊗ e m . (25)
Even when prescribing appropriate body forces as described in the previous section, the localization tensors as they were introduced do not prevent the non-physical strain-gradient effects which show up in the case of a local homogeneous problem. To understand why, consider the case of a homogeneous RVE with elastic tensor C 1 and apply the homogenization scheme with prescribed body forces to balance the QBC. In such a situation the localization tensors are 1 :

A 0 (x) = I ∀x ∈ Ω, A 1 (x) = I ⊗ x ∀x ∈ Ω, ( 26
)
where I is the fourth-order symmetric identity tensor defined by (I

) ijkl = 1 2 (δ ik δ jl + δ il δ jk ).
In that case, and using the standard local Hooke's law, the expressions of the strain and the stress fields are readily obtained as:

ε(x) = I : ε + (I ⊗ x) : ∇ε , σ(x) = C 1 : I : ε + (I ⊗ x) . . . ∇ε . (27)
The effective strain density energy associated with the RVE is expressed by:

E = 1 2 ⟨σ(x) : ε(x)⟩ , (28)
and by introducing ( 27) in [START_REF] Germain | La méthode des puissances virtuelles en mécanique des milieux continus, premiere partie: théorie du second gradient[END_REF], we obtain:

E = 1 2 ε : I : C 1 : I : ε + ε : I : C 1 : I ⊗ x . . . ∇ε + 1 2 ∇ε . . . x ⊗ I : C 1 : I ⊗ x . . . ∇ε. (29)
Considering energy equivalence between micro and macro scales (Hill-Mandel lemma's extension [START_REF] Gologanu | Recent extensions of gurson's model for porous ductile metals[END_REF][START_REF] Auffray | A complete description of bidimensional anisotropic strain-gradient elasticity[END_REF]), and comparing ( 29) with [START_REF] Dirrenberger | Effective elastic properties of auxetic microstructures: anisotropy and structural applications[END_REF] we deduce that for a homogeneous RVE, M and G reduce to:

M ijklm = C 1 ijkl ⟨x m ⟩ , G ijklmp = C 1 jklm ⟨x i x p ⟩ (30)
which are expected to be zero in the present situation. However, M can only be zero if the RVE is of centrosymmetric shape and if the origin of the local frame is located at the center of the RVE. For its part, G, is always different from zero. In both cases, the values of M and G depend on the choice of the origin related to the frame of the RVE, which is not satisfying. In conclusion, the localization tensors as defined by the relation [START_REF] Forest | Mechanics of generalized continua: construction by homogenizaton[END_REF] are not appropriate to remove spurious strain-gradient effects in the case of homogeneous RVEs. To correct this, and following the analysis developed in [START_REF] Monchiet | Strain-gradient homogenization: a bridge between asymptotic expansion and quadratic boundary condition methods[END_REF], we propose to use the modified relation:

ε(x) = A 0 (x) : ε + Ã1 (x) . . . ∇ε (31)
in which

Ã1 (x) = A 1 (x) -A 0 (x) ⊗ x. (32)
The modified definition for Ã1 amounts to remove from A 1 the influence of linear part of the elastic field. By doing so, the second term of (31) only contains purely strain-gradient contribution. Using the same procedure as above, we recover the stress as:

σ(x) = C 1 (x) : A 0 (x) : ε + Ã1 (x) . . . ∇ε . (33)
Introducing ( 31) and ( 33) in ( 28), we obtain:

E = 1 2 ε : A 0 (x) T : C(x) : A 0 (x) : ε +ε : A 0 (x) T : C(x) : Ã1 (x) . . . ∇ε + 1 2 ∇ε . . . Ã1 (x) T : C(x) : Ã1 (x) . . . ∇ε. (34)
Considering energy equivalence between micro and macro scales, and comparing [START_REF] Kanit | Determination of the size of the representative volume element for random composites: statistical and numerical approach[END_REF] with [START_REF] Dirrenberger | Effective elastic properties of auxetic microstructures: anisotropy and structural applications[END_REF] the following expressions for the effective tensors are obtained as:

C = A 0 (x) T : C(x) : A 0 (x) , (35)
M = A 0 (x) T : C(x) : Ã1 (x) , (36)
G = Ã1 (x) T : C(x) : Ã1 (x) . (37)
Using these new definitions for the effective tensors, it can easily be verified that, for a homogeneous RVE, the following results are verified:

M = O, G = O. (38)
It is worth noting that the present model is consistent with the results obtained by asymptotic expansion techniques [START_REF] Tran | A micromechanics-based approach for the derivation of constitutive elastic coefficients of strain-gradient media[END_REF], where the same expressions for the effective higher-order tensors [START_REF] Kouznetsova | Multi-scale second order computational homogenization of multi-phase materials: a nested finite element solution strategy[END_REF], [START_REF] Lesičar | A second-order two-scale homogenization procedure using c 1 macrolevel discretization[END_REF] are obtained, and where the localization tensors A 0 (x) and A 1 (x) are also solutions of an equilibrium problem with body forces involving the effective elastic tensor C.

The local equilibrium in the case C 0 (x) = C has been obtained by Boutin [START_REF] Boutin | Microstructural effects in elastic composites[END_REF] and also considered in Tran et al. [START_REF] Tran | A micromechanics-based approach for the derivation of constitutive elastic coefficients of strain-gradient media[END_REF] and Monchiet et al. [START_REF] Monchiet | Strain-gradient homogenization: a bridge between asymptotic expansion and quadratic boundary condition methods[END_REF]. Note that the introduction of the body force in the asymptotic homogenization approach is required to comply with the anti-periodicity of the traction on two opposite sides of the unit cell.

The definition [START_REF] Jänicke | Two-scale modelling of micromorphic continua[END_REF] for the correction term to the fundamental elastic solution has been recently derived in Monchiet et al. [START_REF] Monchiet | Strain-gradient homogenization: a bridge between asymptotic expansion and quadratic boundary condition methods[END_REF]. In this paper, a connection has been established between higher order homogenization approaches based on asymptotic expansions and Quadratic Boundary Condition on the boundary. The principle of the approach consists in replacing the periodicity condition by the condition u = 0 on the boundary for the fluctuation at any order of the expansion series. By doing so, it has been found that the first higher order term of the strain is given by:

ε 1 = ε 1 - ∂ε 0 ∂x k y k ( 39
)
where ε 1 is the local strain in the RVE with prescribed macroscopic gradient of strain and ε 0 is the local strain due to prescribed macroscopic strain. It must be emphasis that, by computing the difference between ε 1 and ∂ε 0 ∂x k y k , we eliminate the fundamental elastic response and we only keep the flexural effects. Replacing ε 1 by A 1 (x) . . . ∇ε and ε 0 by A 0 (x) . . . ε one recovers Eq. ( 32).

5 Macroscopic strain-gradient problem and C 1 Finite Element discretization

In this section, we define the elastostatic problem associated with the homogeneous anisotropic Mindlin strain-gradient elasticity model together with its Finite Element discretization. These results will be used in section 6 for comparing fully-meshed computations with their strain-gradient homogenized counterparts.

From strong to weak form

Let Ω ⊂ R 2 a domain associated with the homogeneous strain-gradient medium and ∂Ω its boundary (see Fig. 1 (c)), the momentum balance for the straingradient medium is given by (see e.g. [START_REF] Mindlin | Micro-structure in linear elasticity[END_REF][START_REF] Germain | La méthode des puissances virtuelles en mécanique des milieux continus, premiere partie: théorie du second gradient[END_REF]):

∇ • σ -∇ • ∇ • S = 0 in Ω, (40)
with boundary conditions:

u = u * on ∂Ω u , (41)
σ • n -∇ • S • n + ∇ S • S • n = F * on ∂Ω F , (42)
S : (n ⊗ n) = 0 on ∂Ω, (43)
S : (n 1 ⊗ n 2 ) = 0 on ∂∂Ω, (44)
(see e.g. [START_REF] Liu | An energy formulation of continuum magneto-electro-elasticity with applications[END_REF] for a justification), and where u * and F * denote prescribed displacements and tractions, respectively, ∂Ω u and ∂Ω F are the Dirichlet and Neumann parts of the boundary ∂Ω, and ∇ S • indicates the surface divergence operator. As we assume no prescribed hyperstress over the structure, there is nor additional condition involving a jump in the normal along the surface [START_REF] Mindlin | Micro-structure in linear elasticity[END_REF] neither condition depending on the curvature of ∂Ω [START_REF] Maranganti | A novel atomistic approach to determine straingradient elasticity constants: Tabulation and comparison for various metals, semiconductors, silica, polymers and the (ir) relevance for nanotechnologies[END_REF] .

The weak form associated with problem (40) is given by:

Ω σ : ε(δu) + S . . . ∇ε(δu)dΩ = ∂Ω F F * • δudΓ. (45)
Using ( 21) and ( 22) we obtain:

Ω C : ε(u) + M . . . ∇ε(u) : ε(δu)dΩ + Ω M T : ε(u) + G . . . ∇ε . . . ∇ε(δu)dΩ = ∂Ω F F * • δudΓ. (46)

FEM Discretization of the strain-gradient problem

In the following, we present the developments associated with the Finite Element discretization of the above problem. For the sake of simplicity, the developments are restricted to 2D. A mesh of triangular elements is constructed over the domain Ω, containing N e elements. It is worth noting that Eq. ( 40) is a fourth-order system of partial differential equations with respect to u, which require C 1 continuity for the components of the displacement field. In the present work, we have adopted C 1 triangular Argyris elements [START_REF] Dominguez | A simple Matlab implementation of the Argyris element[END_REF] with 21 degrees of freedom per element for discretizing the displacements.

The 2D vector form associated with the components of ∇ε and of the hyperstress tensor S can be written as2 :

[∇ε] =

⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ∇ε 111 ∇ε 221 2∇ε 122 ∇ε 222 ∇ε 112 2∇ε 121 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ∂ 2 u 1 ∂x 2 1 ∂ 2 u 2 ∂x 1 ∂x 2 ∂ 2 u 1 ∂x 2 2 + ∂ 2 u 2 ∂x 1 ∂x 2 ∂ 2 u 2 ∂x 2 2 ∂ 2 u 1 ∂x 1 ∂x 2 ∂ 2 u 1 ∂x 1 ∂x 2 + ∂ 2 u 2 ∂x 2 1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , [S] = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ S 111 S 221 S 122 S 222 S 112 S 121 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ . (47)
In matrix form, the constitutive equations ( 21) and ( 22) are given by:

[σ] = C[ε] + M[∇ε], (48)
S = M T [ε] + G[∇ε]. (49)
with

C = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ C 1111 C 1122 C 1112 C 1122 C 2222 C 2212 C 1112 C 2212 C 1212 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , (50)
M = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ M 11111 M 11221 M 11122 M 11222 M 11112 M 11121 M 22111 M 22221 M 22122 M 22222 M 22112 M 22121 M 12111 M 12221 M 12122 M 12222 M 12112 M 12121 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , (51)
and

G = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ G 111111 G 111221 G 111122 G 111222 G 111112 G 111121 G 221221 G 221122 G 221222 G 221112 G 221121 G 122122 G 122222 G 122112 G 122121 G 222222 G 222112 G 222121 Sym. G 112112 G 112121 G 121121 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , (52)
where we have taken into account the index symmetries of these tensors as defined by the relation [START_REF] Dominguez | A simple Matlab implementation of the Argyris element[END_REF].

The vectors [ε] and [δε] are related to nodal displacements and test functions in the element e through:

[ε] = Bu e , [δε] = Bδu e , (53)
where

B = D (1) ; D (2) ; ... ; D (n) , with D (I) = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ∂N I (X) ∂X 1 0 0 ∂N I (X) ∂X 2 ∂N I (X) ∂X 2 ∂N I (x) ∂X 1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (54)
where n denotes the number of nodes of the element. Then the associated vectors ∇ε(u) and ∇ε(δu) are related to displacement and test functions through:

∇ε(u) = Bu e , ∇ε(δu) = Bδu e , (55)
where B = D(1) ; D(2) ; ... ; D(n) , ,with

D(I) = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ∂ 2 N I (x) ∂x 2 1 0 0 ∂ 2 N I (x) ∂x 1 ∂x 2 ∂ 2 N I (x) ∂x 2 2 ∂ 2 N I (x) ∂x 1 ∂x 2 0 ∂ 2 N I (x) ∂x 2 2 ∂ 2 N I (x) ∂x 1 ∂x 2 0 ∂ 2 N I (x) ∂x 1 ∂x 2 ∂ 2 N I (x) ∂x 2 1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ . ( 56
)
Introducing the above discretization in [START_REF] Monchiet | Strain-gradient homogenization: a bridge between asymptotic expansion and quadratic boundary condition methods[END_REF] we obtain

where N T is the matrix of shape functions such that u(x) = N(x)u e . Using the arbitrariness of the trial function δu e it yields the linear system:

Ku = F (58)
with

K = Ω B T CB + B T M B + BT M T B + BT G B dΩ, (59)
and

F = ∂Ω F N T F * dΓ. ( 60
)

Numerical evaluation of the effective tensors

The discrete form associated with [START_REF] Huet | Application of variational concepts to size effects in elastic heterogeneous bodies[END_REF] is given by:

[ε(x)] = A 0 (x) [ε] + A 1 (x) -A 0 x (x) ∇ε (61)
with

A 0 (x) = B(x)U, A 1 (x) = B(x)W, A 0 x (x) = B(x)W x (x), (62)
where B(x) is the matrix of shape functions derivatives in one element of the RVE and U, W and W x (x) are defined as:

U = u 1 ; u 2 ; u 3 , W = u 4 ; u 5 ; u 6 ; u 7 ; u 8 ; u 9 , W x (x) = xu 1 ; yu 1 ; xu 2 ; yu 2 ; xu 3 ; yu 3 . ( 63
)
The displacement fields u i are the vector columns containing the nodal displacement solution of the problems (11)-( 13) formulated with the BC contained in (1, 0, 0) (0, 0, 0, 0, 0, 0) u 2 (0, 1, 0) (0, 0, 0, 0, 0, 0) u3 (0, 0, 1/2) (0, 0, 0, 0, 0, 0) u4 (0, 0, 0) (1, 0, 0, 0, 0, 0) u 5 (0, 0, 0) (0, 1, 0, 0, 0, 0) u 6 (0, 0, 0) (0, 0, 1, 0, 0, 0) u 7 (0, 0, 0) (0, 0, 0, 1, 0, 0) u 8 (0, 0, 0) (0, 0, 0, 0, 1, 0) u 9 (0, 0, 0) (0, 0, 0, 0, 0, 1) Table 2 Sets of elementary solution and corresponding activated strain and strain gradient components.

Introducing ( 61) in ( 35), ( 36) and (37), we obtain:

C = U T B T (x)C(x)B(x)U , (64)
M = U T B T (x)C(x)B(x) [W -W x ] , (65)
G = [W -W x ] T B T (x)C(x)B(x) [W -W x ] . ( 66
)
6 Numerical examples

In the following examples, the different RVEs presented in Fig. 2 will be investigated. In what follows, the different phases are supposed filled by an isotropic elastic material, and we denote by (E m , ν m ) and (E i , ν i ) the Young's moduli and the Poisson's ratios of the matrix and inclusions, respectively.

The RVE with the circular inclusion (see Fig. Unless otherwise specified, the length of the RVE is chosen as ℓ = 1 mm.

The proposed numerical homogenization scheme will now be validated in different situations. Firstly, its results will be compared to those obtained using asymptotic analysis. Secondly, the shape of the effective tensors will be analyzed for different anisotropic RVEs, and finally the scaling of the effective parameters will be studied. For all of these studies, the homogenization scheme ( 11)-( 13) will be applied considering the Effective Body Force Correction for which C 0 = C, and the RVE will consist in a cluster of 4 × 4 unit cells.

Comparison with asymptotic solution

We first validate the present model by comparing it with the solution obtained from asymptotic expansion analysis as obtained in [START_REF] Tran | A micromechanics-based approach for the derivation of constitutive elastic coefficients of strain-gradient media[END_REF]. In the mentioned work, results obtained by asymptotic analysis have been computed on periodic microstructures using Fast Fourier Transform. For this comparison, we consider the RVE with circular inclusions (see Fig. 2 (a)) and a volume fraction varying as a function of the ratio between the radius of the inclusions R and the length of the unit cell ℓ. In this study, the material parameters are E m = 1, ν m = 0.45 for the matrix and E i = 10, ν i = 0.3 for the inclusion. It is worth noting that here the elastic parameters are different than in [START_REF] Tran | A micromechanics-based approach for the derivation of constitutive elastic coefficients of strain-gradient media[END_REF]. Then, the reference solution presented in Fig. 3 and referred to as Asymptotic analysis has been re-computed using the FFT procedure described in [START_REF] Tran | A micromechanics-based approach for the derivation of constitutive elastic coefficients of strain-gradient media[END_REF] but with the above elastic parameters for direct comparison. Results provided in Fig. 3 show a good agreement between both models. It is worth noting that results obtained from the asymptotic analysis correspond to periodic cells, and thus to converged values of the effective coefficients with respect to the number of unit cells.

Symmetry classes

In this study, we show that the matrix forms M and G computed for anisotropic RVEs with the present scheme are in agreement with the theoretically matrix shapes predicted in [START_REF] Auffray | A complete description of bidimensional anisotropic strain-gradient elasticity[END_REF]. We first consider the RVE of Fig. 2 [4] that M = 0 and

G = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ G 111111 G 111221 G 111122 0 0 0 G 221221 G 221122 0 0 0 G 122122 0 0 0 G 111111 G 111221 G 111122 Sym. G 221221 G 221122 G 122122 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (N). ( 67
)
The numerical values (not reported here) of components of M are found to be of the order of 10 -5 which are negligible as compared to the values found in the next cases (see ( 73), (C.2) and (C.5)). The numerical values of C and G are given by:

C = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣
2.221 0.804 0.000 0.804 2.221 0.000 0.000 0.000 0.601

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (MPa), (68)
and

G = 10 -2 × ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣
1.779 -0.627 0.163 0.000 0.000 0.000 -0.627 0.535 -0.023 0.000 0.000 0.000 0.163 -0.023 0.949 0.000 0.000 0.000 0.000 0.000 0.000 1.779 -0.627 0.163 0.000 0.000 0.000 -0.627 0.535 -0.023 0.000 0.000 0.000 0.163 -0.023 0.949

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (N).( 69
)
It can further be checked that both matrix are positive definite, hence ensuring the same property for the overall elastic energy. Comparing (69) and ( 67), it can be verified that the homogenized matrices comply with the expected symmetries. For another example, we consider the RVE of Fig. 2 (b) containing triangular inclusions. For illustration, the elementary local strain states are depicted in Fig. 4. For such RVE (Z π 2 -symmetry), M is given by [START_REF] Auffray | A complete description of bidimensional anisotropic strain-gradient elasticity[END_REF]:

M = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ M 11111 M 11221 M 11122 0 0 0 M 22111 M 22221 M 22122 0 0 0 0 0 0M 12222 M 12112 M 12121 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (N/mm), (70)
and G, which is orthotropic in this situation (D 2 -invariance) is given by We obtain, using the same elastic parameters as in the previous example:

G = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ G 111111 G 111221 G 111122 0 0 0 G 221221 G 221122 0 0 0 G 122122 0 0 0 G 222222 G 222112 G 222121 Sym. G 112112 G 112121 G 121121 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (N). ( 71
) (a) (b) (c) (d) (e) (f) (g) (h) (i)
C = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣
2.822 0.784 0.000 0.784 2.586 0.000 0.000 0.000 0.748

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (MPa), (72)
M = 10 -2 × ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣
-8.432 0.641 -1.119 0.000 0.002 0.000 -0.893 -0.220 -0.414 0.000 0.000 0.000 0.000 0.000 0.001 0.370 -1.988 0.217

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (N/mm), (73)
and

G = 10 -3 × ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣
6.124 -0.876 0.656 0.000 0.000 0.000 -0.876 0.672 -0.065 -0.001 0.000 0.000 0.656 -0.065 1.160 0.000 0.000 0.000 0.000 -0.001 0.000 2.883 -1.071 -0.160 0.000 0.000 0.000 -1.071 1.553 -0.083 0.000 0.000 0.000 -0.160 -0.083 1.501

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (N).( 74
)
Here again, by comparing ( 73) and (70), and ( 74) and ( 71), it can be concluded that M and G possess the expected invariance properties. Finally, the RVE of Figs. 2 (c) (anisotropic RVE) and 2 (d) (arrow-shaped) have been analysed and, for the sake of conciseness, their numerical values postponed to Appendix C.

Size effects

Asymptotic analysis predicts that the values of M and G vary according to ϵ and ϵ 2 [START_REF] Boutin | Microstructural effects in elastic composites[END_REF]. In this example, we show that this predicted size-effects are accurately captured by the proposed model. To that aim, we consider the RVE of Fig. 2 (b) with triangular inclusions and we use the same elastic parameters as in the previous subsection. Here again, the RVE is composed of 4 × 4 unit cells. The dimensions of the RVE are varied according to ϵ = ℓ/ℓ 0 , where ℓ 0 = 1 mm and ℓ is the length of the RVE. It appears clearly from figure 5 that M and G vary according to ϵ and ϵ 2 , respectively, as it should.

Comparison with other schemes

In this last subsection, the homogenization scheme involving Pointwise Body Force Correction is compared with schemes involving other kind of body forces:

Null body force (Standard QBC): C 0 (x) = O.
Pointwise body force correction:

C 0 (x) = C(x).
More specifically, different choices for C 0 in Eq. ( 9) are compared together with their respective influence on the obtained values of M and G.

For this study, a RVE made of several periodic N × N unit cells with circular inclusions is considered (see Fig. 2 (a)), the elastic properties of the matrix It can be observed from Figs. 6 (a-b) that using C 0 = 0 and C 0 = C(x) lead to divergent values of G 111111 with respect to N . These effects have also been observed in other studies, such as in [START_REF] Barboura | Establishment of strain gradient constitutive relations by using asymptotic analysis and the finite element method for complex periodic microstructures[END_REF] in the context of approaches combining asymptotic analysis and finite element calculations. In contrast, the choice C 0 = C leads to convergent values of the effective component.

Secondly, we compare the different choices of C 0 for a variation of phase contrast of properties between the inclusion and the matrix. We still consider the RVE containing circular inclusions with a volume fraction f = 0.3 . This time a 1 × 1 unit cell is considered in order to avoid the above-mentioned divergence problems with respect to N . Results are presented in Figs. 78. It can be noted from Fig. 7 (a) that the solution corresponding to C 0 = 0 does not lead to a vanishing the sixth-order tensor G when the material is homogeneous (E i /E m = ν i /ν m = 1). On the contrary, as illustrated by Fig. 6 (b), the choice of C 0 (x) = C(x) implies a null effective sixth-order tensor G for homogeneous medium but leads to divergent higher-order properties for hard inclusions. Results corresponding to the case in which C 0 = C (Fig. 6 (c)) lead to zero effective sixth-order tensor G in the same situation. However, it can be observed that for soft inclusions, the effective higher-order properties diverge.

It is important to mention that this solution corresponds to the asymptotic solution obtained by Tran et al. [START_REF] Tran | A micromechanics-based approach for the derivation of constitutive elastic coefficients of strain-gradient media[END_REF]. However, the case of soft inclusions was so far not discussed and this issue not reported. The conclusions of this numerical analysis are as follows:

(1) Using the present corrected scheme with body forces involving C 0 (x) = C leads to convergent values of the higher-order properties with respect to the size of the RVE (number of unit cells) while other solutions (no body forces or C 0 (x) = C(x)) lead to properties which diverge with respect to the size of the RVE. (2) Using quadratic boundary conditions without body forces leads to nonzero higher-order effective properties for homogeneous media and to an effective strain-gradient medium even if the material is made of one phase with local Cauchy behavior which is contradictory; (3) The solutions with C 0 (x) = C and C 0 (x) = C(x) both satisfy vanishing higher-order effective properties for homogeneous RVE but lead to divergent properties for soft and hard inclusions, respectively. Regarding the present scheme (C 0 = C), this constitutes an issue shared with asymptotic expansion solutions [START_REF] Tran | A micromechanics-based approach for the derivation of constitutive elastic coefficients of strain-gradient media[END_REF].

Structure calculations

In this last subsection, the approach is applied to structure calculations and we study how the strain-gradient homogenized model is able to capture the higher-order effects induced by the coarseness of the internal architecture.

The result of such a study allows to decide whether it is necessary, or not, to include higher-order descriptors in the continuous description of a heterogeneous continuum. It is expected from asymptotic analysis that strain-gradient effects should have a greater importance as soon as non centro-symmetric microstructure are concerned.

x y 0 L H Fig. 10. Heterogeneous structure: geometry and boundary conditions.

To investigate this point, the structure described on Fig. 10 is considered.

To obtain a strong coupling effect, the unit cell can be either the complete anisotropic one (Fig. 2 (c)) or the arrow-shaped one (Fig. 2 (d)). The length of the structure is L = 5ℓ 0 and H = 2ℓ 0 . We first consider the situation in which the structure is made out of the full anisotropic unit cell. On the one hand, we perform a complete structure calculation in which all phases are meshed, and on the other hand, a computation on the same structure filled out with the strain gradient effective material as determined by the current approach. The complete structure calculation is made using a standard FEM formulation, while the effective computation involves the specific numerical scheme described in section 5.2.

Firstly, we consider the anisotropic RVE of Fig. It can be observed that there is no significant difference in the deformation mode between the strain-gradient model and the Cauchy medium. In this case the effect of the classical anisotropy elasticity is dominant as compared to strain-gradient effects4 . Secondly, we consider that the structure is generated from the unit cell depicted in Fig. 2 Even if the effects are extremely small, it can be observed in Figs. 12 (a) and (b) the bending of the structure under the action of the strain in the ydirection due to the strain-gradient effects. It can be concluded from Fig. 12 (b) that the strain-gradient model is able to reproduce the global bending of the structure associated with strain in the y-axis, while the Cauchy model is not (see Fig. 12 (c)).

Conclusion

In this work, we have proposed a computational homogenization framework to model anisotropic strain-gradient effects in elastic structures. The contributions of this work are as follows: (a) we have introduced corrections in the homogenization scheme to ensure that the strain-gradient effects vanish when the material is homogeneous with a local Cauchy behavior, which is usually not the case in the approaches reported in the literature and to avoid divergence of higher-order properties with respect to the size of the RVE, which is found when using classical quadratic boundary conditions; (b) we have validated the present scheme with respect to asymptotic expansion solutions and available results on higher-order tensors symmetries for different classes of anisotropy; (c) we have provided a C 1 Finite Element Framework to solve the strain-gradient macroscopic model and have compared some direct numerical simulations with the homogenized model to show its validity. Finally, we have pointed out that the present framework is valid for stiff inclusions but that for soft inclusions the effective higher-order properties diverge with respect to the contrast of properties. This issue remains open and shall be the topic of future extensions.

Concerning the importance of strain gradient effects in the overall description, and hence the necessity to incorporate them, or not, in an effective model, the study reveals that (a) in statics , (b) for the microstructures considered in the present paper, the correction remains small. The correction due to the introduction of the sixth-order tensor does not seem to be of much influence, and strain-gradient effects become perceptible only for non centro-symmetric structure for which the fifth-order elasticity tensor plays a role. Nevertheless as shown by the different computations its presence seem to be necessary but not sufficient and a specific micro-structure such as the arrow-shaped one is necessary to maximize its contribution. Then, definition of discriminant benchmark tests for numerical validation of the effective sixth-order tensor remains to be addressed. To enhance gradient effects in statics its seems necessary to (a) decrease the influence of the standard elasticity and to tend toward a degenerated first order elasticity such as in the pantographic structure [START_REF] Seppecher | Linear elastic trusses leading to continua with exotic mechanical interactions[END_REF][START_REF] Dell'isola | Pantographic metamaterials: an example of mathematically driven design and of its technological challenges[END_REF][START_REF] Abdoul-Anziz | Strain gradient and generalized continua obtained by homogenizing frame lattices[END_REF], (b) increase the contrast between phases. However, it is important to note that theses comments are restricted to elastostatic. Indeed, as evidenced for instance in [START_REF] Rosi | Anisotropic and dispersive wave propagation within strain-gradient framework[END_REF][START_REF] Rosi | Continuum modelling of frequency dependent acoustic beam focusing and steering in hexagonal lattices[END_REF], when it comes to dynamics and wave propagation strain gradient contributions which have vanishingly low effects in statics can be of great importance.

A Appendix: Relationship between second gradient of displacements and strain-gradient tensors

The relationship between G and ∇ε is provided as follows:

G ijk = ∂ 2 u i ∂x j ∂x k = 1 2
∂ 2 u i ∂x j ∂x k + ∂ 2 u i ∂x j ∂x k + ∂ 2 u j ∂x i ∂x k - ∂ 2 u j ∂x i ∂x k + ∂ 2 u k ∂x i ∂x j - ∂ 2 u k ∂x i ∂x j = 1 2 ∂ 2 u i ∂x j ∂x k + ∂ 2 u j ∂x i ∂x k + ∂ 2 u i ∂x k ∂x j + ∂ 2 u k ∂x i ∂x j - ∂ 2 u j ∂x k ∂x i - ∂ 2 u k ∂x j ∂x i = 1 2 ∂ 2 u i ∂x j ∂x k + ∂ 2 u j ∂x i ∂x k + 1 2 ∂ 2 u i ∂x k ∂x j + ∂ 2 u k ∂x i ∂x j - 1 2
∂ 2 u j ∂x k ∂x i - ∂ 2 u k ∂x j ∂x i
= ∇ε ijk + ∇ε ikj -∇ε jki .

(A.1)

B Appendix: Quadratic boundary conditions

We can show that the displacement field compatible with a linear strain field in the form

ε ij (x) = ∇ε ijk x k (B.1)
is given by:

u i = 1 2 G ijk x j x k , (B.2)
as shown below. Starting from

ε ij = 1 2 ∂u i ∂x j + ∂u j ∂x i (B.3)
and using (A.1), we have:

∂u i ∂x j = 1 2
G ipq (δ pj x q + x p δ qj ) = 1 2 G ijq x q + G ipj x p (B.4)

∂u j ∂x i = 1 2
G jpq (δ pi x q + xpδ qi ) = 1 2 G jiq x q + G jpi x p . (B.5)

Note that from (3) G ijp ̸ = G jip but G ijp = G ipj . Then: and

ε ij = 1 2 G ijp x p + G
G = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣
12.353 -0.076 0.086 0.000 0.000 0.000 -0.076 0.002 0.000 0.000 0.000 0.000 0.086 0.001 0.008 0.000 0.000 0.000 0.000 0.000 0.000 0.119 -0.006 -0.003 0.001 0.000 0.000 -0.006 0.078 -0.041 0.000 0.000 0.000 -0.003 -0.041 0.055

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (N). (C.6)

  GL(d): the group of all linear invertible transformations of R d , i.e. F ∈ GL(d) iff det(F) ̸ = 0; • O(d): the orthogonal group, which is the group of all isometries of R d i.e. Q ∈ O(d) iff Q ∈ GL(d) and Q -1 = Q T , where the superscript T denotes the transposition; • SO(d): the special orthogonal group, i.e. the subgroup of O(d) consisting of transformations satisfying det(Q) = 1.

Fig. 1 .

 1 Fig. 1. (a) heterogeneous structure; (b) equivalent strain-gradient homogeneous structure; (c) RVE.

Ω(

  δu e ) T B T CBu e dΩ + Ω (δu e ) T B T M Bu e dΩ + Ω (δu e ) T BT M T Bu e dΩ + Ω (δu e ) T BT G Bu e dΩ = ∂Ω F(δu e ) T N T F * dΓ,

Fig. 2 .

 2 Fig. 2. (a) RVE with circular inclusions; (b) RVE with triangular inclusions; (c) Asymetric RVE; (d) arrow-shaped RVE.the RVE of Fig.2 (d), which is also Z π 2 , is expected to produce large values for the coupling tensor M. The arrow-shape is such that under uniaxial strain along the x-direction the two upper bars are stretched and create a vertical displacement of the vertical bar to the bottom leading to strain-gradient in the y-direction. The geometric description of the different RVEs is provided below.

Fig. 3 .

 3 Fig.3. Comparison between the present model and effective moduli obtained from asymptotic analysis and FFT method[START_REF] Tran | A micromechanics-based approach for the derivation of constitutive elastic coefficients of strain-gradient media[END_REF].

Fig. 4 .

 4 Fig. 4. ε 11 (x) strain field in deformed (×0.2) configurations for RVE with triangular inclusions: (a) [ε] = [1; 0; 0], ∇ε = 0; (b) [ε] = [0; 1; 0], ∇ε = 0, (c) [ε] = [0; 0; 1/2], ∇ε = 0, (d) [ε] = 0, ∇ε = [1; 0; 0; 0; 0; 0]; (e) [ε] = 0, ∇ε = [0; 1; 0; 0; 0; 0]; (f) [ε] = 0, ∇ε = [0; 0; 1; 0; 0; 0], (g) [ε] = 0, ∇ε = [0; 0; 0; 1; 0; 0], (h) [ε] = 0, ∇ε = [0; 0; 0; 0; 1; 0], (i) [ε] = 0, ∇ε = [0; 0; 0; 0; 0; 1] .

Fig. 5 .

 5 Fig. 5. Size-dependent effective strain-gradient properties for the RVE with triangular inclusions.are E m = 1 Mpa and ν = 0.3. In a first case, we chose the properties of the inclusions as E i = 10 3 MPa, ν i = 0.3. the retained volume fraction is f = 0.3.

Fig. 6 .

 6 Fig. 6. RVE with stiff circular inclusions (E i /E m = 10 3 ), convergence of the values of G 111111 for: (a): C 0 = C(x); (b) C 0 = 0; (c) C 0 = C.

Fig. 7 .Fig. 8 .

 78 Fig. 7. Values of G 111111 versus contrast of phase properties E i /E m , C 0 = 0.

Fig. 9 .

 9 Fig. 9. Values of G 111111 versus contrast of phase properties E i /E m , C 0 = C.

  2 (c). The material properties of the different phases are provided in Appendix C together with the tensors describing the effective properties ((C.1)-(C.4)).Concerning the boundary conditions, on the left boundary (x = 0) displacement are blocked, i.e. u x = u y = 0, while on the right boundary (x = L) the longitudinal displacement is prescribed u x = 2.5 × 10 -2 mm (u y is left free).The results of the different simulations are depicted on Figs. 11. The subfigures (a)-(b) represent the deformed configuration (exaggerated 50 times) of the structure for the reference and the homogeneous model. It can be observed that the deformed configuration is qualitatively well reproduced by the homogeneous model. However, in Fig.11 (c) we plot the solution of the homogenized model when we remove the strain-gradient effects (effective Cauchy medium).

Fig. 11 .

 11 Fig. 11. (a) Reference solution; (b) present homogeneous strain-gradient model; (c) Cauchy medium (ε 22 strain field); displacements exaggerated 50 times.

  (d). The material properties of the different phases are provided in Appendix C. Here, concerning the boundary conditions, on the left boundary (x = 0) displacement is blocked i.e. u x = u y = 0, , while on the right boundary (x = L) the longitudinal displacement is prescribed u x = 2.5 × 10 -2 mm and the transverse one blocked, u y = 0.Comparisons between the reference model (direct numerical simulation), the present homogenized model and the effective Cauchy medium are reported on Fig.12.

Fig. 12 .

 12 Fig. 12. (a) Reference solution; (b) present homogeneous strain-gradient model; (c) Cauchy medium (ε 22 strain field); displacements exaggerated 50 times in the x-direction and 100 times in the y-direction to better appreciate the bending of the structure.

Fig. 2 ( 3 )

 23 c), the elastic properties of the matrix are E m = 1 MPa and ν m = 0.3, the elastic properties of the inclusions are E i = 10 2 MPa and ν i = 0.3. For such RVE, there is no symmetry and the matrices M and G are expected to be fully populated with independent constants. For this case, the values are provided below.For the RVE of Fig.2(d) (arrow-shape), the elastic properties of the matrix are E m = 1 MPa and ν m = 0.4, the elastic properties of the inclusions are E i = 10 3 MPa and ν i = 0.3. The values are provided below.

Table 2 :

 2 Name (ε 11 , ε 22 , ε 12 ) (∇ε 111 , ∇ε 221 ∇ε 122 , ∇ε 222 , ∇ε 112 , ∇ε 121 )

	u 1

  jpi x p . (B.6) Using (A.1) and ∇ε ijp = ∇ε jipε ij = 1 2 ∇ε ijp + ∇ε ipj -∇ε jpi + ∇ε jip + ∇ε jpi -∇ε ipj x p = ∇ε ijp x p .(B.7) ∇ε pjk (δ ji x k + x j δ ki ) = 1 2 ∇ε pik x k + ∇ε pji x j (B.10) ∇ε ipk x k + ∇ε ijp x j + ∇ε pik x k + ∇ε pji x j ∇ε ijp x j + ∇ε jpi x j ̸ = ∇ε ipk x k . (B.13) C Appendix: numerical values of specific RVEs The numerical values for the RVE of Figs. 2 (c) and 2 (d) are reported. For each case, a 4 × 4 unit cells RVE has been used. For the anisotropic RVE of

	∂u p ∂x i	=	1 2
	and
	ε ip =	1 2	∂u i ∂x p	+	∂u p ∂x i	(B.11)
	=	1 4						(B.12)
	=	1 2		∇ε ipk x k +	1 4
	On the contrary, the choice
	u i =		1 2	∇ε ijk x j x k	(B.8)
	does not lead to the strain field (B.1), as shown in the following:
	∂u i ∂x p	=	1 2	∇ε ijk (δ jp x k + x j δ kp ) =	1 2	∇ε ipk x k + ∇ε ijp x j ,	(B.9)

It should be observed, that for the body force free scheme, which corresponds to C 0 (x) = 0 in[START_REF] Barthelat | Architectured materials in engineering and biology: fabrication, structure, mechanics and performance[END_REF], the expression of A 1 (x) is not correct.

The notations retained here are different from the ones used in the references[START_REF] Auffray | Derivation of anisotropic matrix for bidimensional strain-gradient elasticity behavior[END_REF][START_REF] Auffray | A complete description of bidimensional anisotropic strain-gradient elasticity[END_REF], hence some cautions are needed when using the results coming from theses papers.

The overall symmetry is the intersection of the symmetry group of the inclusion and the one of the matrix, since the first one is D 3 -invariant while the other is D

invariant, only the mirror line of symmetry lies in their intersection, hence resulting

Such a conclusion is restricted to a static analysis. For elastodynamics this vanishingly low correction can produce wave propagation that strongly differs from the one described using Cauchy equivalent continuum[START_REF] Rosi | Anisotropic and dispersive wave propagation within strain-gradient framework[END_REF].