Julien Yvonnet
email: julien.yvonnet@univ-eiffel.fr

Qi-Chang He

Pengfei Li

A data-driven harmonic approach to constructing anisotropic damage models with a minimum number of internal variables

Keywords: Anisotropic damage, data-driven, phase field, craks, image-based models, induced anisotropy, harmonic damage models

A data-driven approach is proposed to construct anisotropic damage models with a minimal number of internal variables from numerical simulations on Representative Volume Elements (RVEs) of quasi-brittle materials. The approach resorts in particular to a harmonic analysis of damage. The orientation distribution functions of two elastic moduli are numerically determined while accounting for the effects of the nucleation and propagation of microcracks by the phase-field method. Given these two functions, the effective elastic tensor of a material without or with microcracks is uniquely determined. The expansions into two Fourier series of the relative variations of these two functions related to an undamaged reference state and to a damage state make appear damage internal variables naturally. The number and natures of these variables can be optimized by truncating the Fourier series according to the degree of approximation desired. Thus, 2D and 3D anisotropic damage models can be constructed without resorting to usual assumptions made in damage mechanics. This construction holds for complex microstructures including image-based ones and for arbitrary loading history. Two-and three-dimensional applications are provided to evaluate the accuracy of the damage models constructed and to show the potential of the approach proposed.

Introduction

Construction of anisotropic damage models is an important topic in engineering, due to the fact that the distribution of cracks in many materials, such as composites, bone tissues, geomaterials or concrete, often exhibits privileged directions. Kachanov [START_REF] Kachanov | Time of the rupture process under creep conditions[END_REF] and Rabotnov [START_REF] Rabotnov | Creep rupture[END_REF] formulated early theories of isotropic damage, which were then extensively developed by Lemaitre and coworkers [START_REF] Lemaitre | Aspect phénoménologique de la rupture par endommagement[END_REF][START_REF] Lemaitre | Mechanics of solid materials[END_REF]. Isotropic damage models are useful in some situations but incapable of correctly describing the damage of materials whose microstructure is initially strongly anisotropic or becomes strongly anisotropic due to loading [START_REF] Chow | An anisotropic theory of elasticity for continuum damage mechanics[END_REF][START_REF] Hayhurst | Creep rupture under multi-axial states of stress[END_REF][START_REF] Lee | An anisotropic damage criterion for deformation instability and its application to forming limit analysis of metal plates[END_REF]. If the importance of establishing anisotropic damage models has been recognized for a long time [START_REF] Cordebois | Damage induced elastic anisotropy[END_REF][START_REF] Chaboche | Le concept de contrainte effective appliqué à l'élasticité et à la viscoplasticité en présence d'un endommagement anisotrope[END_REF][START_REF] Kattan | Decomposition of damage tensor in continuum damage mechanics[END_REF], the relevant results reported in the literature are still quite limited.

The construction of anisotropic damage models has up to now run into difficulty due to the lack of suitable tools for the treatment of microcracks in materials and particularly in the ones with complex microstructures. It is now established that the phase field method to fracture [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF][START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF][START_REF] Amor | Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments[END_REF][START_REF] Miehe | A phase field model for rateindependent crack propagation: Robust algorithmic implementation based on operator splits[END_REF][START_REF] Ambati | A review on phase-field models of brittle fracture and a new fast hybrid formulation[END_REF] can overcome this difficulty. Indeed, this method is able to deal efficiently with the initiation, propagation and coalescence of multiple cracks in complex and realistic microstructures of heterogeneous materials such as those resulting directly from experimental images (see e.g. [START_REF] Nguyen | Large-scale simulations of quasi-brittle microcracking in realistic highly heterogeneous microstructures obtained from micro ct imaging[END_REF][START_REF] Liu | An efficient matrix-free preconditioned conjugate gradient based multigrid method for phase field modeling of fracture in heterogeneous materials from 3d images[END_REF]).

Another central and highly challenging problem is the choice of internal variables [START_REF] Cauvin | Damage mechanics: basic variables in continuum theories[END_REF][START_REF] Lemaitre | Anisotropic damage law of evolution[END_REF][START_REF] Voyiadjis | Some basic issues of isotropic and anisotropic continuum damage mechanics[END_REF]. The lack of uniformity and rigor in making this choice was pointed out, in particular, in a critical review of Rabier [START_REF] Rabier | Some remarks on damage theory[END_REF]. Different kinds of damage variables, ranging from scalar ones [START_REF] Lemaitre | A course on damage mechanics[END_REF][START_REF] Zhu | Micromechanical model for simulating the fracture process of rock[END_REF] to tensorial ones [START_REF] Chaboche | Development of continuum damage mechanics for elastic solids sustaining anisotropic and unilateral damage[END_REF][START_REF] Chen | On damage strain energy release rate Y[END_REF][START_REF] Murakami | Mechanical modeling of material damage[END_REF], have been proposed. Recently, Olsen-Kettle [START_REF] Olsen-Kettle | Bridging the macro to mesoscale: Evaluating the fourthorder anisotropic damage tensor parameters from ultrasonic measurements of an isotropic solid under triaxial stress loading[END_REF] has extended the work of Cauvin and Testa [START_REF] Cauvin | Damage mechanics: basic variables in continuum theories[END_REF] to the general case of anisotropic damage by using a fourth-order damage tensor. As argued in [START_REF] Kanatani | Distribution of directional data and fabric tensors[END_REF][START_REF] Rabier | Some remarks on damage theory[END_REF], only second or higher evenorder tensors can be adopted as variables to describe anisotropic damage. Indeed, the models [START_REF] Oller | Isotropic damage models and smeared crack analysis of concrete[END_REF][START_REF] Mazars | The unilateral behaviour of damaged concrete[END_REF][START_REF] Faria | A strain-based plastic viscous-damage model for massive concrete structures[END_REF][START_REF] Wu | An energy release rate-based plastic-damage model for concrete[END_REF] with one or more scalar damage variables are generally incapable of representing anisotropic damage. The models employing vectorial damage variables [START_REF] Krajcinovic | The continuous damage theory of brittle materials, part 1: general theory[END_REF][START_REF] Talreja | A continuum mechanics characterization of damage in composite materials[END_REF] are generally incorrect. Applications relying on second-order tensors damage models were reviewed in [START_REF] Wu | On the mathematical and thermodynamical descriptions of strain equivalence based anisotropic damage model[END_REF]. Models with a fourth-order [START_REF] Chaboche | Le concept de contrainte effective appliqué à l'élasticité et à la viscoplasticité en présence d'un endommagement anisotrope[END_REF][START_REF] Ortiz | A constitutive theory for the inelastic behavior of concrete[END_REF][START_REF] Simo | Strain-and stress-based continuum damage models: I. formulation[END_REF][START_REF] Yazdani | Combined plasticity and damage mechanics model for plain concrete[END_REF][START_REF] Lemaitre | Mechanics of solid materials[END_REF][START_REF] Lubarda | Damage tensors and the crack density distribution[END_REF][START_REF] Govindjee | Anisotropic modelling and numerical simulation of brittle damage in concrete[END_REF][START_REF] Meschke | An anisotropic elastoplastic-damage model for plain concrete[END_REF] or an eighth-order damage tensor [START_REF] Chaboche | Continuum damage mechanics: PartI and II[END_REF][START_REF] Cauvin | Damage mechanics: basic variables in continuum theories[END_REF] can be found in the literature. Obviously, the use of an eighth-order tensor damage variable leads to a huge level of complexity both in modelling and in identification.

One appealing approach to constructing a general anisotropic damage model was initiated by Ladèveze [START_REF] Ladevèze | Sur une théorie de l'endommagement anisotrope[END_REF] and substantially extended by He and Curnier [START_REF] He | A more fundamental approach to damaged elastic stress-strain relations[END_REF][START_REF] He | Characterising a 2D elasticity tensor by two orientation distribution functions[END_REF]. By this approach, the anisotropic damage state of a material is completely characterized by two damage orientation distribution functions; the elasticity tensor of the damaged material for a given damage state is totally determined in terms of the elasticity tensor of the damaged material and of the two damage orientation distribution functions. This approach allows avoiding the empirical notions and assumptions made in damage mechanics, such as the notion of effective stress (or strain) and the hypothesis of strain (or stress) equivalence, and can describe any anisotropic damage in the particular two-dimensional (2D) or general three-dimensional (3D) situation. In addition, with the help of generalized spherical harmonics, the two damage orientation distribution functions can be developed into two Fourier series where the coefficients in form of even-order tensors act naturally as damage internal variables. Truncating these series according to a desired degree of accuracy offers a consistent and straightforward way of reducing the number of damage internal variables. In this manner, isotropic damage models can be recovered and anisotropic damage models of different degrees of complexity can be constructed. Note that a related approach using fabric tensors and involving the evolution laws for damage tensors was proposed in [START_REF] Voyiadjis | Damage mechanics with fabric tensors[END_REF] with application to microcrack distributions.

One main difficulty met in using the aforementioned Ladevèze-He-Curnier approach [START_REF] Ladevèze | Sur une théorie de l'endommagement anisotrope[END_REF][START_REF] He | A more fundamental approach to damaged elastic stress-strain relations[END_REF][START_REF] He | Characterising a 2D elasticity tensor by two orientation distribution functions[END_REF] resides in determining the two damage orientation distribution functions in the general case where the microstructure of the material is complex and the loadings applied to the material are not proportional. Indeed, up to now, in such situations it has experimentally or numerically been very tough to quantify the effects of microcracks on the macroscopic elastic properties of the material and to deal with the nucleation and propagation of microcracks in the material.

In the present work, the approach proposed by Ladevèze, He and Curnier [START_REF] Ladevèze | Sur une théorie de l'endommagement anisotrope[END_REF][START_REF] He | A more fundamental approach to damaged elastic stress-strain relations[END_REF][START_REF] He | Characterising a 2D elasticity tensor by two orientation distribution functions[END_REF], named harmonic analysis-based damage models in this paper, is combined with a data-driven method to construct anisotropic damage models. The two damage orientation distribution functions are determined from data obtained numerically by calculations on Representative Volume Elements (RVE) of microstructures with microcracks and by computational homogenization; a numerical procedure is suggested to extract the optimal set of damage internal variables while accounting for microstructure, induced anisotropy, nucleation and propagation of microcracks; the evolution of the effective (or homogenized) elastic stiffness tensor in terms of the nucleation and propagation of microcracks produced by loading is computed by using the phase-field method [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF][START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF][START_REF] Amor | Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments[END_REF][START_REF] Miehe | A phase field model for rateindependent crack propagation: Robust algorithmic implementation based on operator splits[END_REF][START_REF] Ambati | A review on phase-field models of brittle fracture and a new fast hybrid formulation[END_REF]. The method thus elaborated has the following salient features. First, it allows optimizing the choice of internal variables to characterize a given damage state. Second, the optimized choice of damage internal variables accounts not only for the initial microstructure of a material but also for the evolution of microstructure due to the microcrack nucleation and propagation related to the whole history of the loadings undergone. Finally, the method is valid for any 2D and 3D microstructures and for any load-ing history. The present work contribute both to the development of damage mechanics and to the one of data-driven computational mechanics which, popularized by Kirchdoerfer and Ortiz [START_REF] Kirchdoerfer | Data-driven computational mechanics[END_REF], employs directly experimental data to construct a constitutive law in a finite element simulation without resorting to any empirical model. Remark that the basic idea of data-driven computational mechanics has been adopted recently in other topics such as fracture mechanics [START_REF] Carrara | Data-driven fracture mechanics[END_REF] and history-dependent materials [START_REF] Ladevèze | Data-driven computation for history-dependent materials[END_REF]. The use of data-driven techniques in multiscale modelling and computational homogenization was initiated in [START_REF] Yvonnet | Numerically explicit potentials for the homogenization of nonlinear elastic heterogeneous materials[END_REF][START_REF] Le | Computational homogenization of nonlinear elastic materials using neural networks[END_REF] concerning nonlinear elasticity, and extended to phenomena ranging from elasto-plasticity to nonlinear conduction (see a recent survey of datadriven multiscale approaches in [START_REF] Lu | A stochastic FE2 data-driven method for nonlinear multiscale modeling[END_REF]). However, to the best of our knowledge, data-driven mechanics has not yet been applied and developed to construct anisotropic damage models.

The remainder of this paper is organized as follows. In section 2, a brief review on the computational homogenization of cracked RVEs, which is one of the central ingredients of the proposed approach, is provided. In section 3, the one-to-one correspondence between the effective elastic tensor and a pair of elastic modulus orientation distribution functions together with their Fourier series representations is recalled. In section 4, we describe the harmonic analysis-based damage models and a method to extract a minimal set of internal variables from RVE calculations. The full data-driven procedure is summarized in section 5. Finally, numerical examples are provided in section 6 to show the capabilities, convergence properties and accuracy of the method, with applications to RVE-based microstructures involving image-based models.

The main notations used in this paper are specified as follows. Vectors and second-order tensors are denoted by bold letters; fourth-order tensors are denoted by double case letters. Simple, double and twice double contractions of indices are carried out as follows:

(A • v) i = A ij v j , A : B = A ij B ij , (A : B) ij = A ijkl B kl , A :: B = A ijkl B ijkl ,
where v denotes a vector, A and B represent two second-order tensors, and A and B are two fourthorder tensors; the following tensor products are defined as:

(u ⊗ v) ij = u i v j , (A ⊗ B) ijklmp = A ijkl B mp , (A⊗B) ijkl = (A ik B jl + A il B jk ).
Finally, ∇(.) denotes the gradient operator and ∇ • (.) denotes the divergence operator.

In the present work, strains are assumed to be small and only quasi-brittle materials are considered.

Computational homogenization of cracked RVEs

The basics of linear computational homogenization on a Representative Volume Element (RVE) (see e.g. [START_REF] Michel | Effective properties of composite materials with periodic microstructure: a computational approach[END_REF][START_REF] Zohdi | An introduction to computational micromechanics[END_REF][START_REF] Yvonnet | Computational Homogenization of Heterogeneous Materials with Finite Elements[END_REF]) are first reviewed. These calculations will be used in the sequel to compute the effective elastic tensor C(t) at each load step (denoted here as time t) during the preliminary calculations and will serve as data to extract the internal variables. Let us consider an RVE occupying a domain Ω ⊂ R dim , with dim the space dimension (see Fig. 1). We first define the following micro-macro relationship between the micro strain tensor ε and the macro strain tensor ε ( In the following, the quantities denoted by (.) will be associated with macroscopic quantities):

ε = 1 |Ω| Ω εdΩ, (1)
which holds only if the RVE does not contain sharp cracks, imperfectly bounded interfaces or holes and where ε = ∇u + ∇ T u . In the case where holes or sharp cracks are present within the RVE, the previous definition must be replaced by

ε = 1 2|Ω| Ω (u ⊗ x + x ⊗ u) dΓ, (2)
where u denotes the displacements vector and x the position vector of a point within the RVE. The relation between the micro stress tensor σ and the macro stress tensor σ is given by

σ = 1 |Ω| Ω σdΩ (3)
in the absence of sharp cracks or cohesive interfaces, or

σ = 1 2|Ω| Ω [(σ • n) ⊗ x + x ⊗ (σ • n)]dΓ (4)
if the RVE does contain sharp cracks or cohesive interfaces. The effective properties can be computed by prescribing, e.g., the following boundary conditions on the surface ∂Ω of the RVE

u(x) = ε • x + ũ(x) on ∂Ω, (5)
where ũ represents a periodic fluctuation over Ω. Solving the linear problem

∇ • σ = 0 on ∂Ω (6)
with ∇ • (.) denotes the divergence operator and σ(x) = C(x, t) : ε(x), the localization tensor at time t can be written as

A ijkl (x, t) = ε (kl) ij (x, t) (7)
where ε

(kl) ij
represents the strain solution of the above problem when ε = 1 2 (e k ⊗ e l + e l ⊗ e k ).

We can then define at each time t the following linear relationship

ε(x, t) = A(x, t) : ε(t) (8)
and

σ(x, t) = C(x, t) : A(x, t) : ε. (9)
Then, averaging over Ω, we obtain the macro stress-strain relationship:

σ(t) = C(t) : ε(t) (10)
with

C(t) = 1 |Ω| Ω C(x, t) : A(x, t)dΩ (11)
in the absence of sharp cracks or cohesive interfaces, or

C ijkl (t) = 1 2|Ω| [C ikst (x)A stpq (x, t)n k x j + C jkst (x)A stpq (x, t)n k x i ] dΓ (12)
if the RVE contains sharp cracks or cohesive interfaces. A more detailed description of the procedure to compute the effective elastic tensor with finite elements can be found e.g. in [START_REF] Yvonnet | Computational Homogenization of Heterogeneous Materials with Finite Elements[END_REF]. In the present work, only continuous (or regularized) description of cracks is used in the phase field method [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF][START_REF] Miehe | A phase field model for rateindependent crack propagation: Robust algorithmic implementation based on operator splits[END_REF]. Then, only definition (11) will be employed.

Assuming the evolution of C(t) known through computations on the RVE, a procedure for describing the damage of C(t) with a minimal number of internal variables is introduced in the next section. To this end, the notion of scalarvalued orientation distribution functions is developed in the next section. 3 Determination of the fourth-order elasticity tensor from two scalar orientation distribution functions

Scalar-valued orientation distribution functions

Let f be a macroscopic scalar property of the material. As defined in [START_REF] He | A more fundamental approach to damaged elastic stress-strain relations[END_REF], at one point of the macro structure, f depends generally on the orientation, specified by the unit vector n (see Fig. 1). Alternatively, f can be defined as a function of an angle θ in 2D in a polar coordinate system or by two angles (θ, φ) in 3D in a spherical coordinate system (see Fig. 19). We have

f : S -→ R (13)
with S denotes the unit circle (2D) or unit sphere (3D). The function f must satisfy the radial symmetry condition

f (n) = f (-n) ∀n ∈ S, (14)
because f must not be affected by the choice of n or -n for defining an orientation.

Assuming f to be square-integrable on S, f can be classically expanded into Fourier series in the 2D case as

f (θ) = a 0 + k (a k cos(kθ) + b k sin(kθ)) (15)
with

a 0 = 1 2π π -π f (θ)dθ, (16)
a k = 1 π π -π f (θ)cos(kθ)dθ, b k = 1 π π -π f (θ)sin(kθ)dθ (17)
where the integrals can be evaluated numerically. The foregoing expansion can be written in the following equivalent form [START_REF] Vilenkin | Fonctions spéciales et théories de la représentation des groupes[END_REF][START_REF] Jones | Spherical harmonics and tensors for classical field theory[END_REF][START_REF] Onat | Representation of mechanical behavior in the presence of changing internal structure[END_REF]:

f (θ) = g + G : F(θ) + G :: F(θ) + ... . (18)
where the second-and fourth-order tensors F(θ) and F(θ) denote generalized spherical harmonics, the scalar, second-and fourth-order tensors g, G and G are coefficients. The explicit expressions of F, F, g, G and G are specified in section 3.2.1 for the 2D case and in section 3.2.2 for the 3D case.

In [START_REF] Cauvin | Damage mechanics: basic variables in continuum theories[END_REF], {1, F, F, ...} is a complete orthogonal basis for square-integrable functions defined on S. Note that the orthogonality of the basis functions {1, F, F, ...} implies that

S FdS = 0, S FdS = O, (19)
S F ⊗ FdS = S F ⊗ FdS = O 6 ( 20
)
where dS = dθ in 2D and dS = sin(θ)dθdϕ in 3D with θ and ϕ denoting the Euler angles in a spherical coordinate system (see Fig. 19). Above, O 6 represents the sixth-order zero tensor. Additionally, F is symmetric and traceless (ST), i.e.

F T = F, T r(F) = 0, ( 21
)
and F is completely symmetric and traceless (CST), i.e.

F ijkl = F jikl = F klij = F ikjl (22)
and

F ijkk = 0. ( 23
)
Note that from (18) G is also ST and G is also CST.

Construction of the elasticity tensor from two orientation distribution functions

In [START_REF] He | A more fundamental approach to damaged elastic stress-strain relations[END_REF], He and Curnier formulated an approach to describe damaged elasticity using orientation functions. The approach, initially proposed by Ladevèze [START_REF] Ladevèze | Sur une théorie de l'endommagement anisotrope[END_REF][START_REF] Ladevèze | On an anisotropie damage theory[END_REF], is based on a theorem on the Fourier series representation of any elasticity tensor. It was shown that an elasticity tensor can be completely described by only two scalar distribution functions depending on the spatial orientation. Precisely, given an elasticity tensor C, the two orientation distribution functions giving the elongation and bulk moduli are defined by

η(N) = N : C : N = n i n j n k n l C ijkl (24)
κ(N) = I : C : N = δ ij n k n l C ijkl (25)
with N = n ⊗ n and I denotes the second-order identity tensor. In turn, these two orientation distribution functions determine C completely and uniquely.

Expanding η and κ in the Fourier series as in [START_REF] Cauvin | Damage mechanics: basic variables in continuum theories[END_REF] and identifying the relevant coefficient terms, we can reconstruct the tensor C from η and κ as described as follows.

2D case

In 2D, we express η and κ as functions of the angle θ (see Fig. 1 (a)). Classically,

n 1 = cos(θ), n 2 = sin(θ), (26)
and

F(θ) = N - 1 2 I. (27)
The fourth-order tensor F is given in 2D as [START_REF] He | Characterising a 2D elasticity tensor by two orientation distribution functions[END_REF]

F(θ) = N ⊗ N - 1 6 (I ⊗ N + N ⊗ I + I⊗N + N⊗I) + 1 24 (I ⊗ I + I⊗I) . (28)
In addition, g, G and G are specified in 2D by [START_REF] He | Characterising a 2D elasticity tensor by two orientation distribution functions[END_REF] g

= 1 2π π -π f (θ)dθ, G = 2 π π -π f (θ)F(θ)dθ, (29)
G = 8 π π -π f (θ)F(θ)dθ. ( 30
)
Expanding ( 24) and ( 25) in 2D gives

η(θ) = C 1111 n 4 1 + C 2222 n 4 2 + (2C 1122 + 4C 1212 )n 2 1 n 2 2 +4C 1112 n 3 1 n 2 + 4C 2212 n 1 n 3 2 (31)
and

κ(θ) = C 1122 + C 1111 n 2 1 + C 2222 n 2 2 + 2(C 1112 + C 2212 )n 1 n 2 . ( 32
)
We introduce the 2D matrices F and Ĝ of the fourth-elasticity tensors F and G such that F ::

G = F : Ĝ F(θ) =        n 4 1 -n 2 1 + 1 8 n 2 1 n 2 2 -1 8 n 3 1 n 2 -1 2 n 1 n 2 n 4 2 -n 2 2 + 1 8 n 1 n 3 2 -1 2 n 1 n 2 Sym. n 2 1 n 2 2 -1 8        (33)
and Ĝ denotes a (3 × 3) matrix. The elasticity tensor C having the two minor symmetries and a major symmetry, i.e. C ijkl = C jikl = C klij , it admits the following irreducible form (see e.g. [START_REF] Backus | A geometrical picture of anisotropic elastic tensors[END_REF][START_REF] Spencer | A note on the decomposition of tensors into traceless symmetric tensors[END_REF]):

C = λI ⊗ I + µI⊗I + I ⊗ A + A ⊗ I + I⊗A + A⊗1 + Z. ( 34
)
Then expressing η and κ as Fourier series in the form [START_REF] Cauvin | Damage mechanics: basic variables in continuum theories[END_REF]:

η(θ) u + V : F(θ) + W :: F(θ), (35)
κ(θ) v + V : F(θ), (36)
introducing [START_REF] Talreja | A continuum mechanics characterization of damage in composite materials[END_REF] in ( 24) and ( 25) and identifying the different terms with those of ( 35) and [START_REF] Ortiz | A constitutive theory for the inelastic behavior of concrete[END_REF], we obtain:

λ = v -u, µ = u - v 2 , A = 1 6 V, Z = W, (37)
with

u = 1 2π π -π η(θ)dθ, v = 1 2π π -π κ(θ)dθ, (38)
V = V = 2 π π -π η(θ)F(θ)dθ, W = 8 π π -π η(θ)F(θ)dθ. ( 39
)
For convenience and Finite Element applications purpose, a matrix form of the above expression is provided as

C =        C 1111 C 1122 C 1112 C 2222 C 2212 Sym. C 1212        (40)
and

C =        (λ + 2µ) λ 0 (λ + 2µ) 0 Sym. µ        +        6A 11 (A 11 + A 22 ) 3A 12 6A 22 3A 12 Sym. (A 11 + A 22 )        +        Z 1111 Z 1122 Z 1112 Z 2222 Z 2212 Sym. Z 1212        . (41)
From ( 18) and ( 37), then A is ST and Z is CST. Then

A 22 = -A 11 , (42)
Z 1122 = -Z 1111 , (43)
Z 2222 = Z 1111 , (44)
Z 2212 = -Z 1112 , (45)
Z 1212 = -Z 1111 . (46)
The 6 independent constants are set as

c 1 = λ, c 2 = µ, c 3 = A 11 , c 4 = A 12 , c 5 = Z 1111 , c 6 = Z 1112 . ( 47
)
Then the elasticity tensor can be re-written in the matrix form as

C =        (c 1 + 2c 2 ) c 1 0 (c 1 + 2c 2 ) 0 Sym. c 2        +        6c 3 0 3c 4 -6c 3 3c 4 Sym. 0        +        c 5 -c 5 c 6 c 5 -c 6 Sym. -c 5        .(48)

3D case

In 3D, the components of the unit vector n are expressed as a function of the Euler angles θ and ϕ in a spherical coordinate system (see Fig. 19):

n 1 = sin(θ)cos(ϕ), n 2 = sin(θ)sin(ϕ), n 3 = cos(θ). (49)
Then,

F(θ, ϕ) = N - 1 3 I, (50)
and the fourth-order tensor F is given by (see e.g. [START_REF] He | A more fundamental approach to damaged elastic stress-strain relations[END_REF])

F(θ) = N ⊗ N - 1 7 (I ⊗ N + N ⊗ I + I⊗N + N⊗I) + 1 35 (I ⊗ I + I⊗I) . ( 51
)
The coefficients of the Fourier series can be calculated according to [START_REF] Kanatani | Distribution of directional data and fabric tensors[END_REF]

g = 1 4π π θ=0 2π ϕ=0 f (θ, ϕ)sin(θ)dθdϕ, (52)
G = 15 8π π θ=0 2π ϕ=0 f (θ, ϕ)F(θ, ϕ)sin(θ)dθdϕ, (53)
and

G = 315 32π π θ=0 2π ϕ=0 f (θ, ϕ)F(θ, ϕ)sin(θ)dθdϕ. (54)
Expanding [START_REF] Chaboche | Development of continuum damage mechanics for elastic solids sustaining anisotropic and unilateral damage[END_REF] in 3D gives

η(θ, ϕ) = C : M 1 (θ, ϕ) (55)
with

C =                        C 1111 C 1122 C 1133 C 1112 C 1113 C 1123 C 2222 C 2233 C 2212 C 2213 C 2223 C 3333 C 3312 C 3313 C 3323 C 1212 C 1213 C 1223 C 1313 C 1323 C 2323 Sym.                        (56)
and

M 1 (θ, ϕ) =                        n 4 1 n 2 1 n 2 2 n 2 1 n 2 3 2n 3 1 n 2 2n 3 1 n 3 2n 2 1 n 2 n n 4 2 n 2 2 n 2 3 2n 1 n 3 2 2n 1 n 2 2 n 3 2n 3 2 n 3 n 4 3 2n 1 n 2 n 2 3 2n 1 n 3 3 2n 2 n 3 3 4n 2 1 n 2 2 4n 2 1 n 2 n 3 4n 1 n 2 2 n 4n 2 1 n 2 3 4n 1 n 2 n Sym. 4n 2 2 n 2 3                        . ( 57
)
Expressing [START_REF] Chen | On damage strain energy release rate Y[END_REF] in 3D, we obtain

κ(θ, ϕ) = C : M 2 (θ, ϕ) (58)
with

M 2 (θ, ϕ) =                        n 2 1 1 2 (n 2 1 + n 2 2 ) 1 2 (n 2 1 + n 2 3 ) n 1 n 2 n 1 n 3 n 2 n 3 n 2 2 1 2 (n 2 2 + n 2 3 ) n 1 n 2 n 1 n 3 n 2 n 3 n 2 3 n 1 n 2 n 1 n 3 n 2 n 3 0 0 0 Sym.                        . ( 59
)
We introduce the matrix forms F and Ĝ of the fourth-elasticity tensor F and G in 3D, respectively, such that F :: G = F : Ĝ as

F(θ, φ) =                        F 1111 F 1122 F 1133 F 1112 F 1113 F 1123 F 2222 F 2233 F 2212 F 2213 F 2223 F 3333 F 3312 F 3313 F 3323 F 1212 F 1213 F 1223 F 1313 F 1323 Sym. F 2323                        ( 60
)
and Ĝ is a 6×6 matrix. The closed-form expression for F is given in Appendix 8.

In 3D, the irreducible form of the elasticity tensor is given by [START_REF] Backus | A geometrical picture of anisotropic elastic tensors[END_REF]:

C = αI ⊗ I + βI⊗I + I ⊗ A + A ⊗ I + I⊗B + B⊗I + Z. (61)
Expanding η and κ into Fourier series as in [START_REF] Cauvin | Damage mechanics: basic variables in continuum theories[END_REF], introducing ( 61) in ( 24) and ( 25) and identifying the different terms with those of the Fourier series of η and κ, the elasticity tensor can be re-constructed from η and κ as [START_REF] He | A more fundamental approach to damaged elastic stress-strain relations[END_REF]

α = 1 2 (b -a), β = 1 4 (3a -b), (62)
A = V -U, B = 1 4 (3U -2V) , (63)
a = 1 4π π θ=0 2π ϕ=0 η(θ, ϕ)sin(θ)dθdϕ, (64) b
= 1 4π π θ=0 2π ϕ=0 κ(θ, ϕ)sin(θ)dθdϕ, (65)
U = 15 8π π θ=0 2π ϕ=0 F(θ, ϕ)η(θ, ϕ)sin(θ)dθdϕ, (66)
V = 15 8π π θ=0 2π ϕ=0 F(θ, ϕ)κ(θ, ϕ)sin(θ)dθdϕ, (67)
Z = 315 32π π θ=0 2π ϕ=0 F(θ, ϕ)η(θ, ϕ)sin(θ)dθdϕ. (68)
For finite element purpose, a matrix form of the elasticity tensor is introduced as

C = αI 1 + βI 2 + A + B + Z (69)
with

I 1 =                        1 1 1 0 0 0 1 1 1 0 0 0 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0                        , I 2 =                       
2 0 0 0 0 0 0 2 0 0 0 0 0 0 2 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0

0 0 0 0 0 1                        , (70)
A =                        2A 11 (A 11 + A 22 ) (A 11 + A 33 ) A 12 A 13 A 23 2A 22 (A 22 + A 33 ) A 12 A 13 A 23 2A 33 A 12 A 13 A 23 0 0 0 0 0 Sym. 0                        , (71)
and

B =                        4B 11 0 0 2B 12 2B 13 0 4B 22 0 2B 12 0 2B 23
4B 33 0 2B 13 2B 23
(B 11 + B 22 ) B 23 B 13 (B 11 + B 22 ) B 12
Sym.

(B 22 + B 33 )                        , ( 72
)
and Z is the matrix form associated with Z and is computed by

Z = 315 32π π θ=0 2π ϕ=0 F(θ, ϕ)η(θ, ϕ)sin(θ)dθdϕ. ( 73
)
4 Harmonic analysis-based damage models and extraction of internal variables

2D formulation

We denote by η(θ, t) and κ(θ, t) the elongation and bulk modulus orientation distribution functions at time t, and by η 0 (θ) ≡ η(θ, t = 0) and κ 0 (θ) ≡ κ(θ, t = 0) the corresponding one of the undamaged RVE. Following [START_REF] He | A more fundamental approach to damaged elastic stress-strain relations[END_REF], the two following damage orientation distribution functions d(θ, t) and h(θ, t) are introduced:

η(θ, t) = [1 -d(θ, t)] η 0 (θ) (74) κ(θ, t) = [1 -h(θ, t)] κ 0 (θ). (75)
Note that alternative forms may be adopted, like e.g. η(θ, t) = [1 -d(θ, t)] 2 η 0 (θ) and κ(θ, t) = [1 -h(θ, t)] 2 κ 0 (θ) to match popular phase field formulations [START_REF] Miehe | A phase field model for rateindependent crack propagation: Robust algorithmic implementation based on operator splits[END_REF] of gradient damage model. However, for the sake of simplicity, use is made only of ( 74

h(θ, t) = h 0 (t) + H(t) : F(θ) + ... (76)
with

d 0 (t) = 1 2π π -π d(θ, t)dθ, h 0 (t) = 1 2π π -π h(θ, t)dθ, (78)
and

D = 2 π π -π d(θ, t)F(θ)dθ, H = 2 π π -π h(θ, t)F(θ)dθ (79)
D = 8 π π -π d(θ, t)F(θ)dθ. (80)
These coefficients are computed numerically and stored as data for each time step. Then, for a given time t, we can reconstruct η(θ, t) and κ(θ, t) according to (74)-( 75) and the elasticity tensor by ( 37)- [START_REF] Carrara | Data-driven fracture mechanics[END_REF].

At this stage, the coefficients d 0 (t), h 0 (t), D(t) and D(t) constitute the set of internal variables. However, in practice, the expansions (81) and (82) can be truncated and reduced thanks to the symmetry properties of the tensors. Given the structure ( 35)-( 36), we propose the following approximation for (81)-(82):

d(θ, t) = d 0 (t) + D(t) : F(θ) + D(t) :: F(θ) + ... (81) h(θ, t) = h 0 (t) + D(t) : F(θ) + ... ( 82
)
As D is traceless (see ( 22)-( 23)), then,

D 22 = -D 11 . ( 83
)
As D is CST, then,

D 1122 = -D 1111 , (84)
D 2222 = D 1111 , (85)
D 2212 = -D 1112 , (86)
D 1212 = -D 1111 . (87)
Then, a maximum number of 6 independent internal variables is sufficient to describe a general anisotropic 2D damage state. For later use, let us set

α 1 = d 0 , α 2 =

3D formulation

Similarly, at each time step t, d(θ, ϕ, t) and h(θ, ϕ, t) are computed according to (74)-( 75) and are then expressed as Fourier series in the form

d(θ, ϕ) d 0 + D : F(θ, ϕ) + D :: F(θ, ϕ), (89)
h(θ, ϕ) h 0 + H : F(θ, ϕ), (90)
with

d 0 = 1 4π π θ=0 2π ϕ=0 d(θ, ϕ)sin(θ)dθdϕ, (91)
h 0 = 1 4π π θ=0 2π ϕ=0 h(θ, ϕ)sin(θ)dθdϕ, (93)
H = 15 8π π θ=0 2π ϕ=0 h(θ, ϕ)F(θ, ϕ)sin(θ)dθdϕ. (94)
These coefficients are computed numerically and stored as data for each time step. Note that these convolutions can be evaluated efficiently by pre-computing the associated linear operators (see [START_REF] Yvonnet | Harmonic analysis -based damage models: an efficient numerical framework[END_REF]). Then, for a given time t, η(θ, ϕ, t) and κ(θ, ϕ, t) can be reconstructed according to (89)-( 90) and ( 74)-( 75) and the elasticity tensor can be reconstructed from η and κ by ( 62)-( 69).

Here again, the whole set of internal variables associated with the coefficients d 0 (t), D(t), H(t), h 0 (t) and D(t) can be reduced by truncating the Fourier series (89)-(89). Additionally, because of the symmetry properties of the different tensors, a maximum of 21 internal variables can be defined as follows.

The tensors D and H are traceless. Then

D ii = 0, H ii = 0. (96)
The tensors D and H then contain 10 independent constants. As D is completely symmetric, then

D ijkl = D jikl = D ikjl = D ijlk . ( 97
)
This allows establishing the relationships

D 1212 = D 1122 , D 1213 = D 1123 , (98)
D 1223 = D 2213 , D 1313 = D 1133 , (99)
D 1323 = D 3312 , D 2323 = D 2233 . (100)
We can then eliminate 6 parameters among the 21 ones contained in D. In addition, D is traceless

D iikl = 0, (101)
which allows eliminating 6 more ones according to

D 3312 = -D 1112 -D 2212 , D 3313 = -D 1113 -D 2213 , (102)
D 3323 = -D 1123 -D 2223 , D 1133 = -D 1111 -D 1122 , (103)
D 2233 = -D 1122 -D 2222 , D 3333 = D 1111 + D 2222 + 2D 1122 . (104)
Finally, D can be expressed using only 9 independent parameters which are added to the 12 other internal variables. Finally, a maximum of 21 internal variables fully defines the damage state of an elastic anisotropic medium with microcracks. Note that in practice, this number can be much reduced, because some internal variables may play a negligible role, as shown in some following numerical examples.

For the sake of clarity, we re-name the 21 internal variables as

α 1 = d 0 , α 2 = h 0 , (105)
α 3 = D 11 , α 4 = D 12 , α 5 = D 13 , α 6 = D 22 , α 7 = D 23 , (106)
α 8 = H 11 , α 9 = H 12 , α 10 = H 13 , α 11 = H 22 , α 12 = H 23 , (107)
α 13 = D 1111 , α 14 = D 1122 , α 15 = D 2222 , (108)
α 16 = D 1112 , α 17 = D 1113 , α 18 = D 1123 , (109)
α 19 = D 2212 , α 20 = D 2213 , α 21 = D 2223 . (110)
The tensor D can then be recast in the matrix form such as D :: F = D : F:

D =                        (α 13 ) (α 14 ) (-α 13 -α 14 ) (α 16 ) (α 17 ) (α 18 ) (α 15 ) (-α 14 -α 15 ) (α 19 ) (α 20 ) (α 21 ) (α 13 + α 15 + 2α 14 ) (-α 16 -α 19 ) (-α 17 -α 20 ) (-α 18 -α 21 ) (α 14 ) (α 18 ) (-α 16 -α 19 ) (-α 13 -α 14 ) (-α 16 -α 19 ) (-α 14 -α 15 )                        .(111)
Finally, a simple criterion is provided to select the internal variables according to their cumulated relative values during the whole history of P sampling paths. Defining a pseudo time (load increment) t ∈ [0; T ] with T the final simulation time as, a selection criterion is adopted according to

α i max k (α k ) ≤ T ol, (112)
with

α i = P p=1 T 0 |α p i (t)|dt (113)
where α p i (t) denotes the i -th internal variable evolution sampled during the p -th sampling path and where T ol is a numerical prescribed tolerance.

Data-driven procedure

In the following, the different steps of the proposed data-driven procedure to construct the damage model from RVE calculations, are summarized. At each time step during the preliminary calculations, a numerical method must be used to simulate the crack propagation process within the RVE. In this work, the phase field method [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF][START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF][START_REF] Amor | Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments[END_REF][START_REF] Miehe | A phase field model for rateindependent crack propagation: Robust algorithmic implementation based on operator splits[END_REF][START_REF] Ambati | A review on phase-field models of brittle fracture and a new fast hybrid formulation[END_REF] is used, owing to its capability to model the initiation, propagation and interaction of cracks in complex microstructures [START_REF] Nguyen | A phase field method to simulate crack nucleation and propagation in strongly heterogeneous materials from direct imaging of their microstructure[END_REF][START_REF] Nguyen | Direct comparisons of 3D crack networks propagation in cementitious materials between phase field numerical modeling and in-situ microtomography experimental images[END_REF][START_REF] Nguyen | Large-scale simulations of quasi-brittle microcracking in realistic highly heterogeneous microstructures obtained from micro ct imaging[END_REF]. However, this choice is not restrictive and other methods can be used alternatively. For the sake of conciseness, we do not present in this paper a description about the phase field method and its implementation, and refer to the above references and the vast literature on this method. In what follows, the procedure is detailed only in the 2D case but can easily be extended to 3D by using the results of section 4.2.

Off-line calculations

The off-line calculations refer to preliminary calculations on the RVE to compute the effective tensor evolution C(t) and the evolution of the damage variables α(t) evolution. The steps of the procedure are given as follows.

(1) Define the RVE geometry and construct the RVE mesh.

(2) Define p = 1, 2, ..., P loading paths and associated macroscopic strains ε p (t) over a discrete set of loading steps t 1 , t 2 , ..., t N load .

• LOOP over all loading paths p (3) Compute the evolution of C(t)

LOOP over all loading steps t k Apply the boundary conditions [START_REF] Lee | An anisotropic damage criterion for deformation instability and its application to forming limit analysis of metal plates[END_REF] with ε p (t k ) and solve the fracture problem to obtain the local elasticity field C(x, t k ) within the RVE.

Compute the effective elasticity tensor C(t k ) using [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF] and store.

END LOOP (4) Compute internal variables evolution

Compute η 0 (θ) and κ 0 (θ) using ( 31)- [START_REF] Wu | An energy release rate-based plastic-damage model for concrete[END_REF].

LOOP over all loading steps k; given C(t k ): Compute d 0 (t k ), h 0 (t k ), D(t k ) and D(t k ) using ( 81)-( 82). Store the internal variables α

1 (t k ) = d 0 (t k ), α 2 (t k ) = h 0 (t k ), α 3 (t k ) = D 11 (t k ), α 4 (t k ) = D 12 (t k ), α 5 (t k ) = D 1111 (t k ) and α 6 (t k ) = D 1112 (t k ). END LOOP • END LOOP (5) Select relevant internal variables
• Compute α i according to (113).

• Select the relevant internal variables using the criterion (112) and the associated internal variables indices in a set S k such that α = α j |j ∈ S k .

(114)

On-line calculations

The on-line calculations refer to the reconstruction of the effective elastic tensor C(t) from the sole knowledge of the internal variables evolution α(t) without performing any simulation on the RVE. The steps of the procedure are described as follows.

• LOOP over all time steps t k • Given α(t k )1 :

• Compute u, v, V and Z using ( 38)- [START_REF] Lubarda | Damage tensors and the crack density distribution[END_REF]. This requires Defining D( α) and D( α) by ( 83)-(87). Computing d(θ, t), h(θ, t) by ( 81)-(82).

Computing η(θ, t), κ(θ, t) by ( 74)-(75). • Compute λ, µ and A by [START_REF] Simo | Strain-and stress-based continuum damage models: I. formulation[END_REF].

• Compute C(t k ) by ( 47)- [START_REF] Carrara | Data-driven fracture mechanics[END_REF].

• END LOOP

In the above, the integrals on S are evaluated numerically. 6 Numerical examples

2D microstructures

First, 2D RVEs are considered, as depicted in Fig. 2. All RVEs have the same dimensions L × L = 1 × 1 mm 2 . The first RVE (RVE A) contains a weak layer oriented along the x -axis, centered in a square domain. The width of the layer is h = 0.05 mm. The second RVE contains a circular pore of radius r = 0.2 mm, which is centered into a square domain. The RVE C contains an inclined weak layer oriented by an angle equal to π/4 with respect to the main frame and of width h = 0.1 mm. Finally, the RVE D corresponds to a realistic microstructure of concrete material, whose geometry is given in Fig. 2 (d). RVEs A and C have been chosen as simple configurations, as they preserve the symmetry of the effective elastic tensor under uniaxial tension. However, RVE C induces a more pronounced anisotropic damage rotated with respect to the main axis. Finally, RVE D induces a complex anisotropic damage, whose evolution is difficult to predict. This case is presented to illustrate the potential of the method for image-based modeling.

The mechanical properties of the different phases are indicated in Table 1. For all RVEs, a mesh composed of 100 × 100 regular 4-node elements has been used.

Inclusion

Matrix depicted for the different RVEs. These graphs indicate that the most significant internal variables, considering all loading paths, are in all cases the variables α 1 , α 2 , α 3 and α 4 , even though the ordering of these internal variables is changed for RVE C. Then, the criterion (112) can be applied to restrict the number of internal variables. For example, using T OL = 10 -1 leads to N = 4 internal variables in all cases.

RVE E (GPa) ν σ c (Gpa) E (GPa) ν σ c (Gpa) A,C 10
The approach is then validated by testing the accuracy of the reconstructed elastic tensor C along a path not included in the off-line calculations. The loading is defined by prescribing boundary conditions [START_REF] Lee | An anisotropic damage criterion for deformation instability and its application to forming limit analysis of metal plates[END_REF] with

ε p (t) = ε(t)n V ⊗n V ,
where n V is the loading direction, characterized by an angle θ V = π/3 (see the corresponding final crack pattern in Fig. 4 (b)). As a reference solution, the effective elastic tensor C ref (t) is computed for each loading step during the crack evolution. The corresponding internal variables α(t) are then extracted by using the procedure described in section 5.1. Then, the approximated value of C(N ) as a function of the number of internal variables N is evaluated by the proposed model, as described in section 5. agreement with the reference solution.

To better appreciate the convergence of the approximation [START_REF] Carrara | Data-driven fracture mechanics[END_REF] with respect to the number of internal variables, we compute the error defined as

Err = C(N ) -C ref C ref (115)
where C(N ) and C ref are evaluated at the last step of the simulation and . is the Frobenius norm. Results are provided in Fig. 11. We can note that for N = 4, which is the number of internal variables selected by the criterion (112) with T OL = 0.1, an error of roughly 10 -2 is found for all RVEs except for RVE C for which the error is slightly higher.

Finally, we have checked the influence of the number P of loading paths on the accuracy of the solution for the last case (RVE D). Precisely, the number of loading paths has an influence on the number of selected internal variables using the criterion (112). For RVE D and a tolerance T ol = 10 -1 , a single loading path P = 1 led to N = 3, which correspond to the results in Fig. 10 (second figure, N = 3), or Fig. 11 for N = 3. For P ≥ 2, the number of internal variables was detected as N = 4 which correspond to the results in Fig. 10 (third figure, N = 4), or Fig. 11 for N = 3. It can be concluded that in this example, a single loading path seems insufficient to accurately capture the influence of all internal variables, while two loading paths are enough to obtain a very satisfying accuracy. Then, even though such discussion is casedependent, it seems that several loading paths are required to capture the influences of the different internal variables for other general loads.

3D microstructures

In the section, the potential of harmonic analysis-based damage models is investigated to describe complex anisotropic damage in the 3D situation. Four RVEs are considered, as depicted in Fig. 12, and they are essentially 3D extensions of the 2D RVEs presented in the previous section. All RVEs have dimensions 1 × 1 × 1 mm 3 . The RVE A concerns a cubic domain containing a weak layer of width h = 0.04 mm and oriented along the x -y plane. The RVE B contains a spherical pore of radius r = 0.2 mm, centered in the cubic domain. The RVE C contains a weak plane which is misoriented with respect to the main axis. The equation of the mid-plane for the layer is z = 0.3x+0.3y+0.25. The width of the layer is h = 0.05 mm. Finally, RVE D is obtained from a micro tomography image of a 3D printed polymer-glass fiber composite, obtained at Navier Laboratory (courtesy M. Bornert, T.X. Le). This microstructure contains fibers which are randomly oriented, but with an average direction in the y-direction, and roughly contained within the x -y plane. Such microstructure induces a strong anisotropy, whose type cannot be easily predicted by only observing the microstructure. RVEs A, B and C are discretized by regular 50 × 50 × 50 8-node cubic finite elements, while RVE D are discretized by 697,659 tetrahedral elements. The properties of the different phases are specified in Table 2. Here, to avoid too large computational simulation times, the preliminary (offline) calculations consist in crack simulations along P = 3 different loading paths defined by prescribing boundary conditions [START_REF] Lee | An anisotropic damage criterion for deformation instability and its application to forming limit analysis of metal plates[END_REF] with ε p (t) = ε(t)n p ⊗n p , where n p is the loading direction characterized by the Euler angles θ p and ϕ p (see Fig. 19), with {θ 1 = 0, ϕ 1 = 0}, {θ 2 = π/2, ϕ 2 = 0} and {θ 3 = π/2, ϕ 3 = π/2}. For each load, 100 increments with ∆ε = 10 -5 have been used. Here, for the sake of simplicity, only linear boundary conditions u(x) = εx = 0 have been employed. As an illustration, the orientation damage distribution functions in polar plots are provided in Figs. by prescribing boundary conditions [START_REF] Lee | An anisotropic damage criterion for deformation instability and its application to forming limit analysis of metal plates[END_REF] with ε p (t) = ε(t)n V ⊗ n V where n V is the loading direction characterized by the angles θ V = π/4 and ϕ V = π/4. The corresponding internal variables α(t) are then extracted using the procedure described in section 5.1. The approximation of C(N ) as a function of the number of internal variables N is evaluated by the proposed model, as described in section 5.2.

Results are provided in Figs. 14, 15, 16 and 17 for RVEs A,B,C and D, respectively. The following observations are provided. First, the number N = 6 associated with the error criterion (112) with T OL = 10 -1 leads in all cases to a very satisfying approximation of all the effective tensor components for RVEs A,B and C. For RVE D, N = 10 is sufficient to obtaining a very good accuracy. Finally, the solutions seem to converge around N = 10 in all cases, which is much lower than the total number N = 21. This illustrates the potential of harmonic analysis-based damage models combined with the present data- driven approach in describing anisotropic damage with a minimal number of internal variables for arbitrary 3D microstructures with a low approximation error.

Finally, we compute the error (115) as a function of N in Fig. 18. We note that for RVE A and B, a convergence is quickly reached for a number of internal variables around N = 9, while for RVEs C and D which exhibit more complex anisotropic behaviors, a slower convergence with respect to N is noticed.

Conclusion

An approach combining harmonic analysis-based damage models and a datadriven procedure has been proposed to construct models for the anisotropic damage of quasi-brittle materials by performing numerical calculations on representative volume elements (RVEs). This approach, based on the relevant elasticity and damage orientation distribution functions and their expansions into Fourier's series, avoids resorting to the usual empirical assumptions made in damage mechanics. It allows extracting damage internal variables in an optimal way according to the degree of accuracy desired, and directly gives damaged macroscopic constitutive laws involving the chosen internal variables.

In addition, combined with the phase-field method, our approach is capable of dealing with complex microstructures, such as those obtained by microcomputed tomography (micro-CT), and remain valid for non-proportional loadings. The nucleation and propagation of microcracks as well as their effects on the effective (or homogenized) behavior can also be taken into account.

In combination with the phase-field method, the approach elaborated in the present work also renders it possible to establish data-driven evolution laws for the chosen damage internal variables within the framework of the thermodynamics of irreversible processes. This will be done and reported in forthcoming studies. 8 Appenfix A: 3D components of F

In this appendix we provide the expressions for the different terms of the fourth-order spectral harmonic function F.

F 1111 =

Figure 1 .

 1 Figure 1. (a) RVE; (b) orientation function η(θ); (c) orientation function κ(θ).

  )-(75) in the present work. At each time t, these quantities can be expanded into Fourier series in the forms d(θ, t) = d 0 (t) + D(t) : F(θ) + D(t) :: F(θ) + ...

h 0 , α 3 =

 3 D 11 , α 4 = D 12 , α 5 = D 1111 , α 6 = D 1112 .(88)

  ϕ)F(θ, ϕ)sin(θ)dθdϕ,

Figure 2 .

 2 Figure 2. 2D geometries of the considered RVEs: A) layered mater; B) porous material; C) layered material with π/4 angle orientation of the layers; D) image-based model of concrete microstructure.

Figure 4 .

 4 Figure 4. Cracks at final load for the different loading paths: (a) loading paths along directions characterized by an angle θ used for data sampling; (b) verification path.

2 .Figure 5 .

 25 Figure 5. Evolution of internal variables for the 2D RVEs, for load path characterizd by a loading direction θ = π/2.

Figure 6 .

 6 Figure 6. Relative weights of internal variables for the 2D RVEs.

N=5Figure 7 .Figure 8 .

 78 Figure 7. 2D RVE A: comparison between direct numerical simulation (reference) and approximation of the effective elasticity tensor using N internal variables.

N=6Figure 9 .N=6Figure 10 .

 910 Figure 9. 2D RVE C: comparison between direct numerical simulation (reference) and approximation of the effective elasticity tensor using N internal variables.

Figure 11 .

 11 Figure 11. Error of the approximated effective elastic tensor at final step of the loading history with respect to the number of internal variables for the 2D RVEs.

Figure 12 .

 12 Figure 12. 3D RVEs: geometries; A: Weak layer; B) porous RVE; C) Weak layer with mismatch orientation with respect to the main frame; D) RVE obtained from a micro-tomography image of a 3D printed glass-fiber polymer.

Figure 13 .

 13 Figure 13. Relative weights of internal variables for the 3D RVEs.

Figure 14 .Figure 15 .

 1415 Figure 14. 3D RVE A: comparison between direct numerical simulation (reference) and approximation of the effective elasticity tensor using N internal variables.

Figure 16 .Figure 17 .

 1617 Figure 16. 3D RVE C: comparison between direct numerical simulation (reference) and approximation of the effective elasticity tensor using N internal variables.

Figure 18 .

 18 Figure 18. Error of the approximated effective elastic tensor at final step of the loading history with respect to the number of internal variables for the 3D RVEs.

Table 2 :

 2 Material parameters for the 3D RVEs.

		Inclusion			Matrix
	RVE E (GPa) ν σ c (Gpa) E (GPa) ν σ c (Gpa)
	A,D	10.4	0.3	0.01	52	0.3	0.03
	B	10 -3	0	10 -4	52	0.3	0.03
	C	10 4	0.3	10 4	52	0.3	0.03

  In this appendix, we provide plots of the orientations functions d and h in both 2D and 3D cases, in Figs.20 and 21, respectively.

								z
								q	n
	n 4 1 -	6 7	n 2 1 +	3 35	(116)
	F 1122 = n 2 1 n 2 2 -	1 7	(n 2 1 + n 2 2 ) +	1 35	j	y	(117)
	F 1133 = n 2 1 n 2 3 -	1 7	(n 2 1 + n 2 3 ) + x	1 35	(118)
	F 1112 = n 3 1 n 2 -	3 7	n 1 n 2	Figure 19. Sperical frame.	(119)
	F 1113 = n 3 1 n 3 -F 1212 = n 2 1 n 2 2 -	3 7 1 7	n 1 n 3 (n 2 1 + n 2 2 ) +	1 35	(120) (131)
	F 1123 = n 2 1 n 2 n 3 -F 1213 = n 2 1 n 2 n 3 -F 2222 = n 4 2 -6 7 n 2 2 + 1 n 2 n 3 1 n 2 n 3 7 7 3 F 1223 = n 1 n 2 2 n 3 -1 n 1 n 3 7 35 F 2233 = n 2 2 n 2 3 -1 7 (n 2 2 + n 2 3 ) + F 1313 = n 2 1 n 2 3 -1 7 (n 2 1 + n 2 3 ) + F 2212 = n 1 n 3 2 -3 n 1 n 2 F 1323 = n 1 n 2 n 2 3 -1 n 1 n 2 7 7 F 2213 = n 1 n 2 2 n 3 -7 n 1 n 3 1 F 2323 = n 2 2 n 2 3 -1 7 (n 2 2 + n 2 3 ) +	1 1 35 35 1 35	(132) (121) (133) (122) (134) (123) (135) (124) (125) (136)
	F 2223 = n 3 2 n 3 -9 Appendix B: Plots of orientation distribution functions 3 7 n 2 n 3	(126)
	F 3333 = n 4 3 -	6 7	n 2 3 +	3 35	(127)
	F 3312 = n 1 n 2 n 2 3 -	1 7	n 1 n 2	(128)
	F 3313 = n 1 n 3 3 -	3 7	n 1 n 3	(129)
	F 3323 = n 2 n 3 3 -	3 7	n 2 n 3	(130)

To avoid confusion, we specify that in the present work, we do not provide evolution equations for the internal variables. Then, the reconstruction of C(t k ) requires the knowledge of the evolution of the internal variables for a given loading path, obtained by full-field simulations. Of course, a further objective will be to obtain a model which does not require the use of full-field simulations in the on-line stage. This point is the topic of ongoing investigations.

Acknowledgements

We thank Michel Bornert and Thi Xiu Le for providing the segmented image used to construct the 3D microstructural mesh of the realistic 3D image-based microstructure.

The internal variables are then extracted from the procedure described in section 5.1. As an illustration, the damage orientation distribution functions at final loading step n = 80 for RVEs A and C are depicted in Appendix 9, Fig. 20.

The evolutions of the internal variables for the different RVEs and for the loadinng with θ = π/2 are depicted in Fig. 5.

In Fig. 6, the evolution of the relative cumulated internal variables (113) is