# Anisotropic damage state modeling based on harmonic decomposition and discrete simulation of fracture

Flavien Loiseau, Cécile Oliver Oliver-Leblond, Thomas Verbeke, Rodrigue

Desmorat

► To cite this version:

Flavien Loiseau, Cécile Oliver Oliver-Leblond, Thomas Verbeke, Rodrigue Desmorat. Anisotropic damage state modeling based on harmonic decomposition and discrete simulation of fracture. Engineering Fracture Mechanics, 2023, 293, pp.109669. 10.1016/j.engfracmech.2023.109669. hal-04245338

HAL Id: hal-04245338

https://hal.science/hal-04245338v1

Submitted on 17 Oct 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Highlights

Anisotropic damage state modeling based on harmonic decomposition and discrete simulation of fracture

F. Loiseau, C. Oliver-Leblond, T. Verbeke, R. Desmorat

* This study proposes an anisotropic damage state modeling for quasi-brittle materials.

* It aims at representing around 76 000 (discrete) beam-particle computations of bi-dimensional effective (damaged) elasticity tensors.

* It proposes a methodology to formulate an anisotropic damage model based on a covariant reconstruction formula of bi-dimensional orthotropic elasticity tensors.
# Anisotropic damage state modeling based on harmonic decomposition and discrete simulation of fracture

F. Loiseau<sup>a,*</sup>, C. Oliver-Leblond<sup>a</sup>, T. Verbeke<sup>a</sup> and R. Desmorat<sup>a</sup>

<sup>a</sup>Université Paris-Saclay, CentraleSupélec, ENS Paris-Saclay, CNRS, LMPS - Laboratoire de Mécanique Paris-Saclay, 91190, Gif-sur-Yvette, France.

## ARTICLE INFO

**Keywords:**

Anisotropic damage

Discrete Element method

Effective elasticity tensor

Harmonic decomposition

Microcracking

## ABSTRACT

This study proposes an anisotropic damage state modeling based on (i) a discrete element model for quasi-brittle material and (ii) a decomposition of the elasticity tensor in covariants. A procedure is proposed to measure the evolutions of effective (damaged) elasticity tensors computed by a beam-particle model. Various multiaxial damaging loadings allow us to constitute a dataset of around 76 000 effective elasticity tensors. We then identify, in a tensorial manner, the anisotropic damage state for the whole dataset. A detailed analysis of the dataset, using the distance to orthotropy as a guideline, justifies representing the induced micro-cracking by a single second-order damage variable, even in the final stages with strong micro-cracks interaction. To formulate the damage state coupling, we use a reconstruction formula of orthotropic elasticity tensors in terms of invariants and (tensor) covariants. Thanks to this formula, some parts of the effective elasticity tensors $\tilde{E}$ (such as the dilatation part) are modeled exactly from the single damage variable. Constitutive equations are proposed for the remaining parts of $\tilde{E}$ (such as its generalized shear modulus and fourth-order harmonic part) using physical assumptions from micro-mechanics and a sparse data driven approach. The proposed anisotropic damage state coupling accurately models the damaged elasticity tensors in multiaxial loading, proportional or non-proportional, up to high damages. The present study firstly highlights the need for an anisotropic damage model for quasi-brittle materials and, secondly, offers a methodology to formulate the damage state coupling by explicit formulas introducing at most two dedicated parameters: the (optional) nonlinear shear-damage coupling parameter *m* and the harmonic prefactor *h*.

## Introduction

Predicting the behavior of quasi-brittle materials – such as concrete or mortar – is essential to guarantee the integrity of civil engineering structures. When undamaged, these quasi-brittle materials often exhibit an isotropic elastic behavior. Severe mechanical loading leads to the nucleation and growth of micro-cracks and the loss of load-bearing capabilities (Bažant and Gambarova, 1984; Bažant and Oh, 1985; Landis, 1999). The orientation of the micro-cracks due to loading produces induced mechanical anisotropy (Mazars, 1984; Bažant and Prat, 1988a,b; Mazars et al., 1990; Ramtani et al., 1992; Lubarda and Krajcinovic, 1993; Fichant et al., 1999). Further increase of the loading gradually leads to the concentration and the coalescence of micro-cracks, which result in structural failure by macroscopic cracking (Lemaitre, 1992).

Continuous damage models account for the material mechanical degradation during loading, including damage initiation (Lemaitre, 1992; Krajcinovic, 1996; Lemaitre and Desmorat, 2005; Murakami, 2012). The modeling of macroscopic damage usually starts with the choice of a thermodynamics variable, the damage variable, which represents the micro-cracking state of the material. The tensorial nature of the damage variable has been discussed in classical literature (Vakulenko and Kachanov, 1971; Chaboche, 1979; Leckie and Onat, 1981; Chaboche, 1984; Lemaitre and Chaboche, 1985; Murakami, 1988) as well as in recent works (Cormery and Welemane, 2010; Desmorat and Desmorat, 2016; Dormieux and Kondo, 2016; Fassin et al., 2019; Oliver-Leblond et al., 2021). The main question was whether or not the damage state could be represented by a fourth-order tensor (Chaboche, 1978, 1979; Krajcinovic, 1985; Kachanov, 1993), by two second-order tensors (Ladevèze, 1983, 1995; Desmorat and Desmorat, 2016), or by a single second-order tensor (Murakami and Ohno, 1978; Cordebois and Sidoroff, 1980, 1982; Murakami, 1988;

<figure><img src="image_2.png" /></figure>

<footnote>*Corresponding author.</footnote>

<footnote>flavien.loiseau@ens-paris-saclay.fr (F. Loiseau); cecile.oliver@ens-paris-saclay.fr (C. Oliver-Leblond); thomas.verbeke@ens-paris-saclay.fr (T. Verbeke); rodrigue.desmorat@ens-paris-saclay.fr (R. Desmorat) ORCID(s): 0000-0002-6305-7309 (F. Loiseau); 0000-0002-8656-6800 (C. Oliver-Leblond); 0000-0003-3934-5738 (R. Desmorat)</footnote>
Kachanov, 1993; Halm and Dragon, 1996; Papa and Taliercio, 1996; Lemaitre et al., 2000; Lemaitre and Desmorat, 2005; Desmorat et al., 2018) instead of a scalar variable.

The simplest modeling choice is the one of a scalar damage variable (Kachanov, 1958; Rabotnov, 1969; Lemaitre, 1984). Thanks to their simplicity, isotropic damage models are often used to compute the degradation of concrete structures (Mazars, 1984; Grassl and Jirásek, 2006; Richard and Ragueneau, 2013). These models assume that the material behavior remains isotropic when damaging. In accordance with damage measurements (Ramtani et al., 1992; Lemaitre et al., 2000), micro-mechanics studies of micro-cracked media, such as the ones of Chaboche (1984), Lubarda and Krajcinovic (1993) and Kachanov (1993), show indeed that damage is not isotropic and has to be represented by a damage tensor physically linked to the crack density:

* a fourth-order tensor in the 3D case,

* a second-order tensor in the simpler 2D case with lubricated non-interacting cracks.

We aim to generalize the latter 2D result to the case of strongly interacting micro-cracks, i.e., up to their coalescence, to model the total failure of a specimen (an Area Element in the present work). A first possibility would be to use a nonlinear homogenization scheme (Kachanov, 1993; Ponte Castañeda and Willis, 1995; Dormieux and Kondo, 2016). But these schemes rely on Fracture Mechanics at the microscale. This theory that does not deal with crack initiation, nor does accurately represent multiple cracks interaction and coalescence. Following Rinaldi and Lai (2007) and Delaplace and Desmorat (2008), we prefer instead to rely on discrete simulations of multiple cracking by lattice models (Hrennikoff, 1941; Kawai, 1978; Herrmann and Roux, 1990; Schlangen and van Mier, 1992; Bolander et al., 1996) — more precisely, on a beam-particle model (Delaplace, 2008; Vassaux et al., 2016). Indeed, these discrete models are based on the brittle failure of beams (at the micro-scale). They allow for both the modeling of micro-cracks initiation and the natural representation of micro-cracks coalescence. They represent well the material behavior of quasi-brittle materials such as concrete or mortar (Bažant et al., 1990; Delaplace et al., 1996; van Mier et al., 2002; Challamel et al., 2015; Oliver-Leblond, 2019).

We aim at formulating an accurate anisotropic damage state coupling for quasi-brittle materials as a twin modeling of discrete (beam-particle) bi-dimensional fracture. The present work is a continuation of the one of Oliver-Leblond et al. (2021), who did introduce a suitable macroscopic damage variable with such an approach. But, contrary to these authors, who did assess for the representativity of their damage variable in a few loading cases only, we here build a large dataset of around 76 000 bi-dimensional effective elasticity tensors. The dataset is constituted of tensors from 21 mechanical loadings, proportional or not, and 36 virtual specimens with different micro-structures. Repeated computations for different micro-structures allow us to obtain statistically representative results. Performing complex multiaxial loadings allow to obtain multiple fracture patterns representative of quasi-brittle specimens. The harmonic decomposition of bi-dimensional elasticity tensors (Blinowski et al., 1996; Vianello, 1997) is completed by an analysis their harmonic part (Desmorat and Desmorat, 2015) and by computations of their distance to orthotropy. It is used to derive a meaningful macroscopic second-order damage variable, as well as its full coupling with the scalar (invariant) and tensorial (covariant) components of the elasticity tensor harmonic decomposition.

Section 1 introduces the harmonic decomposition and the reconstruction of a bi-dimensional elasticity tensor by means of its covariants, which serves as the basis of the proposed modeling. It also defines the distance of an elasticity tensor to a symmetry class, which will be used to analyze our whole dataset. Section 2 begins with the presentation of the considered discrete model, then details a measurement method for the effective elasticity tensors computed by discrete virtual testing (of Area Element cracked by severe mechanical loading). This procedure is applied to multiple complex mechanical loadings, proportional or not, to constitute a dataset of 76 356 effective (i.e., damaged) elasticity tensors. Based on this large dataset, sections 3 and 4 are dedicated to modeling the effect of micro-cracking, i.e., of anisotropic damage, on the (invariant) shear modulus $\mu$ and the (covariant) fourth-order harmonic part $\mathbf{H}$ of the elasticity tensors. Section 5 summarizes and combines the results from sections 1, 3, and 4 to propose an anisotropic damage state coupling by means of only two specific material parameters: the nonlinear shear-damage coupling parameter $m$ and the harmonic prefactor $h$. Finally, the proposed coupling is assessed in section 6.

## Notations and definitions

Let $d = 2$ be the dimension, and $O(2)$ be the orthogonal group. The (left) action of an orthogonal transformation $Q \in O(2)$ on a second-order tensor $\mathbf{t}$ or a fourth-order tensor $\mathbf{T}$ is

$$ (Q \star \mathbf{t})_{ij} = Q_{ik}Q_{jl}t_{kl}, \quad (Q \star \mathbf{T})_{ijkl} = Q_{ip}Q_{jq}Q_{kr}Q_{ls}T_{pqrs}. \quad (1) $$
An elasticity tensor is a positive-definite fourth-order tensor **E** having the index symmetries $E_{ijkl} = E_{jikl} = E_{ijkl} = E_{klij}$. The vector space

$$ \mathbb{E}\mathrm{la}(\mathbb{R}^2) = \{\mathbf{E} \mid E_{ijkl} = E_{jikl} = E_{ijkl} = E_{klij}\}, \qquad (2) $$

of bi-dimensional elasticity tensors, is of dimension 6. A covariant of a bi-dimensional elasticity tensor **E** is a tensorial function **C**(**E**) such that

$$ \mathbf{C}(\mathbf{Q} \star \mathbf{E}) = \mathbf{Q} \star \mathbf{C}(\mathbf{E}), \quad \forall \mathbf{Q} \in \mathcal{O}(2). \qquad (3) $$

An invariant of a *d*-dimensional elasticity tensor **E** is a covariant of order zero of **E**, i.e., a function $I(\mathbf{E})$ such that

$$ I(\mathbf{Q} \star \mathbf{E}) = I(\mathbf{E}), \quad \forall \mathbf{Q} \in \mathcal{O}(2). \qquad (4) $$

A harmonic tensor is a traceless totally symmetric tensor. A second-order harmonic tensor is a so-called deviatoric tensor. We denote by $\mathbb{H}^n(\mathbb{R}^2)$ the vector space of harmonic bi-dimensional tensors of order $n$. Remark that $\mathbb{H}^0(\mathbb{R}^2)$ is isomorphic to $\mathbb{R}$.

# 1. Reconstruction of a bi-dimensional orthotropic elasticity tensor by means of its covariants

## 1.1. Harmonic decomposition

The harmonic decomposition of a 2D elasticity tensor **E** is its equivariant decomposition into harmonic tensors (Blinowski et al., 1996; Vianello, 1997)

$$ \mathbf{E} = (\mu, \kappa, \mathbf{d}', \mathbf{H}), \qquad (5) $$

such that

$$ \mathbf{Q} \star \mathbf{E} = (\mu, \kappa, \mathbf{Q} \star \mathbf{d}', \mathbf{Q} \star \mathbf{H}), \quad \forall \mathbf{Q} \in \mathcal{O}(2), \qquad (6) $$

the harmonic components of **E** being the invariants $\mu, \kappa \in \mathbb{H}^0(\mathbb{R}^2)$ (the generalized shear and bulk moduli, respectively), the second-order covariant $\mathbf{d}' = \mathbf{d}'(\mathbf{E}) \in \mathbb{H}^2(\mathbb{R}^2)$ and the fourth-order covariant $\mathbf{H} = \mathbf{H}(\mathbf{E}) \in \mathbb{H}^4(\mathbb{R}^2)$. The covariant $\mathbf{d}'(\mathbf{E})$ is the deviatoric part $\mathbf{d}' = \mathbf{d} - \frac{1}{2}(\mathrm{tr} \, \mathbf{d}) \mathbf{1}$ of the dilatation tensor

$$ \mathbf{d} = \mathrm{tr}_{34} \mathbf{E} = \mathbf{E}:\mathbf{1}, \qquad (7) $$

and the covariant $\mathbf{H}(\mathbf{E})$ is the fourth-order harmonic part of **E**. We note **1** the second-order identity tensor. The scalar components of the harmonic decomposition of **E** are

$$ \mu = \frac{1}{8} (2 \mathrm{tr} \, \mathbf{v} - \mathrm{tr} \, \mathbf{d}), \quad \kappa = \frac{1}{4} \mathrm{tr} \, \mathbf{d}, \qquad (8) $$

where $\mathbf{v} = \mathbf{v}(\mathbf{E})$ is the second-order covariant

$$ \mathbf{v} = \mathrm{tr}_{13} \mathbf{E}. \qquad (9) $$

It is such that $\mathbf{v}' = \mathbf{d}'$.

The scalars

$$ I_2(\mathbf{E}) = I_2(\mathbf{d}') = \|\mathbf{d}'\|^2 = \mathbf{d}':\mathbf{d}', \quad \|\mathbf{H}\|^2 = \mathbf{H}:\mathbf{H}, \quad \text{and} \quad K_3(\mathbf{E}) = \mathbf{d}:\mathbf{H}:\mathbf{d}, \qquad (10) $$

are three other invariants of the elasticity tensor (Vianello, 1997).

A first reconstruction formula of an elasticity tensor **E** by means of its covariant is the explicit harmonic decomposition itself,

$$ \mathbf{E} = 2\mu\mathbf{J} + \kappa\mathbf{1} \otimes \mathbf{1} + \frac{1}{2}(\mathbf{1} \otimes \mathbf{d}' + \mathbf{d}' \otimes \mathbf{1}) + \mathbf{H}, \qquad \left\{ \begin{array}{l} \mathbf{Iso} = 2\mu\mathbf{J} + \kappa\mathbf{1} \otimes \mathbf{1}, \\ \mathbf{Dil} = \frac{1}{2}(\mathbf{1} \otimes \mathbf{d}' + \mathbf{d}' \otimes \mathbf{1}), \end{array} \right. \qquad (11) $$

which defines **Iso** as the isotropic part of **E**, **Dil** as its dilatation part and $\mathbf{H} = \mathbf{E} - \mathbf{Iso} - \mathbf{Dil}$ as its fourth-order harmonic part. All three tensors **Iso**, **Dil** and **H** are fourth-order covariants of **E**. The reconstruction formula (11) applies to all 2D elasticity tensors, possibly fully anisotropic (i.e., biclinic). In Eq. (11), **I** is the fourth order identity tensor (of components $I_{ijkl} = \frac{1}{2}(\delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk})$) and $\mathbf{J} = \mathbf{I} - \frac{1}{2}\mathbf{1} \otimes \mathbf{1}$ is the fourth-order deviatoric projector.
## 1.2. Harmonic square

To obtain a reconstruction formula dedicated to orthotropic elasticity tensors, we have to introduce the harmonic product $\mathbf{h}_1 * \mathbf{h}_2 \in \mathbb{H}^4(\mathbb{R}^2)$ of two second-orders harmonic tensors $\mathbf{h}_1$, and $\mathbf{h}_2$. In 2D, it is defined as the fourth-order harmonic tensor (see Olive et al. (2017))

$$ \mathbf{h}_1 * \mathbf{h}_2 = \mathbf{h}_1 \otimes \mathbf{h}_2 - \frac{1}{2}(\mathbf{h}_1 : \mathbf{h}_2) \mathbf{J}, \quad \mathbf{J} = \mathbf{I} - \frac{1}{2} \mathbf{I} \otimes \mathbf{I}. \qquad (12) $$

where $\mathbf{J}$ is the so-called deviatoric projector. Recall then that any 2D fourth-order harmonic tensor $\mathbf{H} \in \mathbb{H}^4(\mathbb{R}^2)$ can be written as a harmonic square (Desmorat and Desmorat, 2015),

$$ \mathbf{H} = 2\Lambda \mathbf{e} * \mathbf{e}, \quad \mathrm{tr} \, \mathbf{e} = 0, \quad \|\mathbf{e}\| = 1, \qquad (13) $$

where $\mathbf{e} \in \mathbb{H}^2(\mathbb{R}^2)$ is a unit second-order deviatoric eigentensor associated with a non-zero eigenvalue $\Lambda$ or $-\Lambda$ of the harmonic tensor $\mathbf{H}$. This is not a so-called reconstruction formula for $\mathbf{H}$ since $\mathbf{e}$ is not a covariant of $\mathbf{H}$.

A more interesting formula is obtained when $\mathbf{H} = \mathbf{H}(\mathbf{E})$ is the fourth order harmonic part of an (exactly) orthotropic elasticity tensor. We have in that case, in 2D still (Oliver-Leblond et al., 2021),

$$ \mathbf{H} = \frac{2K_3(\mathbf{E})}{I_2^2(\mathbf{E})} \mathbf{d}'(\mathbf{E}) * \mathbf{d}'(\mathbf{E}). \qquad (14) $$

Remark that the covariants—including invariants—involved are covariants of the elasticity tensor, which is orthotropic, not of the harmonic tensor $\mathbf{H}$, which has the square symmetry (Verchery, 1982; Vianello, 1997; Vannucci, 2005). This means that we have the following reconstruction formula (by means of its covariants) for a 2D orthotropic elasticity tensor

$$ \mathbf{E} = 2\mu\mathbf{J} + \kappa\mathbf{1} \otimes \mathbf{1} + \frac{1}{2}(\mathbf{1} \otimes \mathbf{d}' + \mathbf{d}' \otimes \mathbf{1}) + \frac{2K_3}{I_2^2} \mathbf{d}' * \mathbf{d}', \qquad (15) $$

where $\mu$, $\kappa$, $I_2$ and $K_3$ are invariants of the elasticity tensor (defined by Eq. (8) and Eq. (10)), and $\mathbf{d}$ is a second-order covariant of $\mathbf{E}$ (defined by Eq. (7)). Remark that for this formula to hold, since $I_2(\mathbf{d}' = 0) = 0$ it is necessary that $\mathbf{d}' \neq 0$. In other words, it is necessary that the dilatation tensor $\mathbf{d} = \mathbf{d}(\mathbf{E})$ —which inherits the symmetry of $\mathbf{E}$, Olive et al. (2022)— and therefore $\mathbf{E}$, are orthotropic.

## 1.3. Definition of a tensorial damage variable

We can furthermore assume that the elasticity tensor $\mathbf{E}$ of a quasi-brittle material evolves during loading, for instance, due to damage, and that it has the initial isotropic value

$$ \mathbf{E}_0 = 2\mu_0\mathbf{J} + \kappa_0\mathbf{1} \otimes \mathbf{1}. \qquad (16) $$

The initial dilatation tensor is then isotropic (i.e., spherical),

$$ \mathbf{d}_0 = \mathbf{E}_0 : \mathbf{1} = 2\kappa_0 \mathbf{1}. \qquad (17) $$

A damage variable, noted $\mathbf{D}$ in the present work, represents the state of micro-cracking of a quasi-brittle material (Lemaitre and Chaboche, 1985). It is set as zero when the effective elasticity tensor $\mathbf{E}$ is the initial elasticity tensor $\mathbf{E}_0$. The eigenvalues of the damage tensor are usually bounded by 1. The present work aims to determine the general coupling with damage $\mathbf{E} = \mathbf{E}(\mathbf{D})$.

By Eq. (8), the bulk modulus $\kappa$ is exactly reconstructed from the dilatation tensor $\mathbf{d}$. This has led Oliver-Leblond et al. (2021) to define the dimensionless damage variable $\mathbf{D}$ as the second-order tensor

$$ \mathbf{D} = \mathbf{1} - \frac{\mathbf{d}}{2\kappa_0}. \qquad (18) $$

Remark that since the initial dilatation tensor is given by $\mathbf{d}_0 = \mathrm{tr}_{12} \mathbf{E}_0 = 2\kappa_0 \mathbf{1}$, we have the equalities

$$ \mathbf{D} = \frac{1}{2\kappa_0} \mathrm{tr}_{12}(\mathbf{E}_0 - \mathbf{E}) = (\mathbf{d}_0 - \mathbf{d}) \cdot \mathbf{d}_0^{-1} = \mathbf{d}_0^{-1} \cdot (\mathbf{d}_0 - \mathbf{d}). \qquad (19) $$
This mapping provides a bijection between the damage variable **D** and the dilatation tensor **d**. Indeed, the latter and the bulk modulus are related to **D** as

$$ \mathbf{d} = 2\kappa_0 (\mathbf{1} - \mathbf{D}), \quad \kappa = \frac{1}{4} \operatorname{tr} \mathbf{d} = \kappa_0 \left(1 - \frac{1}{4} \operatorname{tr} \mathbf{D}\right). \qquad (20) $$

## 1.4. Distance to isotropy – Distance to orthotropy

A key question is then whether or not the damage of quasi-brittle materials can be represented in 2D by the single second-order damage tensor **D**. An underlying question is whether or not the damaged elasticity tensor of an Area Element (AE) of a quasi-brittle material can be close to orthotropy (since **D** is either isotropic or orthotropic). One will also have to check to what extent the initial (elastic) AE is isotropic.

The determination of the symmetry class of a measured elasticity tensor is a difficult problem (Gazis et al., 1963; François, 1995; François et al., 1998; Moakher and Norris, 2006; Diner et al., 2011). A cause is that the measurement orientation might not correspond to the principal direction of the expected symmetry class, which prevents direct identification by comparison to normal (Kelvin) forms. Furthermore, experimental measurements provide a noised approximation of the material’s elastic properties (Roux et al., 1985; Migliori et al., 1993). The measured elasticity tensor will generically be biclinic in 2D (and triclinic in 3D).

Those issues can be mitigated by calculating the distance to the expected elasticity symmetry class, more precisely, the distance to the considered symmetry stratum (which is the set of all tensors which have the same symmetry class (Auffray et al., 2014; Abramian et al., 2020)). This usually consists, first, in finding elasticity tensor **E**\* in the symmetry stratum $\Sigma$ which is the closest to the measured elasticity tensor **E**, and second, in calculating the distance between **E** and **E**\*. The relative distance to the symmetry stratum $\Sigma$ is then

$$ \Delta_{\Sigma}(\mathbf{E}) = \min_{\mathbf{E}^* \in \Sigma} \frac{\|\mathbf{E} - \mathbf{E}^*\|}{\|\mathbf{E}\|}. \qquad (21) $$

In the case of isotropy, the harmonic decomposition provides, by orthogonal projection (Vianello, 1997), the closest isotropic tensor **E**\* to **E** as its isotropic part **Iso** (defined by Eq. (11)),

$$ \Delta_{Iso}(\mathbf{E}) = \min_{\mathbf{E}^* \text{ isotropic}} \frac{\|\mathbf{E} - \mathbf{E}^*\|}{\|\mathbf{E}\|} = \frac{\|\mathbf{E} - \mathbf{Iso}(\mathbf{E})\|}{\|\mathbf{E}\|}, \quad \mathbf{E}^* = \mathbf{Iso}(\mathbf{E}). \qquad (22) $$

The calculation of the distance to 2D elastic orthotropy requires more mathematical development (Vianello, 1997; Antonelli et al., 2022). An upper bound of this distance is obtained thanks to the orthotropic reconstruction formula (15) (Oliver-Leblond et al., 2021).

$$ \Delta_{\text{Ort}}(\mathbf{E}) = \min_{\mathbf{E}^* \text{ orthotropic}} \frac{\|\mathbf{E} - \mathbf{E}^*\|}{\|\mathbf{E}\|} \le \frac{\|\mathbf{E} - \mathbf{E}^{\text{up}}\|}{\|\mathbf{E}\|}, \qquad (23) $$

where

$$ \mathbf{E}^{\text{up}} = 2\mu\mathbf{J} + \kappa\mathbf{1} \otimes \mathbf{1} + \frac{1}{2}(\mathbf{1} \otimes \mathbf{d}' + \mathbf{d}' \otimes \mathbf{1}) + \frac{2K_3}{I_2^2} \mathbf{d}' * \mathbf{d}'. \qquad (24) $$

This upper bound, based on the covariants of the elasticity tensor, is easier to calculate than the exact distance (the corresponding formulas are recalled in Appendix A).

## 2. Discrete virtual testing

This part aims at presenting the discrete virtual testing procedure based on discrete simulations of the cracking of Area Elements. The discrete model used is a hybrid beam-particle model which combines two types of models (Meguro and Hakuno, 1989; Kun and Herrmann, 1996; D'Addetta et al., 2002; Bolander et al., 2021): the lattice models (Hrennikoff, 1941; Kawai, 1978; Herrmann and Roux, 1990) and the particular models (Cundall and Strack, 1979; Bažant et al., 1990). The model, which directly originates from the work of (Delaplace et al., 1996; Delaplace, 2008), is first presented, followed by the procedure for the systematic measurement of elasticity tensors. Then, the evolutions of effective elasticity tensors measured during various loadings allow us to generate a large dataset of elasticity tensors.
<figure><img src="image_3.png" /><figcaption>**Figure 1:** Representation of the mesh of the beam-particle model. (1) is the positioning of the particle centers in the grid. (2) is the generation of the particle border by the Voronoi tessellation of the particle centers. (3) adds the beam network based on the Delaunay triangulation of the particle centers.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 2:** Representation of a beam and notations</figcaption></figure>

## 2.1. Beam-particle model

The discrete hybrid model considered is the beam-particle one developed by (Vassaux, 2015; Vassaux et al., 2016). As shown by Oliver-Leblond (2019), it accurately represents the failure process encountered in quasi-brittle materials. For this study, a new linear solver (Davis, 2011) has been incorporated in the implementation. Here, the 2D version of the model is used.

An elementary area is modeled using a beam-particle model. The left part of Figure 1 illustrates the procedure to generate the mesh of the model. The beam-particle specimen is composed of a set of rigid particles. Particle centers are randomly placed in each cell of a grid with a cell size $\bar{l}_b$ which corresponds to the average beam length. Particle boundaries are obtained from the Voronoi tessellation of the particle centers. The dual graph of the Voronoi tessellation is the Delaunay triangulation. It associates a segment to each pair of neighboring particles. Those segments are used as the geometric support for a beam network. This beam network models the cohesion of the material. Each beam $(p, q)$, linking particles $p$ and $q$, has an Euler-Bernoulli behavior and is parametrized by:

* its length $l_{(p,q)}$,

* its section $A_{(p,q)}$,

* its Young modulus $E_{(p,q)}$, and

* its coefficient of inertia $\alpha_{(p,q)} = 64I_{(p,q)}\pi/A_{(p,q)}^2$.
<table><caption>**Table 1**

Parameters of beam-particle model</caption><thead><tr><th>Name</th><th>Symbol</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td>Average beam length</td><td>l&#773;<sub>b</sub></td><td>0.002</td><td>m</td></tr><tr><td>Young's modulus</td><td>E<sub>b</sub></td><td>60</td><td>GPa</td></tr><tr><td>Coefficient of inertia</td><td>&alpha;<sub>b</sub></td><td>0.85</td><td>-</td></tr><tr><td>Scale factor (extension)</td><td>&lambda;<sub>&epsilon;cr</sub></td><td>5.0 &times; 10<sup>-4</sup></td><td>-</td></tr><tr><td>Scale factor (rotation)</td><td>&lambda;<sub>&theta;cr</sub></td><td>2.8 &times; 10<sup>-3</sup></td><td>-</td></tr><tr><td>Shape factor</td><td>k</td><td>1.0</td><td>-</td></tr></tbody></table>

The geometric parameters $l_{(p,q)}$ and $A_{(p,q)}$ are obtained from the mesh's geometry; they thus depend on the beam. Note that the cell size of the grid $\bar{l}_b$ corresponds to the average beam length. The mechanical parameters $E_{(p,q)}$ and $\alpha_{(p,q)}$ are chosen equal to $E_b$ and $\alpha_b$ for all beams, which are identified to match a cement mortar macroscopic elastic behavior.

Fracture properties are introduced by adding a brittle failure criterion $P_{(p,q)}$ to each beam, so beam rupture occurs when

$$P_{(p,q)} = \frac{\varepsilon_{(p,q)}}{\varepsilon_{(p,q)}^{cr}} + \frac{|\theta_{(p)} - \theta_{(q)}|}{\theta_{(p,q)}^{cr}} > 1. \quad (25)$$

where:

* $\varepsilon_{(p,q)} = \frac{\|u_{(p)} - u_{(q)}\|}{l_{(p,q)}}$ is the extension of the beam,

* $u_{(p)}$ is the displacement of particle $p$,

* $\theta_{(p)}$ is the rotation of particle $p$,

* $\varepsilon_{(p,q)}^{cr}$ is the breaking threshold in extension of beam $(p, q)$, and

* $\theta_{(p,q)}^{cr}$ is the breaking threshold in rotation.

Those quantities are illustrated in Figure 2.

Delaplace and Desmorat (2008) took the same failure thresholds for all beams. Following (Rossi and Richer, 1987; de Arcangelis and Herrmann, 1989; Herrmann et al., 1989; Herrmann and Roux, 1990; D'Addetta et al., 2002; Vassaux, 2015), random distributions for the failure thresholds are here considered, as they are more suitable for the modeling of cement (Schlangen and van Mier, 1992; van Mier et al., 2002).

Here, both breaking thresholds $\varepsilon_{(p,q)}^{cr}$ and $\theta_{(p,q)}^{cr}$ are randomly drawn from a Weibull distribution,

$$f(x) = \frac{k}{\lambda} \left(\frac{x}{\lambda}\right)^{k-1} e^{-\left(\frac{x}{\lambda}\right)^k} \quad (26)$$

where $\lambda$ is the scale factor and $k$ is the shape factor. The spatial variability of the breaking thresholds is supposed to be identical for both thresholds: $k_{\varepsilon cr} = k_{\theta cr} = k$. This means that the fracture is controlled by three parameters: the shape factor $k$, the scale factor in extension $\lambda_{\varepsilon cr}$, and the scale factor in rotation $\lambda_{\theta cr}$. Those three parameters are identified by fitting the non-linear macroscopic behavior.

**Remark 1.** The beam-particle model considered can also represent the microcrack closure effects by adding contact and friction between the particles when a beam is broken. Contact and friction are not accounted for in this work.

The parameters of the beam-particle model, given in Table 1, correspond to a quasi-brittle material such as cement with the following properties: a Young modulus $E_0 = 36.35$ GPa, a Poisson ratio $v_0 = 0.22$, a tensile strength $f_t = 5$ MPa. Note that the elastic properties are equivalent to a bulk modulus $\kappa_0 = 30.0$ GPa and a shear modulus $\mu_0 = 19.4$ GPa.
<figure><img src="image_5.png" /><figcaption>Figure 3: Periodic beam-particle mesh where the red boxed particles are guiding particles and green circled ones are guided particles. Guided particles have the same number as their guiding particle.</figcaption></figure>

## 2.2. Periodic Boundary Conditions (PBC)

Periodic Boundary Conditions are used to compute the effective elasticity tensor of a square specimen. The PBC are imposed by adding a layer of guided particles on the top and right part of the square mesh, as shown in Figure 3. Each guided particle has the same geometry as the associated guiding particle on the bottom and left part of the mesh. The movement of a guided particle $q$ is constrained to follow its guiding particle $p$ through,

$$ u_{(q)} = u_{(p)} + \varepsilon_{\text{imp}} (x_{(q)} - x_{(p)}) \qquad (27) $$

where $\varepsilon_{\text{imp}}$ is the imposed strain. This relation is enforced in the linear system via Lagrange multipliers.

**Remark 2.** The linear solver used in the previous versions of the beam-particle model was based on the Cholesky decomposition (Chen et al., 2008). When using Lagrange multipliers, the matrix of the linear system is no longer definite positive; thus, the Cholesky decomposition is no longer applicable. We use a linear solver based on the QR decomposition (Davis, 2011) to circumvent this issue.

## 2.3. Measurement of an effective elasticity tensor

To obtain the effective —i.e., damaged— elasticity tensor of a discrete specimen, the average strain and stress must be defined from the particle displacements and forces. Numerous definitions of the average strain tensor in discrete media are discussed by Bagi (2006). For the present study, the average strain tensor is defined as

$$ \bar{\varepsilon} = \frac{1}{S} \sum_{p=1}^{N} \left( \frac{\mathbf{u}_{(p)} + \mathbf{u}_{(p+1)}}{2} \odot \mathbf{n}_{(p,p+1)} \right) l_{(p,p+1)} \qquad (28) $$

where $S$ is the surface of the space cell system<fnref n="1" /> proposed by Bagi (1996) (see Figure 1), $\mathbf{u}_{(p)}$ is the displacement of the particle $p$, $\mathbf{n}_{(p,p+1)}$ is the outward-pointing normal to the beam linking particles $p$ and $p+1$, and $l_{(p,p+1)}$ is the length of the same beam, and where $\odot$ is the symmetrized tensorial product. For two vectors $\mathbf{a}$ and $\mathbf{b}$ it is such that $\mathbf{a} \odot \mathbf{b} = \frac{1}{2}(\mathbf{a} \otimes \mathbf{b} + \mathbf{b} \otimes \mathbf{a})$. The sum is carried over the $N$ particles on the boundary of the specimen.

The definition of the average stress tensor is a symmetrization of the definition proposed by (Bagi, 1996),

$$ \bar{\sigma} = \frac{1}{S} \sum_{p=1}^{N} f_{(p)} \odot \mathbf{x}_{(p)} \qquad (29) $$

<footnote n="1">¹The surface of the space cell system is the surface of the specimen that can be deformed; thus, the surface bounded by the lattice network. It is illustrated on Figure 1.</footnote>
where $f_{(p)}$ is the force applied to particle $p$ and $x_{(p)}$ is the position of particle $p$.

**Remark 3.** The assumption that the strain and stress are symmetric tensors is necessary for the present small strain framework. Due to the rotation of the particles, the beam-particle model shall be modeled by a generalized continuum. Indeed, some studies use higher-order continuum models to represent discrete media (Pradel and Sab, 1998; Ehlers et al., 2003; Dos Reis and Ganghoffer, 2012; Rezakhani and Cusatis, 2016).

Knowing three linearly independent strain tensors $\varepsilon^{(i)}$ and the associated stress tensors $\sigma^{(i)}$, in Kelvin notation,

$$ [\varepsilon^{(i)}] = \begin{bmatrix} \varepsilon_{xx}^{(i)} \\ \varepsilon_{yy}^{(i)} \\ \sqrt{2}\varepsilon_{xy}^{(i)} \end{bmatrix}, \quad [\sigma^{(i)}] = \begin{bmatrix} \sigma_{xx}^{(i)} \\ \sigma_{yy}^{(i)} \\ \sqrt{2}\sigma_{xy}^{(i)} \end{bmatrix} \qquad (30) $$

the elasticity tensor in Kelvin notation can be obtained as the symmetrized 3 × 3 matrix product

$$ [E] = \{ [\sigma^{(1)}], [\sigma^{(2)}], [\sigma^{(3)}] \} \cdot [\varepsilon^{(1)}], [\varepsilon^{(2)}], [\varepsilon^{(3)}]^{-1} \}^{\mathrm{S}} \qquad (31) $$

To obtain three linearly independent strain tensors, elastic periodic strain loadings are applied to the virtual specimen

$$ [\varepsilon_{\text{imp}}^{(1)}] = \begin{bmatrix} \varepsilon \\ 0 \\ 0 \end{bmatrix}, \quad [\varepsilon_{\text{imp}}^{(2)}] = \begin{bmatrix} 0 \\ \varepsilon \\ 0 \end{bmatrix}, \quad [\varepsilon_{\text{imp}}^{(3)}] = \begin{bmatrix} 0 \\ 0 \\ \sqrt{2}\varepsilon \end{bmatrix}, \qquad (32) $$

where $\varepsilon$ is sufficiently small (chosen such that the loading remains elastic). Algorithm 1 details the procedure to measure the evolution of the elasticity tensor during a mechanical loading.

**Algorithm 1:** Measurement of the evolution of effective elasticity tensor

Generate a virtual specimen (micro-structure);

Apply a damaging loading;

Create an empty list of effective elasticity tensor;

foreach load step do

Extract the cracks;

Add the cracks to uncracked virtual specimen;

foreach measurement loading i do

Apply measurement loading $[\varepsilon_{\text{imp}}^{i}]$;

Compute average strain $[\bar{\varepsilon}^{(i)}]$ Eq. (28);

Compute average stress $[\bar{\sigma}^{(i)}]$ Eq. (29);

end

Compute effective elasticity tensor Eq. (31);

Store effective elasticity tensor in the list;

end

## 2.4. Dataset of 76 356 effective elasticity tensors

This part is dedicated to generating a large dataset of effective elasticity tensors. To constitute the dataset, 36 virtual specimens with different micro-structures (but with the same macroscopic properties) are submitted to 21 mechanical loadings, uniaxial or multiaxial. Each specimen (also called Area Element) is a square of 0.2m×0.2m, with an average beam length $\bar{l}_b = 0.002$ m. Thus, the specimens are composed of 100 × 100 particles. Each mechanical loading is discretized into 100 loading steps. The database contains 36×21 = 756 evolutions of elasticity tensors, each containing 101 elasticity tensors, leading to a total of 76 356 elasticity tensors. Note that some elasticity tensors appear multiple times in the dataset (when the specimen is not yet damaged).

Different types of boundary conditions have been applied to generate this dataset:
<figure><img src="image_6.png" /><figcaption>**Figure 4:** Illustration of localized cracking patterns in the dataset.</figcaption></figure>

* Kinematic Uniform Boundary Conditions (KUBC), where strain is imposed on the whole boundary of the specimen,

* Periodic Boundary Conditions (PBC) as described in subsection 2.2, and

* Experimental Boundary Conditions (EXPE), where displacement is imposed on parts of the boundary.

Loadings are also separated into proportional and non-proportional loadings. The latter lead to a rotation of the principal strain/stress directions during the loading. It leads to a misalignment of the loading direction with first-stage micro-cracks, and a change of the damaging direction. The names, descriptions and parameters of the applied loadings are detailed in Appendix B. Examples of final micro-cracking patterns are provided in Figure 4. The measurement procedure is applied for an EXPE bi-tension loading in Figure 5 for illustration purposes.

**Remark 4.** Instead of using physical loadings to generate micro-cracks, one could try to create random cracking patterns by randomly breaking beams. This method does not account for the interactions of cracks during the loading.
<figure><img src="image_7.png" /><figcaption>**Figure 5:** Illustration of the method to measure the evolution of the effective elasticity tensor in the EXPE bi-tension loading.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 6:** Initial distance to isotropy for each micro-structure in the dataset.</figcaption></figure>

We have preferred to apply true mechanical loadings, as those constrain the micro-cracking patterns to "physically-reachable" ones.

## 2.5. Isotropy of the initial elasticity tensor $E_0$

For the present modeling, we suppose that the undamaged elasticity tensor is isotropic. To check that this is also the case for the undamaged (uncracked) elasticity tensors in the dataset, we compute the distance to isotropy, by Equation 22, for the 36 micro-structures in the dataset. The results are provided in Figure 6. The mean value of the initial relative distance to isotropy is $\Delta_{I_{\text{so}}}^{\text{mean}} = 0.0102$, whereas the standard deviation is $\Delta_{I_{\text{so}}}^{\text{mean}} = 0.0018$. These relative distances are all between $\Delta_{I_{\text{so}}}^{\text{min}} = 0.007$ and $\Delta_{I_{\text{so}}}^{\text{max}} = 0.016$. This is sufficient to consider that the elasticity tensors are initially isotropic.
<figure><img src="image_9.png" /><figcaption>Figure 7: Histograms of relative distances to symmetry classes over the dataset</figcaption></figure>

## 2.6. Distances to isotropy and to orthotropy of the effective elasticity tensors

The distances of the elasticity tensors to isotropy or to orthotropy can also be used to justify the tensorial nature of the damage variable (Oliver-Leblond et al., 2021). These distances have been computed for each elasticity tensor in the dataset. The corresponding histograms are plotted in Figure 7a and Figure 7b.

In Figure 7a, the distribution of the relative distance to isotropy (Equation 22) shows that a large part of effective elasticity tensors in the dataset is far from being isotropic. This means that a scalar (isotropic) damage variable is insufficient to represent the loss of stiffness (due to micro-cracking).

In the Figure 7b, the distribution of the relative distance to orthotropy (Equation 23) shows that most of the tensors in the dataset are close being orthotropic. This means that the effective elasticity tensor can be modeled as remaining orthotropic during the mechanical loadings of the whole dataset. From the reconstruction formula (15) and in accordance with (Desmorat and Desmorat, 2016), this implies that at most two second-order tensors are required to represent the impact of micro-cracking on the bi-dimensional elasticity tensor.

The question of whether the single second-order tensor **D** within the reconstruction formulas of section 1 is sufficient to represent this coupling has to be addressed. By the damage definition in Eq. (18), the bulk modulus $\kappa$ and the deviatoric part of the dilatation tensor are exactly modeled from the anisotropic damage variable **D**. To fully determine the effective elasticity tensor $\mathbf{E} = \mathbf{E}(\mathbf{D})$, the shear modulus $\mu$ and the harmonic part $\mathbf{H}$ of $\mathbf{E}$ need to be modeled as functions of the damage **D**. It proves essential to note that each term of the harmonic reconstruction formula (11) is orthogonal to each other (see (Blinowski et al., 1996; Vianello, 1997; Desmorat and Desmorat, 2015)). From a modeling point of view, this means that the modeling errors associated with each part **Iso**, **Dil** and **H** of the harmonic decomposition are independent. They can be handled separately. Conversely, this also means that a modeling error on the harmonic part cannot be compensated for by the isotropic part **Iso** or by the dilatation part **Dil**.

## 3. Modeling of the shear modulus-damage state coupling

We recall first the expression of the generalized shear modulus $\mu$ of the elasticity tensor $\mathbf{E}$,

$$ \mu = \frac{1}{8} (2 \operatorname{tr} \mathbf{v} - \operatorname{tr} \mathbf{d}), \quad \mathbf{d} = \mathbf{E}:\mathbf{1}, \quad \mathbf{v} = \operatorname{tr}_{13} \mathbf{E}. \qquad (33) $$

To express $\mu$ as a function of damage tensor **D**, a relation between the trace $\operatorname{tr} \mathbf{v}$ of the Voigt tensor and the damage tensor **D** has to be exhibited.
<figure><img src="image_10.png" /><figcaption>**Figure 8:** Damage $D_v$ as a function of the eigenvalues $D_1$ and $D_2$ ($D_2 > D_1$) of the damage variable **D** over the dataset. Each colored point is the value of $D_v$ for an elasticity tensor of the dataset. The black dots are the projections of the data points onto the planes ($D_1, D_v$), ($D_2, D_v$) and ($D_1, D_2$).</figcaption></figure>

To ease the modeling process, an intermediate scalar variable $D_v$ such that $\text{tr} \mathbf{v} = \text{tr} v_0(1 - D_v)$ is introduced,

$$D_v = \frac{\text{tr} v_0 - \text{tr} v}{\text{tr} v_0}. \qquad (34)$$

It can be interpreted as a damage variable based on $\text{tr}(v)$. However, since it represents the same micro-cracking pattern as **D**, it is not assumed to be an additional thermodynamics (internal) variable. It will be modeled as a function of the damage variable **D**.

Figure 8 contains a scatter plot of $D_v$ as a function of the two eigenvalues $D_1, D_2$ of the damage tensor **D** for each elasticity tensor of the dataset. It is observed that the points ($D_1, D_2, D_v$) of the dataset are grouped around a surface. Thus, $D_v$ can be modeled as a function of damage by a constitutive equation $D_v = D_v^m(\mathbf{D})$. In practice, we approximate $D_v$ by the linear combination of invariants of **D**,

$$D_v^m(\mathbf{D}) = c_1 I_1(\mathbf{D}) + c_2 I_2(\mathbf{D}), \quad I_k(\mathbf{D}) = \text{tr}(\mathbf{D}^k) = D_1^k + D_2^k, \qquad (35)$$

where the upperscript m stands for model, $I_k(\mathbf{D})$ are invariants of **D**, and $c_k$ are the parameters of the model. A second expression

$$D_v^m(\mathbf{D}) = c_1 I_1(\mathbf{D}) + c_2 I_2(\mathbf{D}) + c_3 I_3(\mathbf{D}), \quad I_3(\mathbf{D}) = \frac{1}{2} (3I_1(\mathbf{D}) I_2(\mathbf{D}) - I_1(\mathbf{D})^3), \qquad (36)$$

will also be studied (but the two-term expansion in Eq. (35) will prove sufficient).

The modeling process will be carried out in two steps. The first step consists in deriving and justifying some physical constraints on the material constants $c_k$. Those constraints will limit the number of independent parameters of the model. The second step consists in identifying the parameters $c_k$.

## 3.1. Physical constraints

*Constraint 1 – Undamaged state.* For an undamaged state, the trace of the Voigt tensor keeps its initial value, $\text{tr} v = \text{tr} v_0$,

$$D_v^m(\mathbf{D} = \mathbf{0}) = 0 \qquad (37)$$

This constraint is satisfied by both Eq. (35) and (36) since the invariants of **D** vanish ($I_k(\mathbf{D} = \mathbf{0}) = 0$).

*Constraint 2 – Fully damaged state.* For a fully damaged state, the effective elasticity tensor is a null fourth-order tensor. We have then $\text{tr} v = 0$, which implies $D_v^m(\mathbf{D} = \mathbf{1}) = 1$ and,

$$2 \sum_{k=1}^{n} c_k = 1, \quad n = 2, 3. \qquad (38)$$
<figure><img src="image_11.png" /><figcaption>**Figure 9:** Check of the assumption of total symmetry of the stiffness loss tensor: the dots in black correspond to the value in the dataset and the red dots correspond to the $D_v$ computed via Equation 41</figcaption></figure>

*Constraints 3 – Diffuse micro-cracking.* In the early degradation stage, the micro-cracking is diffuse (in the sense that it is not localized within the RAE). In Figure 8, this stage corresponds to the region where the points ($D_1, D_2, D_v$) are close to a plane. It corresponds to the assumption of non-interacting cracks (as defined by Kachanov (1992)).

Kachanov has shown that in 2D and as long as the micro-cracks are not interacting, the gain in compliance is a totally symmetric fourth-order tensor. This property is not satisfied anymore when the micro-cracking is localized. This observation guided us to check if this property is satisfied by the weakly damaged elasticity tensors of the dataset.

Assuming in the early damage stage that the stiffness loss $\Delta E = E - E_0$ is totally symmetric implies $\text{tr}_{12}(E - E_0) = \text{tr}_{13}(E - E_0)$, i.e., $\text{tr}\,d - \text{tr}\,d_0 = \text{tr}\,v - \text{tr}\,v_0$. This means that we have, at low damage,

$$D_v = \frac{\text{tr}\,d_0 - \text{tr}\,d}{\text{tr}\,v_0}. \qquad (39)$$

Taking the trace of the definition (20), giving $\text{tr}\,d = 2\kappa_0(2 - \text{tr}\,D)$, and using the relation between initial properties, $\text{tr}\,d_0 = 4\kappa_0$ and $\text{tr}\,v_0 = 4\mu_0 + 2\kappa_0$, we get, at low damage still,

$$D_v = \frac{\kappa_0}{2\mu_0 + \kappa_0} I_1(D), \qquad (40)$$

and

$$\left. \frac{\partial D_v}{\partial D} \right|_{D=0} = \frac{\kappa_0}{2\mu_0 + \kappa_0} \mathbf{1} \quad (\text{diffuse micro-cracking assumption}). \qquad (41)$$

Figure 9 provides a check of this assumption in the 76 356 elasticity tensors dataset. For small values of the damage (tr $D < 1$), the values of $D_v$ obtained by the diffuse micro-cracking assumption provide an accurate model for the values from the dataset.

For both constitutive equations (35) and (36), the diffuse micro-cracking constraint leads to

$$\left. \frac{\partial}{\partial D} \right|_{D=0} \left( \sum_{k=1}^{n} c_k I_k(D) \right) = \frac{\kappa_0}{2\mu_0 + \kappa_0} \mathbf{1}, \quad \text{i.e.,} \quad c_1 = \frac{\kappa_0}{2\mu_0 + \kappa_0}. \qquad (42)$$
## 3.2. Parameters $c_k$

For the first two-term expression in Eq. (35), the physical constraints leads to

$$c_1 = \frac{\kappa_0}{2\mu_0 + \kappa_0}, \quad \text{and} \quad c_2 = \frac{1}{2} - c_1, \qquad (43)$$

so that,

$$D_v^m(\mathbf{D}) = \frac{\kappa_0}{2\mu_0 + \kappa_0} (I_1(\mathbf{D}) - I_2(\mathbf{D})) + \frac{1}{2} I_2(\mathbf{D}). \qquad (44)$$

Note that this first modeling does not introduce any material parameter.

Applying the physical constraints to the three-term expression in Eq. (36) leads to

$$c_1 = \frac{\kappa_0}{2\mu_0 + \kappa_0}, \quad c_2 = \frac{1}{2} - c_1 - c_3, \qquad (45)$$

and

$$D_v^m(\mathbf{D}) = \frac{\kappa_0}{2\mu_0 + \kappa_0} (I_1(\mathbf{D}) - I_2(\mathbf{D})) + \frac{1}{2} I_2(\mathbf{D}) + c_3 (I_3(\mathbf{D}) - I_2(\mathbf{D})). \qquad (46)$$

The parameter $c_3$ is determined via a regression over the whole dataset. We get the small value

$$c_3 = 0.0973, \qquad (47)$$

to be compared to the particular case $c_3 = 0$ of the two-term expression in Eq. (44).

# 4. Modeling of the harmonic part-damage state coupling

The relative distance to orthotropy remains small for most elasticity tensors in our large dataset (as shown in Figure 7b). We thus make the simplifying assumption that the effective elasticity tensors E are orthotropic. By Eq. (14), the harmonic part of a 2D orthotropic elasticity tensor can be written as the harmonic square

$$\mathbf{H} = \frac{2K_3}{I_2^2} \mathbf{d}' * \mathbf{d}' = \pm \|\mathbf{H}\| \frac{\mathbf{d}' * \mathbf{d}'}{\|\mathbf{d}' * \mathbf{d}'\|}, \quad \|\mathbf{d}' * \mathbf{d}'\| = \frac{1}{\sqrt{2}} \mathbf{d}' : \mathbf{d}' = \frac{1}{\sqrt{2}} I_2, \qquad (48)$$

depending on the sign of the invariant of the elasticity tensor $K_3(E) = d:H:d$. Based on this parametrization, the modeling of the harmonic part of the elasticity tensors in the dataset can be carried out in two steps: (i) choosing an orientation (the sign) and (ii) modeling the harmonic part prefactor $||H||$ as a function of the tensorial damage variable **D**.

**Remark 5.** Note that orientation must be well-predicted when the relative norm of the harmonic part $||H||/||E||$ is large. When the harmonic part is small (when $||H||/||E|| \ll 1$), a misprediction of the orientation has a small effect on the predicted elasticity tensor.

## 4.1. Orientation of the harmonic part

A first indicator of the orientation of the harmonic part **H** is the sign of the invariant $K_3(E)$ or, in an equivalent manner, the sign of $K_3(E)^{1/3}$ (which is in GPa). Figure 10a shows the histogram of the invariant $K_3(E)^{1/3}$ over the dataset. It shows that the negative values of $K_3(E)^{1/3}$ are small in absolute, whereas the positive values of $K_3(E)^{1/3}$ are an order of magnitude higher. However, a large value of $K_3$ does not mean that the harmonic part is large, whereas a small value of $K_3$ does mean that the harmonic part is small (see Eq. (14) and remark 5).

Another orientation indicator is the norm of the difference between the normalized harmonic part and the normalized harmonic square

$$O_+(\mathbf{E}) = \left\| \frac{\mathbf{H}}{\|\mathbf{H}\|} - \frac{\mathbf{d}' * \mathbf{d}'}{\|\mathbf{d}' * \mathbf{d}'\|} \right\|^2. \qquad (49)$$

The histogram of the orientation indicator $O_+(\mathbf{E})$ over the dataset is plotted in Figure 10b. It mainly exhibits two peaks, one at 0, corresponding to the plus sign in Eq. (48), and one at 4, corresponding to the minus sign. In order to decide which sign is the best for our modeling, the data set is split into two colored parts:
<figure><img src="image_12.png" /><figcaption>Figure 10: Two indicators for the sign of the harmonic part.</figcaption></figure>

* a blue part when the harmonic part is negligible (||H||²/||E||² lower than 2%),

* and a red part when the harmonic part is significant (||H||²/||E||² larger than 2%).

Figure 10b shows that a higher number of elasticity tensors have a significant harmonic part in the positive orientation ($O_+(\mathbf{E}) = 0$) than in the negative orientation ($O_+(\mathbf{E}) = 4$). This observation is consistent with the analysis based on the invariant $K_3$. It also indicates that the preferred orientation is $+\frac{d'*d'}{||d'*d'||}$ (plus sign in Eq. (48)).

For the remaining of this section, we

* use the proportionality of the deviatoric part of the dilatation tensor with damage $d' = -2\kappa_0 \mathbf{D}'$ (by Eq. (20)), and

* model the harmonic part of the effective elasticity tensors with a plus sign,

by setting (with a slight abuse of notation)

$$ \mathbf{H} = \mathbf{H}(\mathbf{D}) = H^m(\mathbf{D}) \frac{\mathbf{D}'*\mathbf{D}'}{\mathbf{D}' : \mathbf{D}'}, \qquad (50) $$

where $H^m = \sqrt{2}||\mathbf{H}||$ is a positive function (to be determined) of the tensorial damage variable $\mathbf{D}$.

## 4.2. Modeling of the harmonic part

To identify the constitutive equation $H^m(\mathbf{D})$, the first step is to check if the norm $||\mathbf{H}||$ of the harmonic part can be represented by a function of the tensorial damage variable $\mathbf{D}$ (through its invariants). Figure 11 shows the harmonic part norm $||\mathbf{H}||$ versus the damage invariants $I_1(\mathbf{D}) = \text{tr} \mathbf{D}$ and $I_2(\mathbf{D}') = \mathbf{D}' : \mathbf{D}'$. Even if the discrepancy is large, especially in the region where $I_1(\mathbf{D}) > 1$, this figure indicates that it should be possible to approximate the norm $||\mathbf{H}||$ by a function of the two damage invariants $I_1(\mathbf{D}) = \text{and } I_2(\mathbf{D}')$.

**Remark 6.** Figure 11 shows that the norm of the harmonic part does not vanish when the damage variable $\mathbf{D}$ is equal to the second-order identity $\mathbf{1}$. Moreover, for the effective elasticity tensors in the dataset such that $I_1(\mathbf{D}) \approx 2$ and $I_2(\mathbf{D}') \approx 0$, different values of the norm $||\mathbf{H}||$ are associated with the same value of damage. This means that the definition (18) of the damage variable is insufficient in these few cases to fully represent the variations of the effective elasticity tensor due to highly interacting micro-cracks. To account for those variations, a second (internal) damage
<figure><img src="image_13.png" /><figcaption>Figure 11: Norm of the harmonic part (in GPa) as a function of the damage invariants $I_1(D)$ and $I_2(D') = D' : D'$.</figcaption></figure>

variable, possibly of higher order, could be used (see (Cormery and Welemane, 2010; Desmorat and Desmorat, 2016) for instance). Yet, the gain of accuracy might not be worth the increased modeling complexity.

As in the previous section, let us frame the state modeling of the function $H^m(D)$ by physical assumptions/constraints. The assumption of initial isotropy imposes that

$$H^m(D = 0) = 0. \qquad (51)$$

We also assume that the effective elasticity tensor $E$ vanishes when the material is fully damaged. This implies that the harmonic part $H$ must vanish when the damage grows to $D = 1$,

$$H^m(D = 1) = 0. \qquad (52)$$

Figure 11 shows that the harmonic part norm $||H||$ remains small when the damage invariant $I_2(D')$ is small. We can assume that

$$H^m(I_1(D), I_2(D') = 0) = 0. \qquad (53)$$

This new condition includes both previous assumptions.

**Remark 7.** The condition (53) is satisfied when the damage is purely hydrostatic (i.e., when $D = \frac{1}{2}(\text{tr } D) 1$). If the damage is hydrostatic, both the dilatation part $\text{Dil}$ and the harmonic part $H$ are obtained null (by the formula (11)). Thus, the modeled 2D effective elasticity tensor $E(D)$ is isotropic during a hydrostatic damaging.

A sparse regression method (the LASSO regression, Tibshirani (1996)) has been applied to a (multivariate) polynomial in the damage invariants $I_1(D)$ and $I_2(D')$. This optimization method aims at fitting a parametrized function with respect to data while penalizing the number of non-zero parameters. Note that, following Gaines et al. (2018), physical assumptions can be accounted for through additional minimization constraints. The constrained LASSO regression for the polynomial modeling of the function $H^m(D) = \sqrt{2}||H(D)||$ (introduced in Eq. (50)) recasts as the minimization problem

$$\min_{h_{n_1 n_2}} \sum_{i=1}^{N_{\text{ela}}} \left( \sqrt{2} ||\mathbf{H}_i|| - \sum_{n_1, n_2} h_{n_1 n_2} I_1(D_i)^{n_1} I_2(D_i')^{n_2} \right)^2 + \alpha \sum_{n_1, n_2} |h_{n_1 n_2}|, \qquad (54)$$

where $N_{\text{ela}} = 76$ 356 is the number of elasticity tensors in the dataset, $h_{n_1 n_2}$ are the coefficients of the sought polynomial, i.e., the parameters, and $\alpha$ is an arbitrary hyper-parameter of the method. The first term corresponds to a
<figure><img src="image_14.png" /><figcaption>**Figure 12:** Illustration of the compromise between the precision and the number of non-zero parameters in the constrained LASSO regression. The dashed line shows the best compromise for us ($\alpha$ = 0.0043 GPa).</figcaption></figure>

classical regression with a least-square error. The second term (with $\alpha$) penalizes the sum of the absolute values of the parameters $h_{n_1n_2}$. In practice, the higher $\alpha$, the fewer the non-zero parameters. The minimization has been carried out for different values of $\alpha$. The coefficient of determination of the regression

$$r^2 = 1 - \frac{\sum_{i=1}^{N_{ela}} (\sqrt{2}\|\mathbf{H}_i\| - \sum_{n_1, n_2} h_{n_1n_2} I_1(\mathbf{D}_i)^{n_1} I_2(\mathbf{D}_i')^{n_2})^2}{\sum_{i=1}^{N_{ela}} (\sqrt{2}\|\mathbf{H}_i\| - \text{mean}(\sum_{n_1, n_2} h_{n_1n_2} I_1(\mathbf{D}_i)^{n_1} I_2(\mathbf{D}_i')^{n_2}))^2} \quad (55)$$

is introduced as an indicator of the accuracy of the regression.

As shown in Figure 12, the best compromise between the model's accuracy (evaluated via the coefficient of determination $r^2$ and the number of non-zero parameters) has been obtained for the one-parameter polynomial expression

$$H^m(\mathbf{D}) = h I_1(\mathbf{D})^4 I_2(\mathbf{D}'), \quad h = h_{41} = 17 \text{ GPa}, \quad (56)$$

where $h$ is the so-called harmonic prefactor. Combined with equations (50) and (12), the previous expression means that the harmonic part of the elasticity tensors of our large dataset is well modeled by the simple constitutive equation

$$\mathbf{H}(\mathbf{D}) = h (\text{tr} \, \mathbf{D})^4 \mathbf{D}' * \mathbf{D}', \quad \mathbf{D}' * \mathbf{D}' = \mathbf{D}' \otimes \mathbf{D}' - \frac{1}{2} (\mathbf{D}' : \mathbf{D}') \mathbf{J}, \quad (57)$$

which expresses the harmonic part of the 2D effective elasticity tensors as a function of the damage variable $\mathbf{D}$ only. Recall that $\mathbf{J} = \mathbf{I} - \frac{1}{2}\mathbf{1} \otimes \mathbf{1}$ is the deviatoric projector.

### 4.3. Neglecting the harmonic part

As an alternative to the polynomial expression (57), it is also worth neglecting the harmonic part $\mathbf{H}$ and simply setting $h = 0$, i.e.,

$$\mathbf{H}(\mathbf{D}) = 0. \quad (58)$$

Indeed, as shown in Figure 13, the harmonic part is a small proportion of the effective elasticity tensor in most cases. Due to these observations, the harmonic part can often be neglected.

**Remark 8.** The modeling assumption $h = 0$, $\mathbf{H}(\mathbf{D}) = 0$, makes the effective elasticity tensor $r_0$-orthotropic in the sense of (Vannucci, 2002).
<figure><img src="image_15.png" /><figcaption>Figure 13: Proportion of each part of the elasticity tensor for each tensor of the dataset.</figcaption></figure>

## 5. Summary of the proposed anisotropic damage state coupling

Let us summarize the constitutive equations obtained in sections 1, 3, and 4. These equations are obtained from both the elasticity tensors reconstruction formulas (11) and (15), the performed virtual beam-particle computations, and cross-identification over the large effective elasticity tensors dataset. Let us also switch to standard Continuum Mechanics notations and mark with a tilde the effective (damaged) quantities, such as the effective elasticity tensor $\tilde{E} = \tilde{E}(D)$, the effective generalized shear and bulk moduli $\tilde{\mu} = \tilde{\mu}(D)$, $\tilde{\kappa} = \tilde{\kappa}(D)$, the effective harmonic tensor $\tilde{H} = \tilde{H}(D)$. We recall and rewrite the state equations obtained by using the relation

$$ \tilde{d} = 2\kappa_0 (1 - D), \qquad (59) $$

between the effective dilatation tensor $\tilde{d} = \tilde{E}:1$ and the second-order damage variable $D$.

The proposed anisotropic damage coupling is based on the damage dependency of each component of the harmonic decomposition

$$ \tilde{E}(D) = 2\tilde{\mu}(D)J + \tilde{\kappa}(D)1 \otimes 1 - \kappa_0 (1 \otimes D' + D' \otimes 1) + \tilde{H}(D), \qquad (60) $$

of the effective bi-dimensional elasticity tensor. The effective shear modulus has been derived as

$$ \tilde{\mu}(D) = \frac{1}{8} (2 \operatorname{tr} \tilde{v} - \operatorname{tr} \tilde{d}), \quad \operatorname{tr} \tilde{v} = \operatorname{tr}_{13} \operatorname{tr} \tilde{E} = \operatorname{tr} v_0 (1 - D_{\tilde{v}}^{\mathrm{m}(D)}) \qquad (61) $$

with for $D_{\tilde{v}}^{\mathrm{m}(D)}$, function of the invariants $\operatorname{tr} D, D:D$ and $I_3 = \operatorname{tr}(D^3)$, either the no-additional parameter function in Eq. (44) or the one-additional parameter function in Eq. (46). Defining the nonlinear shear-damage coupling parameter

$$ m = \frac{1}{2} (\kappa_0 + 2\mu_0) c_3, \qquad (62) $$

these two expressions are unified as the first constitutive equation of Table 2, in which are also recalled the expressions obtained for $\tilde{\kappa}(D)$ and $\tilde{H}(D)$.

The corresponding elasticity law coupled with anisotropic damage derives from a thermodynamics potential, the Helmholtz free energy density $\rho\psi$ function of the strain and the damage, as

$$ \sigma = \rho \frac{\partial \psi}{\partial \varepsilon} = \tilde{E}(D):\varepsilon, \qquad (63) $$

where

$$ \rho\psi = \frac{1}{2}\varepsilon:\tilde{E}(D):\varepsilon = \tilde{\mu}(D)\varepsilon':\varepsilon' + \frac{1}{2}\tilde{\kappa}(D)(\operatorname{tr}\varepsilon)^2 - \kappa_0 D':\varepsilon'\operatorname{tr}\varepsilon + \frac{1}{2}h(\operatorname{tr}D)^4\varepsilon':(D' * D'):\varepsilon', \qquad (64) $$
<table><caption>**Table 2**

Summary of the constitutive equations.</caption><thead><tr><th>Quantity</th><th>Model</th></tr></thead><tbody><tr><td>Shear modulus</td><td>μ̃(D) = μ<sub>0</sub> - 1/4 κ<sub>0</sub> tr D + 1/4 (κ<sub>0</sub> - 2μ<sub>0</sub>) D: D + m (D: D - tr(D<sup>3</sup>))</td></tr><tr><td>Bulk modulus</td><td>κ̃(D) = κ<sub>0</sub> (1 - 1/2 tr D)</td></tr><tr><td>Dilatation tensor</td><td>d̃ = 2κ<sub>0</sub> (1 - D)</td></tr><tr><td>Harmonic part</td><td>H̃(D) = h (tr D)<sup>4</sup> D' * D'</td></tr></tbody></table>

<table><caption>**Table 3**

Summary of the constitutive parameters (and their value for the 76 356 elasticity tensors dataset).</caption><thead><tr><th>Name</th><th>Symbol</th><th>Value in GPa</th></tr></thead><tbody><tr><td>Initial shear modulus</td><td>μ<sub>0</sub></td><td>19.4</td></tr><tr><td>Initial bulk modulus</td><td>κ<sub>0</sub></td><td>30</td></tr><tr><td>Nonlinear shear-damage coupling parameter</td><td>m</td><td>3.35</td></tr><tr><td>Harmonic prefactor</td><td>h</td><td>17</td></tr></tbody></table>

431 is a polynomial of simple and joint invariants of the state variables, i.e., of the strain and damage tensors. It details as

$$ \rho\psi(\varepsilon, \mathbf{D}) = \left( \mu_0 - \frac{1}{4}\kappa_0 \operatorname{tr} \mathbf{D} + \frac{1}{4}(\kappa_0 - 2\mu_0)\mathbf{D}:\mathbf{D} + m(\mathbf{D}:\mathbf{D} - \operatorname{tr}(\mathbf{D}^3)) \right) \varepsilon' : \varepsilon' \\ + \frac{1}{2}\kappa_0 \left(1 - \frac{1}{2}\operatorname{tr}\mathbf{D}\right) (\operatorname{tr}\varepsilon)^2 - \kappa_0 \mathbf{D}':\varepsilon' \operatorname{tr}\varepsilon + \frac{1}{2}h(\operatorname{tr}\mathbf{D})^4 \varepsilon' : (\mathbf{D}' * \mathbf{D}') : \varepsilon'. \quad (65) $$

432 Using definition (12), this state potential can be rewritten as

$$ \rho\psi(\varepsilon, \mathbf{D}) = \left( \mu_0 - \frac{1}{4}\kappa_0 \operatorname{tr} \mathbf{D} + \frac{1}{4}(\kappa_0 - 2\mu_0)\mathbf{D}:\mathbf{D} + m(\mathbf{D}:\mathbf{D} - \operatorname{tr}(\mathbf{D}^3)) - \frac{1}{4}h(\operatorname{tr}\mathbf{D})^4 \mathbf{D}':\mathbf{D}' \right) \varepsilon' : \varepsilon' \\ + \frac{1}{2}\kappa_0 \left(1 - \frac{1}{2}\operatorname{tr}\mathbf{D}\right) (\operatorname{tr}\varepsilon)^2 - \kappa_0 \mathbf{D}':\varepsilon' \operatorname{tr}\varepsilon + \frac{1}{2}h(\operatorname{tr}\mathbf{D})^4 (\mathbf{D}':\varepsilon')^2. \quad (66) $$

433 The thermodynamics force associated with the damage is then the symmetric second-order tensor

$$ \mathbf{Y} = -\rho \frac{\partial \psi}{\partial \mathbf{D}} \\ = \left( \frac{1}{4}\kappa_0 \mathbf{1} + \frac{1}{2}(2\mu_0 - \kappa_0)\mathbf{D} + m(3\mathbf{D}^2 - 2\mathbf{D}) + h(\operatorname{tr}\mathbf{D})^3 \left( \frac{1}{2}(\operatorname{tr}\mathbf{D})\mathbf{D}' + (\mathbf{D}':\mathbf{D}')\mathbf{1} \right) \right) \varepsilon' : \varepsilon' \\ + \frac{1}{4}\kappa_0(\operatorname{tr}\varepsilon)^2 \mathbf{1} + \kappa_0(\operatorname{tr}\varepsilon)\mathbf{D}' - h(\operatorname{tr}\mathbf{D})^4 (\mathbf{D}':\varepsilon')^2 \varepsilon'. \quad (67) $$

434 The parameters of the final model are the initial shear modulus $\mu_0$, the initial bulk modulus $\kappa_0$, the (optional) shear-damage coupling parameter $m = \frac{1}{2}(\kappa_0 + 2\mu_0)c_3$ and the harmonic prefactor $h$. They are summarized, with their value for our dataset, in Table 3.

435 It is worth pointing out that only one damage variable (the second-order tensor D) and only two material parameters (m and h) are introduced in this final state coupling for quasi-brittle materials between bi-dimensional elasticity and anisotropic damage.

## 6. Representativity of the proposed anisotropic damage state coupling

440 This section provides assessments over the 76 356 tensors dataset of several modeling choices—summarized in the previous section—for the coupling bi-dimensional elasticity-anisotropic damage:

441 • simplified modeling with vanishing parameters $m = 0$ and/or $h = 0$.
<table><caption>**Table 4**</caption><thead><tr><th>Percentile (error)</th><th>1 %</th><th>2 %</th><th>5 %</th><th>10 %</th></tr></thead><tbody><tr><td>m = 0</td><td>43.8 %</td><td>60.7 %</td><td>92.4 %</td><td>99.9 %</td></tr><tr><td>m = 3.35 GPa</td><td>44.0 %</td><td>71.6 %</td><td>96.8 %</td><td>99.9 %</td></tr></tbody></table>

Proportion of tensors with a relative error on the isotropic part below 1%, 2%, 5% and 10% for both value of *m* (results independent from the harmonic prefactor *h*).

* full modeling with non zero parameters *m* = 3.35 GPa and/or *h* = 17 GPa.

Each computed (micro-cracked) elasticity tensor **E**<sub>*i*</sub> from the dataset has an isotropic part **Iso**<sub>*i*</sub>, a dilatation part **Dil**<sub>*i*</sub> and a harmonic part **H**<sub>*i*</sub> (determined by harmonic decomposition (11)). For the assessments, the damage **D** = **D**<sub>*i*</sub> is taken as equal to the damage variable **D**<sub>*i*</sub> = **1** − **d**<sub>*i*</sub>/2κ<sub>0</sub> = **1** − **E**<sub>*i*</sub> : **1**/2κ<sub>0</sub> measured for **E**<sub>*i*</sub> in the dataset.

To be quantitative, the error between a tensor **E**<sub>*i*</sub> of the dataset and the anisotropic damage modeling **E**̃ = **E**(**D**<sub>*i*</sub>) is defined as

$$ \|\mathbf{E}_i - \tilde{\mathbf{E}}(\mathbf{D}_i)\| = \sqrt{\|\mathbf{Iso}_i - \tilde{\mathbf{Iso}}(\mathbf{D}_i)\|^2 + \|\mathbf{Dil}_i - \tilde{\mathbf{Dil}}(\mathbf{D}_i)\|^2 + \|\mathbf{H}_i - \tilde{\mathbf{H}}(\mathbf{D}_i)\|^2}, \quad (68) $$

where, by the formulas of previous section,

$$ \tilde{\mathbf{Iso}}(\mathbf{D}) = 2\tilde{\mu}(\mathbf{D}) \mathbf{J} + \tilde{\kappa}(\mathbf{D}) \mathbf{1} \otimes \mathbf{1}, \quad \tilde{\mathbf{Dil}}(\mathbf{D}) = -\kappa_0 (\mathbf{1} \otimes \mathbf{D}' + \mathbf{D}' \otimes \mathbf{1}), \quad \tilde{\mathbf{H}}(\mathbf{D}) = h (\mathrm{tr} \, \mathbf{D})^4 \mathbf{D}' * \mathbf{D}'. \quad (69) $$

Since the fourth-order isotropic (Iso), dilatation (Dil) and harmonic (H) parts are orthogonal, the errors associated with each part are uncorrelated. This implies that the parameter *m* influences the isotropic part only (by shear modulus expression of Table 2), and that the parameter *h* influences the harmonic part only (by harmonic part expression of Table 2). We point out that the dilatation part has no modeling error (thanks to the equality $\tilde{\mathbf{Dil}}(\mathbf{D}_i) = -\kappa_0 (\mathbf{1} \otimes \mathbf{D}'_i + \mathbf{D}'_i \otimes \mathbf{1}) = \mathbf{Dil}(\mathbf{E}_i)$).

## 6.1. Assessment of the isotropic part $\tilde{\mathbf{Iso}}(\mathbf{D})$

Let us first analyze the contribution of the isotropic part to the modeling error. A histogram of the relative error $||\mathbf{Iso}_i - \tilde{\mathbf{Iso}}(\mathbf{D}_i)||/||\mathbf{E}_0||$ on the isotropic part<fnref n="2" /> is plotted in Figure 14a. The relative error range [0, 1] is discretized in 500 intervals for this histogram, and all the following ones. It shows that both modeling, with vanishing value *m* = 0 and with non-zero value *m* = 3.35 GPa, are able to represent well the isotropic part of all micro-cracked tensors in the dataset (with a relative error below 10%). Figure 14b provides the Cumulative Distribution Function associated with this histogram. The numbers of tensors in each interval are cumulated with the relative error going from 0 to 1. Afterward, the cumulated number is divided by the total number of tensors to define the Cumulative Distribution Function (such that it reaches 1 when the error is 1). This plot shows that most isotropic parts are modeled with an error below 10%.

More quantitative results are provided in Table 4, which focuses on different percentiles of the histogram. Each percentile corresponds to a percentage of elasticity tensors **E**<sub>*i*</sub> with a modeling error —on the isotropic part— below the given threshold. For instance, the 5-percentile gives the proportion of tensors that are modeled with an error below 5%. This table is obtained from the Cumulative Distribution of the error by fixing the error (in the x-axis) and reading the associated cumulative proportion of tensors (in the y-axis). The Table 4 shows that the case *m* = 3.35 GPa models 71.6% of the effective elasticity tensor with less than 2% of error, whereas the case *m* = 0 models only 60.7% of them. Thus, the case *m* = 3.35 GPa provides a slightly more accurate model of the isotropic part than the simplified case *m* = 0. However, it requires the introduction of the additional material parameter *m* = 1/2 (κ<sub>0</sub> + 2μ<sub>0</sub>) c<sub>3</sub>. As the gain in accuracy is small compared to the increase of modeling/identification complexity by adding a parameter, we propose to retain *m* = 0 (i.e., Eq. (44)) for the modeling of the coupling of the shear modulus with anisotropic damage.

## 6.2. Assessment of the harmonic part

Let us now analyze the contribution of the harmonic part to the modeling error. Figure 15a shows the histogram of the relative error $||\mathbf{H}_i - \tilde{\mathbf{H}}(\mathbf{D}_i)||/||\mathbf{E}_0||$ on the harmonic part. It indicates that the case *h* = 17 GPa substantially

<footnote n="2">²Duplicated tensors during a loading (two successive tensors with no crack growth) are filtered out in this plot. The total number of remaining tensors is 60 232.</footnote>
<figure><img src="image_16.png" /><figcaption>**Figure 14:** Histograms and Cumulative Distribution Functions of the relative error on the isotropic part over the dataset for both values of *m* (results independent from the harmonic prefactor *h*).</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 15:** Histograms and Cumulative Distributions Functions of the relative error on the harmonic part over the dataset for both values of *h* (results independent from the parameter *m*).</figcaption></figure>

improves the modeling in comparison to neglecting the harmonic part (case *h* = 0). Note that for *h* = 17 GPa, the contribution of the error on the harmonic part is at the same level as the contribution due to the isotropic part. Figure 15b provides the Cumulative Distribution Function of this histogram. It shows that setting *h* = 17 GPa enables to model most harmonic parts of the dataset with an error below 5% whereas it is not the case for *h* = 0.

More quantitative results are provided in Table 5, where different percentile of the histogram are given (as in Table 4). Once again, a percentile corresponds to percentage of elasticity tensors with a modeling error – on the harmonic part – below a given threshold. It shows that setting *h* = 17 GPa enables the modeling of most harmonic
<table><caption>**Table 5**</caption><thead><tr><th>Percentile (error)</th><th>1%</th><th>2%</th><th>5%</th><th>10%</th></tr></thead><tbody><tr><td>h = 0</td><td>29.7%</td><td>44.7%</td><td>72.0%</td><td>94.7%</td></tr><tr><td>h = 17 GPa</td><td>39.3%</td><td>65.0%</td><td>98.4%</td><td>100%</td></tr></tbody></table>

Proportion of tensors with an error on harmonic part below 1%, 2%, 5% and 10% for both values of the harmonic prefactor $h$.

<figure><img src="image_18.png" /><figcaption>**Figure 16:** Histogram and Cumulative Distribution Functions of the total relative error on the effective elasticity tensors over the dataset for $m = 0$ and $h = 17$ GPa.</figcaption></figure>

<table><caption>**Table 6**</caption><thead><tr><th>Percentile (error)</th><th>1%</th><th>2%</th><th>5%</th><th>10%</th><th>20%</th><th>30%</th><th>40%</th></tr></thead><tbody><tr><td>(1 − d)E<sub>0</sub></td><td>5.63%</td><td>10.6%</td><td>21.0%</td><td>34.3%</td><td>49.0%</td><td>70.0%</td><td>97.3%</td></tr><tr><td>Iso(D)</td><td>5.87%</td><td>11.3%</td><td>21.3%</td><td>34.7%</td><td>50.3%</td><td>74.7%</td><td>99.7%</td></tr><tr><td>Iso(D) + Dil(D)</td><td>18.9%</td><td>27.8%</td><td>59.7%</td><td>94.1%</td><td>100%</td><td>100%</td><td>100%</td></tr><tr><td>E(D)</td><td>22.1%</td><td>40.1%</td><td>85.8%</td><td>99.8%</td><td>100%</td><td>100%</td><td>100%</td></tr></tbody></table>

Proportion of tensors with a total error below 1%, 2%, 5%, 10%, 20%, 30% and 40%.

parts of the dataset (98.4%) with an error below 5%, which is far higher than in the case for $h = 0$ (72.0 %). Thus, the case $h = 17$ GPa (using Eq. (56)) clearly improves the accuracy of the modeling.

### 6.3. Assessment of the complete elasticity tensor

Let us complete the previous independent results by analyzing the modeling error on the whole elasticity tensors. Below, we gradually add each ingredient of the modeling (the dilatation part, then the harmonic part) up to the complete modeling $\tilde{E}(D)$. This allows quantifying the exact contribution of the different parts of the anisotropic damage model.

Figure 16 shows the histogram and Cumulative Distribution Function of the relative error on the effective elasticity tensor for the different modelings. Different percentiles of the histogram are provided in Table 6. The first line of the table,

$$E_0(1-d), \quad d = \|E\|/\|E_0\|, \qquad (70)$$
<figure><img src="image_19.png" /><figcaption>**Figure 17:** Histograms and Cumulative Distribution Functions of total error for large values of damage. The histogram and cumulative distribution of total error over the elasticity tensors in the dataset is also reproduced.</figcaption></figure>

<table><caption>**Table 7**</caption><thead><tr><th>Percentile (error)</th><th>1%</th><th>2%</th><th>5%</th><th>10%</th><th>20%</th><th>30%</th><th>40%</th></tr></thead><tbody><tr><td>(1 &minus; d)E<sub>0</sub></td><td>1.66%</td><td>4.43%</td><td>12.6%</td><td>22.5%</td><td>31.7%</td><td>55.4%</td><td>95.9%</td></tr><tr><td>tilde E(D)</td><td>8.31%</td><td>27.1%</td><td>79.1%</td><td>99.7%</td><td>100%</td><td>100%</td><td>100%</td></tr></tbody></table>

Proportion of highly damaged tensors ($D_{max} > 0.7$) with a total error below 1%, 2%, 5%, 10%, 20%, 30% and 40%.

corresponds to an isotropic model (with scalar damage variable *d*) as it can be found in the literature (Mazars, 1984; Lemaitre and Chaboche, 1985; Lemaitre, 1992). The second line of the table $\tilde{\text{Iso}}(D)$ corresponds to the isotropic part only of the proposed modeling ($\tilde{\text{E}} = \tilde{\text{Iso}}(D)$ with $m = 0$). The third line $\tilde{\text{E}} = \tilde{\text{Iso}}(D) + \tilde{\text{Dil}}(D)$ (with $m = 0$ still, and $h = 0$) adds the modeling of the dilatation part, and the fourth line considers the whole damage model $\tilde{\text{E}} = \tilde{\text{Iso}}(D) + \tilde{\text{Dil}}(D) + \tilde{\text{H}}(D)$ with $m = 0$ and $h = 17$ GPa.

Both isotropic models provide a poor estimation of the elasticity tensors of the dataset: the effective elasticity tensors of a micro-cracked quasi-brittle material cannot be recovered satisfactorily with an isotropic damage model. Adding the dilatation part to the model, i.e., modeling anisotropic damage, clearly improves its accuracy as it enables to account for anisotropy induced by micro-cracking. Accounting for the harmonic part (by the second-order damage tensor *D* still) improves the modeling, bringing the proportion of the effective elasticity tensors with less than 5% of error from 59.7% (for the modeling with $h = 0$, $\tilde{\text{H}}(D) = 0$) to 85.8% (for $h = 17$ GPa).

Let us finally check whether it is possible or not to model the influence of micro-cracking on the elasticity tensor up to high levels of damage. This modeling feature is very complicated to gain from a micromechanics homogenization approach (Kachanov, 1993; Ponte Castañeda and Willis, 1995; Cormery and Welemane, 2010; Dormieux and Kondo, 2016).

In Figure 17, we reproduce the histogram and Cumulative Distribution of relative errors with a filter to keep only the elasticity tensors associated with high damage value (i.e., $D_{max} = \max(D_1, D_2) > 0.7$). The differences between the filtered and the unfiltered distributions correspond to the relative errors for tensors with low damage (diffuse micro-cracking).

The histogram in Figure 17 shows that an isotropic modeling does not describe well effective elasticity tensors at high damage values. Table 6 shows that only 12.5% of highly damaged tensors at $D_{max} > 0.7$ are reconstructed with
<figure><img src="image_20.png" /><figcaption>**Figure 18:** Comparison of the evolutions of elasticity tensors obtained via the two modelings ($m = h = 0$ in blue, and $m = 3.35$ GPa, $h = 17$ GPa in green) with the reference data ($\tilde{\mathbf{E}}$, in red) during a shear→tension loading.</figcaption></figure>

an error below 5% (versus only 21% for the whole dataset). Isotropic damage does not provide an accurate modeling in the strong micro-cracks interactions stage.

The complete anisotropic damage modeling $\tilde{\mathbf{E}} = \tilde{\mathbf{Iso}}(\mathbf{D}) + \tilde{\mathbf{Dil}}(\mathbf{D}) + \tilde{\mathbf{H}}(\mathbf{D})$ (with $m = 0$ and $h = 17$ GPa) is able to model 79.1% of the tensors with an error below 5% (and 99.7% of the tensors with an error below 10%). It thus provides an accurate damage state coupling for the effective elasticity tensor, even for highly damaged specimens. This result emphasizes the strength of the proposed approach based on a beam-particle model of the anisotropic damage of quasi-brittle materials.

## 6.4. Illustration of the modeling on a multiaxial non-proportional loading

In order to illustrate the ability of the modeling to represent the anisotropic damage state coupling, let us detail both the beam-particle model response and the anisotropic damage $\tilde{\mathbf{E}} = \tilde{\mathbf{E}}(\mathbf{D})$ response for the PBC shear→tension loading of Table 8 of Appendix B. This loading is one of the non-proportional loadings with the largest harmonic part in our effective elasticity tensors dataset. Two modelings are compared to the reference data: a first one with $m = 0$ and $h = 0$, and a second one with $m = 3.35$ GPa and $h = 17$ GPa.

The computed beam particle elasticity tensors $\mathbf{E}$ and the results of the modeling $\tilde{\mathbf{E}}(\mathbf{D})$ are reported in Figure 18. As the damage $\mathbf{D} = \mathbf{1} - \mathbf{E}:\mathbf{1}/2\kappa_0$ is known, the bulk modulus $\kappa$ and the deviatoric part of the dilatation tensor $\tilde{\mathbf{d}}'$ are exactly modeled. The shear modulus is correctly predicted by both $m = 0$ and $m = 3.35$ GPa cases, close to each other. This loading leads to the significant growth of the harmonic part, which is neglected when $h = 0$. Setting $h = 17$ GPa provides a good modeling of the harmonic part.

Figure 19 displays the relative errors associated with the isotropic and harmonic parts during this loading. It confirms that the dilatation part is exactly modeled. For $h = 0$, the most significant part of the error is due to the modeling of the harmonic part. Note that for $h = 17$ GPa, the modeling of both the isotropic and the harmonic parts contribute at the same level of the total error on the effective elasticity tensor.
<figure><img src="image_21.png" /><figcaption>Figure 19: Evolution of relative errors during an EXPE tensile loading. Plain curves correspond to $m = 3.35$ GPa and $h = 17$ GPa, and dashed curves to $m = 0$ and $h = 0$.</figcaption></figure>

## 7. Conclusion

We have performed 2D beam-particle simulations of Area Elements submitted to various complex loadings and have generated a dataset of 76 356 effective (damaged) elasticity tensors. High levels of anisotropic damage have been reached. Those correspond to discrete computations with strong micro-cracks interactions and multiple coalescences.

By analyzing the distance to orthotropy over the whole dataset, we have first shown that (at least) one second-order damage variable is necessary to represent the effect of micro-cracking on the effective elastic tensors in our dataset. Thanks to a reconstruction formula of the orthotropic elasticity tensor by means of its covariants, we have proposed an anisotropic damage state coupling by a single second-order damage tensor. By the proper definition (18) of the second-order damage variable $\mathbf{D}$, the isotropic part associated with the effective bulk modulus $\tilde{\kappa}$ and the effective dilatation part $\tilde{\mathbf{Dil}} = \frac{1}{2}(\mathbf{1} \otimes \tilde{\mathbf{d}}' + \tilde{\mathbf{d}}' \otimes \mathbf{1})$ are exactly reconstructed from the damage $\mathbf{D}$ (Table 2). The orthogonality property of the isotropic/dilatation/harmonic parts of the harmonic decomposition has allowed us to build independent constitutive equations: a scalar one

$$ \tilde{\mu}(\mathbf{D}) = \mu_0 - \frac{1}{4}\kappa_0 \mathrm{tr} \mathbf{D} + \frac{1}{4}(\kappa_0 - 2\mu_0) \mathbf{D}:\mathbf{D} + m (\mathbf{D}:\mathbf{D} - \mathrm{tr}(\mathbf{D}^3)), $$

for the effective shear modulus $\tilde{\mu}$, and a tensorial one, polynomial in the damage variable

$$ \tilde{\mathbf{H}}(\mathbf{D}) = h (\mathrm{tr} \mathbf{D})^4 \mathbf{D}' * \mathbf{D}', = h (\mathrm{tr} \mathbf{D})^4 (\mathbf{D}' \otimes \mathbf{D}' - \frac{1}{2} (\mathbf{D}' : \mathbf{D}') \mathbf{J}), \quad (71) $$

for the fourth-order harmonic part of the effective elasticity tensor. In addition to the initial elasticity constants $\mu_0$, $\kappa_0$, only two material parameters have been introduced: the nonlinear shear-damage coupling parameter $m$ (which can furthermore be taken equal to zero) and the harmonic prefactor $h$ (identified as $h = 17$ GPa). The proposed anisotropic damage state coupling models 85.8% of the effective elasticity tensors in the dataset with less than 5% of error, including those with strong micro-cracks interactions and multiple coalescences.

## Acknowledgements

This work was performed using HPC resources from the “Mésocentre” computing center of CentraleSupélec, École Normale Supérieure Paris-Saclay and Université Paris-Saclay supported by CNRS and Région Île-de-France (https://mesocentre.universite-paris-saclay.fr/).
Abramian, S., Desmorat, B., Desmorat, R., Kolev, B., and Olive, M. (2020). Recovering the Normal Form and Symmetry Class of an Elasticity Tensor. *Journal of Elasticity*, 142(1):1–33.

Antonelli, A., Desmorat, B., Kolev, B., and Desmorat, R. (2022). Distance to plane elasticity orthotropy by Euler–Lagrange method. *Comptes Rendus. Mécanique*, 350(G2):413–430.

Auffray, N., Kolev, B., and Petitot, M. (2014). On Anisotropic Polynomial Relations for the Elasticity Tensor. *Journal of Elasticity*, 115(1):77–103.

Bagi, K. (1996). Stress and strain in granular assemblies. *Mechanics of Materials*, 22(3):165–177.

Bagi, K. (2006). Analysis of microstructural strain tensors for granular assemblies. *International Journal of Solids and Structures*, 43(10):3166–3184.

Bažant, Z. P. and Gambarova, P. G. (1984). Crack Shear in Concrete: Crack Band Microplane Model. *Journal of Structural Engineering*, 110(9):2015–2035. Publisher: American Society of Civil Engineers.

Bažant, Z. P. and Oh, B. H. (1985). Microplane Model for Progressive Fracture of Concrete and Rock. *Journal of Engineering Mechanics*, 111(4):559–582. Publisher: American Society of Civil Engineers.

Bažant, Z. P. and Prat, P. C. (1988a). Microplane Model for Brittle-Plastic Material: I. Theory. *Journal of Engineering Mechanics*, 114(10):1672–1688. Publisher: American Society of Civil Engineers.

Bažant, Z. P. and Prat, P. C. (1988b). Microplane Model for Brittle-Plastic Material: II. Verification. *Journal of Engineering Mechanics*, 114(10):1689–1702. Publisher: American Society of Civil Engineers.

Bažant, Z. P., Tabbara, M. R., Kazemi, M. T., and Pijaudier-Cabot, G. (1990). Random Particle Model for Fracture of Aggregate or Fiber Composites. *Journal of Engineering Mechanics*, 116(8):1686–1705. Publisher: American Society of Civil Engineers.

Blinowski, A., Ostrowska-Maciejewska, J., and Rychlewski, J. (1996). Two-dimensional Hooke's tensors - isotropic decomposition, effective symmetry criteria. *Archives of Mechanics*, 48(2):325–345. Number: 2.

Bolander, J. E., Eliáš, J., Cusatis, G., and Nagai, K. (2021). Discrete mechanical models of concrete fracture. *Engineering Fracture Mechanics*, 257:108030.

Bolander, J. E., Shiraishi, T., and Isogawa, Y. (1996). An adaptive procedure for fracture simulation in extensive lattice networks. *Engineering Fracture Mechanics*, 54(3):325–334.

Chaboche, J.-L. (1978). Description thermodynamique et phénoménologique de la viscoélasticité cyclique avec endommagement. *Description thermodynamique et phénoménologique de la viscoélasticité cyclique avec endommagement*.

Chaboche, J.-L. (1979). Le concept de contrainte effective appliqué à l'élasticité et à la viscoplasticité en présence d'un endommagement anisotrope. *Col. Euromech*, 115:737–760.

Chaboche, J.-L. (1984). Anisotropic creep damage in the framework of continuum damage mechanics. *Nuclear Engineering and Design*, 79(3):309–319.

Challamel, N., Picandet, V., and Pijaudier-Cabot, G. (2015). From discrete to nonlocal continuum damage mechanics: Analysis of a lattice system in bending using a continualized approach. *International Journal of Damage Mechanics*, 24(7):983–1012. Publisher: SAGE Publications Ltd STM.

Chen, Y., Davis, T. A., Hager, W. W., and Rajamanickam, S. (2008). Algorithm 887: CHOLMOD, Supernodal Sparse Cholesky Factorization and Update/Downdate. *ACM Transactions on Mathematical Software*, 35(3):22:1–22:14.

Cordebois, J.-P. and Sidoroff, F. (1980). Anisotropic Damage in Elasticity and Plasticity. *Journal de mecanique theorique et appliquee*, pages 45–59.

Cordebois, J.-P. and Sidoroff, F. (1982). Damage Induced Elastic Anisotropy. pages 761–774, Dordrecht. Springer Netherlands. Book Title: Mechanical Behavior of Anisotropic Solids / Comportement Mécanique des Solides Anisotropes.

Cormery, F. and Welemane, H. (2010). A stress-based macroscopic approach for microcracks unilateral effect. *Computational Materials Science*, 47(3):727–738.

Cundall, P. A. and Strack, O. D. L. (1979). A discrete numerical model for granular assemblies. *Géotechnique*, 29(1):47–65. Publisher: ICE Publishing.

D'Addetta, G. A., Kun, F., and Ramm, E. (2002). On the application of a discrete model to the fracture process of cohesive granular materials. *Granular Matter*, 4(2):77–90.

Davis, T. A. (2011). Algorithm 915, SuiteSparseQR: Multifrontal multithreaded rank-revealing sparse QR factorization. *ACM Transactions on Mathematical Software*, 38(1):8:1–8:22.

de Arcangelis, L. and Herrmann, H. J. (1989). Scaling and multiscaling laws in random fuse networks. *Physical Review B*, 39(4):2678–2684. Publisher: American Physical Society.

Delaplace, A. (2008). *Modélisation discrète appliquée au comportement des matériaux et des structures*. Mémoire d'habilitation à diriger des recherches, Ecole Normale Supérieure de Cachan.

Delaplace, A. and Desmorat, R. (2008). Discrete 3D model as complimentary numerical testing for anisotropic damage. *International Journal of Fracture*, 148(2):115–128.

Delaplace, A., Pijaudier-Cabot, G., and Roux, S. (1996). Progressive damage in discrete models and consequences on continuum modelling. *Journal of the Mechanics and Physics of Solids*, 44(1):99–136.

Desmorat, B. and Desmorat, R. (2015). Tensorial Polar Decomposition of 2D fourth-order tensors. *Comptes Rendus Mécanique*, 343(9):471–475.

Desmorat, B. and Desmorat, R. (2016). Second order tensorial framework for 2D medium with open and closed cracks. *European Journal of Mechanics - A/Solids*, 58:262–277.

Desmorat, R., Desmorat, B., Olive, M., and Kolev, B. (2018). Micromechanics based framework with second-order damage tensors. *European Journal of Mechanics - A/Solids*, 69:88–98.

Diner, Ç., Kochetov, M., and Slawinski, M. A. (2011). Identifying Symmetry Classes of Elasticity Tensors Using Monoclinic Distance Function. *Journal of Elasticity*, 102(2):175–190.
624 Dormieux, L. and Kondo, D. (2016). *Micromechanics of Fracture and Damage*. Wiley, 1 edition.

625 Dos Reis, F. and Ganghoffer, J. F. (2012). Construction of micropolar continua from the asymptotic homogenization of beam lattices. *Computers & Structures*, 112-113:354–363.

626 Ehlers, W., Ramm, E., Diebels, S., and D’Addetta, G. A. (2003). From particle ensembles to Cosserat continua: homogenization of contact forces towards stresses and couple stresses. *International Journal of Solids and Structures*, 40(24):6681–6702.

627 Fassin, M., Eggersmann, R., Wulfinghoff, S., and Reese, S. (2019). Gradient-extended anisotropic brittle damage modeling using a second order damage tensor – Theory, implementation and numerical examples. *International Journal of Solids and Structures*, 167:93–126.

628 Fichant, S., La Borderie, C., and Pijaudier-Cabot, G. (1999). Isotropic and anisotropic descriptions of damage in concrete structures. *Mechanics of Cohesive-frictional Materials*, 4(4):339–359. Publisher: Wiley.

629 François, M. (1995). *Identification des symétries matérielles de matériaux anisotropes*. phdthesis, Université Pierre et Marie Curie - Paris VI.

630 François, M., Geymonat, G., and Berthaud, Y. (1998). Determination of the symmetries of an experimentally determined stiffness tensor: Application to acoustic measurements. *International Journal of Solids and Structures*, 35(31-32):4091–4106.

631 Gaines, B. R., Kim, J., and Zhou, H. (2018). Algorithms for Fitting the Constrained Lasso. *Journal of Computational and Graphical Statistics*, 27(4):861–871. Publisher: Taylor & Francis _eprint: https://doi.org/10.1080/10618600.2018.1473777.

632 Gazis, D. C., Tadjbakhsh, I., and Toupin, R. A. (1963). The elastic tensor of given symmetry nearest to an anisotropic elastic tensor. *Acta Crystallographica*, 16(9):917–922.

633 Grassl, P. and Jirásek, M. (2006). Damage-plastic model for concrete failure. *International Journal of Solids and Structures*, 43(22):7166–7196.

634 Halm, D. and Dragon, A. (1996). A Model of Anisotropic Damage by Mesocrack Growth; Unilateral Effect. *International Journal of Damage Mechanics*, 5(4):384–402. Publisher: SAGE Publications Ltd STM.

635 Herrmann, H. J., Hansen, A., and Roux, S. (1989). Fracture of disordered, elastic lattices in two dimensions. *Physical Review B*, 39(1):637–648. Publisher: American Physical Society.

636 Herrmann, H. J. and Roux, S., editors (1990). *Statistical Models for the Fracture of Disordered Media*. Random Materials and Processes. North-Holland, Amsterdam.

637 Hrennikoff, A. (1941). Solution of Problems of Elasticity by the Framework Method. *Journal of Applied Mechanics*, 8(4):A169–A175.

638 Kachanov, L. M. (1958). On Creep Rupture Time. *Izv. Acad. Nauk SSSR, Otd. Techn. Nauk*, 8:26–31.

639 Kachanov, M. (1992). Effective Elastic Properties of Cracked Solids: Critical Review of Some Basic Concepts. *Applied Mechanics Reviews*, 45(8):304–335. Publisher: American Society of Mechanical Engineers Digital Collection.

640 Kachanov, M. (1993). Elastic Solids with Many Cracks and Related Problems. In Hutchinson, J. W. and Wu, T. Y., editors, *Advances in Applied Mechanics*, volume 30, pages 259–445. Elsevier.

641 Kawai, T. (1978). New discrete models and their application to seismic response analysis of structures. *Nuclear Engineering and Design*, 48(1):207–229.

642 Krajcinovic, D. (1985). Continuous Damage Mechanics Revisited: Basic Concepts and Definitions. *Journal of Applied Mechanics*, 52(4):829–834.

643 Krajcinovic, D. (1996). *Damage Mechanics*. Elsevier.

644 Kun, F. and Herrmann, H. J. (1996). A study of fragmentation processes using a discrete element method. *Computer Methods in Applied Mechanics and Engineering*, 138(1):3–18.

645 Ladevèze, P. (1983). Sur une théorie de l'endommagement anisotrope. Rapport Interne 34, LMT Cachan.

646 Ladevèze, P. (1995). Modelling and simulation of the mechanical behaviour of CMCs. *Ceram. Trans. Vol. 57*, pages 53–63.

647 Landis, E. N. (1999). Micro–macro fracture relationships and acoustic emissions in concrete. *Construction and Building Materials*, 13(1):65–72.

648 Leckie, F. A. and Onat, E. T. (1981). Tensorial Nature of Damage Measuring Internal Variables. In Hult, J. and Lemaitre, J., editors, *Physical Non-Linearities in Structural Analysis*, International Union of Theoretical and Applied Mechanics, pages 140–155, Berlin, Heidelberg. Springer.

649 Lemaitre, J. (1984). How to use damage mechanics. *Nuclear Engineering and Design*, 80(2):233–245.

650 Lemaitre, J. (1992). *A Course on Damage Mechanics*. Springer-Verlag, Berlin Heidelberg.

651 Lemaitre, J. and Chaboche, J.-L. (1985). *Mécanique des matériaux solides*. Dunod, english translation 1990 'Mechanics of Solid Materials' Cambridge University Press.

652 Lemaitre, J. and Desmorat, R. (2005). *Engineering Damage Mechanics: Ductile, Creep, Fatigue and Brittle Failures*. Springer-Verlag, Berlin Heidelberg.

653 Lemaitre, J., Desmorat, R., and Sauzay, M. (2000). Anisotropic damage law of evolution. *European Journal of Mechanics - A/Solids*, 19(2):187–208.

654 Lubarda, V. A. and Krajcinovic, D. (1993). Damage tensors and the crack density distribution. *International Journal of Solids and Structures*, 30(20):2859–2877.

655 Mazars, J. (1984). *Application de la mécanique de l'endommagement au comportement non-linéaire et à la rupture du béton de structure*. THESE DE DOCTEUR ES SCIENCES, Université Pierre et Marie Curie - Paris VI - Laboratoire de Mécanique et Technologie.

656 Mazars, J., Berthaud, Y., and Ramtani, S. (1990). The unilateral behaviour of damaged concrete. *Engineering Fracture Mechanics*, 35(4):629–635.

657 Meguro, K. and Hakuno, M. (1989). Fracture Analyses of Concrete Structures by the Modified Distinct Element Method. *Doboku Gakkai Ronbunshu*, 1989(410):113–124.

658 Meurer, A., Smith, C. P., Paprocki, M., Čertík, O., Kirpichev, S. B., Rocklin, M., Kumar, A., Ivanov, S., Moore, J. K., Singh, S., Rathnayake, T., Vig, S., Granger, B. E., Muller, R. P., Bonazzi, F., Gupta, H., Vats, S., Johansson, F., Pedregosa, F., Curry, M. J., Terrel, A. R., Roučka, Š., Saboo, A., Fernando, I., Kulal, S., Cimrman, R., and Scopatz, A. (2017). SymPy: symbolic computing in Python. *PeerJ Computer Science*, 3:e103. Publisher: PeerJ Inc.

659 Migliori, A., Sarrao, J. L., Visscher, W. M., Bell, T. M., Lei, M., Fisk, Z., and Leisure, R. G. (1993). Resonant ultrasound spectroscopic techniques for measurement of the elastic moduli of solids. *Physica B: Condensed Matter*, 183(1):1–24.

660 Moakher, M. and Norris, A. N. (2006). The Closest Elastic Tensor of Arbitrary Symmetry to an Elasticity Tensor of Lower Symmetry. *Journal of Elasticity*, 85(3):215–263.

661 Murakami, S. (1988). Mechanical Modeling of Material Damage. *Journal of Applied Mechanics*, 55(2):280–286.
687 Murakami, S. (2012). *Continuum Damage Mechanics*, volume 185 of *Solid Mechanics and Its Applications*. Springer Netherlands, Dordrecht.

688 Murakami, S. and Ohno, N. (1978). A constitutive equation of creep damage in pollicristalline metals. In *IUTAM Colloquium Euromech*, volume 111.

689 Olive, M., Kolev, B., Desmorat, B., and Desmorat, R. (2017). Harmonic Factorization and Reconstruction of the Elasticity Tensor. *Journal of Elasticity*, 132(1):67–101.

690 Olive, M., Kolev, B., Desmorat, R., and Desmorat, B. (2022). Characterization of the symmetry class of an elasticity tensor using polynomial covariants. *Mathematics and Mechanics of Solids*, 27(1):144–190. Publisher: SAGE Publications Ltd STM.

691 Oliver-Leblond, C. (2019). Discontinuous crack growth and toughening mechanisms in concrete: A numerical study based on the beam-particle approach. *Engineering Fracture Mechanics*, 207:1–22.

692 Oliver-Leblond, C., Desmorat, R., and Kolev, B. (2021). Continuous anisotropic damage as a twin modelling of discrete bi-dimensional fracture. *European Journal of Mechanics - A/Solids*, 89:104285.

693 Papa, E. and Taliercio, A. (1996). Anisotropic damage model for the multiaxial static and fatigue behaviour of plain concrete. *Engineering Fracture Mechanics*, 55(2):163–179.

694 Ponte Castañeda, P. and Willis, J. R. (1995). The effect of spatial distribution on the effective behavior of composite materials and cracked media. *Journal of the Mechanics and Physics of Solids*, 43(12):1919–1951.

695 Pradel, F. and Sab, K. (1998). Cosserat modelling of elastic periodic lattice structures. *Comptes Rendus de l'Académie des Sciences - Series IIB - Mechanics-Physics-Astronomy*, 326(11):699–704. Publisher: Elsevier Masson.

696 Rabotnov, Y. N. (1969). Creep Problems in Structural Members. *North-Holland Publishing Company, Amsterdam*.

697 Ramtani, S., Berthaud, Y., and Mazars, J. (1992). Orthotropic behavior of concrete with directional aspects: modelling and experiments. *Nuclear Engineering and Design*, 133(1):97–111.

698 Rezakhani, R. and Cusatis, G. (2016). Asymptotic expansion homogenization of discrete fine-scale models with rotational degrees of freedom for the simulation of quasi-brittle materials. *Journal of the Mechanics and Physics of Solids*, 88:320–345.

699 Richard, B. and Ragueneau, F. (2013). Continuum damage mechanics based model for quasi brittle materials subjected to cyclic loadings: Formulation, numerical implementation and applications. *Engineering Fracture Mechanics*, 98:383–406.

700 Rinaldi, A. and Lai, Y.-C. (2007). Statistical damage theory of 2D lattices: Energetics and physical foundations of damage parameter. *International Journal of Plasticity*, 23(10):1796–1825.

701 Rossi, P. and Richer, S. (1987). Numerical modelling of concrete cracking based on a stochastic approach. *Materials and Structures*, 20(5):334–337.

702 Roux, J., Hosten, B., Castagnede, B., and Deschamps, M. (1985). Caractérisation mécanique des solides par spectro-interférométrie ultrasonore. *Revue de Physique Appliquée*, 20(6):351–358. Publisher: Société Française de Physique.

703 Schlangen, E. and van Mier, J. G. M. (1992). Simple lattice model for numerical simulation of fracture of concrete materials and structures. *Materials and Structures*, 25(9):534–542.

704 Tibshirani, R. (1996). Regression Shrinkage and Selection via the Lasso. *Journal of the Royal Statistical Society. Series B (Methodological)*, 58(1):267–288. Publisher: [Royal Statistical Society, Wiley].

705 Vakulenko, A. A. and Kachanov, M. (1971). Continuum theory of medium with cracks. *Mekhanika tverdogo tela*, 4:159–166.

706 van Mier, J. G. M., van Vliet, M. R. A., and Wang, T. K. (2002). Fracture mechanisms in particle composites: statistical aspects in lattice type analysis. *Mechanics of Materials*, 34(11):705–724.

707 Vannucci, P. (2002). A Special Planar Orthotropic Material. *Journal of elasticity and the physical science of solids*, 67(2):81–96.

708 Vannucci, P. (2005). Plane Anisotropy by the Polar Method*. *Meccanica*, 40(4):437–454.

709 Vassaux, M. (2015). *Comportement mécanique des matériaux quasi-fragiles sous sollicitations cycliques : de l'expérimentation numérique au calcul de structures*. phdthesis, École normale supérieure de Cachan - ENS Cachan.

710 Vassaux, M., Oliver-Leblond, C., Richard, B., and Ragueneau, F. (2016). Beam-particle approach to model cracking and energy dissipation in concrete: Identification strategy and validation. *Cement and Concrete Composites*, 70:1–14.

711 Verchery, G. (1982). Les Invariants des Tenseurs d'Ordre 4 du Type de l'élasticité. In Boehler, J.-P., editor, *Mechanical Behavior of Anisotropic Solids / Comportement Mécanique des Solides Anisotropes*, pages 93–104, Dordrecht. Springer Netherlands.

712 Vianello, M. (1997). An integrity basis for plane elasticity tensors. *Archives of Mechanics*, 49(1):197–208. Number: 1.

713 William, K., Pramono, E., and Sture, S. (1989). Fundamental Issues of Smeared Crack Models. In Shah, S. P. and Swartz, S. E., editors, *Fracture of Concrete and Rock*, pages 142–157, New York, NY. Springer.

## A. Distance to bi-dimensional elastic orthotropy

This appendix provides a summary of the calculation of the exact distance of a 2D elasticity tensor to orthotropy.

We refer to (Antonelli et al., 2022) for the details.

The first step is to parametrize the orthotropic tensors $E^*$ from Eq. (23), in Kelvin notation,

$$ E^* = Q \star A, \quad Q = \begin{pmatrix} \cos \theta & -\sin \theta \\ \sin \theta & \cos \theta \end{pmatrix}, \quad [A] = \begin{bmatrix} A_{1111} & A_{1122} & 0 \\ A_{1122} & A_{2222} & 0 \\ 0 & 0 & 2A_{1212} \end{bmatrix}, \qquad (72) $$
where $Q$ is a rotation matrix, $A$ is the normal form of orthotropic 2D elasticity tensors, and $[A]$ stands for its Kelvin notation. With this parametrization, the minimization problem becomes

$$d_{\mathcal{O}_{rt}}(\mathbf{E}) = \min_{\theta, A} \|\mathbf{E} - Q \star A\|, \quad \Delta_{\mathcal{O}_{rt}}(\mathbf{E}) = \frac{d_{\mathcal{O}_{rt}}(\mathbf{E})}{\|\mathbf{E}\|}, \qquad (73)$$

where $\|\mathbf{E}\| = \sqrt{E_{ijkl}E_{ijkl}}$. By applying the harmonic decomposition to $\mathbf{E}$ and $\mathbf{A}$, and using the orthogonality of its isotropic/dilatation/harmonic parts, we have

$$d_{\mathcal{O}_{rt}}(\mathbf{E})^2 = \min_{\theta, A} (\|\mathbf{Iso}(\mathbf{E}) - \mathbf{Iso}(\mathbf{A})\|^2 + \|\mathbf{Dil}(\mathbf{E}) - Q \star \mathbf{Dil}(\mathbf{A})\|^2 + \|\mathbf{H}(\mathbf{E}) - Q \star \mathbf{H}(\mathbf{A})\|^2). \qquad (74)$$

For the isotropic part, the closest isotropic tensor is the isotropic tensor itself, $\mathbf{Iso}(\mathbf{A}) = \mathbf{Iso}(\mathbf{E})$ (by Eq. (22)). Furthermore:

* For the dilatation part term, we can introduce the deviatoric part of the dilatation tensor

$$\mathbf{d}'(\mathbf{E}) = \begin{bmatrix} d'_{11} & d'_{12} \\ d'_{11} & -d'_{11} \end{bmatrix}, \quad \mathbf{d}'(\mathbf{A}^*) = d^* \begin{bmatrix} 1 & 0 \\ 0 & -1 \end{bmatrix}, \quad Q \star \mathbf{d}'(\mathbf{A}) = d^* \begin{bmatrix} \cos 2\theta & \sin 2\theta \\ \sin 2\theta & -\cos 2\theta \end{bmatrix}, \qquad (75)$$

where

$$d'_{11} = \frac{1}{2} (E_{1111} - E_{2222}), \quad d'_{12} = E_{1112} + E_{2212}. \qquad (76)$$

* For the harmonic part term, we can introduce the harmonic parts (in Kelvin notation)

$$[\mathbf{H}(\mathbf{E})] = \begin{bmatrix} H_{1111} & -H_{1111} & \sqrt{2}H_{1112} \\ -H_{1111} & H_{1111} & -\sqrt{2}H_{1112} \\ \sqrt{2}H_{1112} & -\sqrt{2}H_{1112} & -2H_{1111} \end{bmatrix}, \quad [\mathbf{H}(\mathbf{A})] = H^* \begin{bmatrix} 1 & -1 & 0 \\ -1 & 1 & 0 \\ 0 & 0 & -2 \end{bmatrix}, \qquad (76)$$

$$[Q \star \mathbf{H}(\mathbf{A})] = H^* \begin{bmatrix} \cos 4\theta & -\cos 4\theta & \sqrt{2} \sin 4\theta \\ -\cos 4\theta & \cos 4\theta & -\sqrt{2} \sin 4\theta \\ \sqrt{2} \sin 4\theta & -\sqrt{2} \sin 4\theta & -2 \cos 4\theta \end{bmatrix},$$

where

$$H_{1111} = \frac{1}{8} (E_{1111} - 2E_{1122} - 4E_{1212} + E_{2222}), \quad H_{1112} = \frac{1}{2} (E_{1112} - E_{2212}).$$

We have then, by Eq. (11),

$$d_{\mathcal{O}_{rt}}(\mathbf{E})^2 = \min_{\theta, d^*, H^*} (\|\mathbf{d}'(\mathbf{E}) - Q \star \mathbf{d}'(\mathbf{A})\|^2 + \|\mathbf{H}(\mathbf{E}) - Q \star \mathbf{H}(\mathbf{A})\|^2). \qquad (77)$$

The minimizations with respect to $d^*$ and to $H^*$ give

$$d^* = d'_{11} \cos 2\theta + d'_{12} \sin 2\theta, \quad H^* = H_{1111} \cos 4\theta + H_{1112} \sin 4\theta. \qquad (78)$$

The minimization with respect to $\theta$ ends then up to the equation

$$A \cos 8\theta + B \sin 8\theta + C \cos 4\theta + D \sin 4\theta = 0, \qquad (79)$$

where (Antonelli et al., 2022):

$$\begin{align}
A &= 2(E_{1111} - 2E_{1122} - 4E_{1212} + E_{2222})(E_{1112} - E_{2212}) = 32H_{1111}H_{1112}, \nonumber \\
B &= 4(E_{1112} - E_{2212})^2 - \frac{1}{4}(E_{1111} - 2E_{1122} - 4E_{1212} + E_{2222})^2 = 16(H_{1112}^2 - H_{1111}^2), \nonumber \\
C &= 2(E_{1111} - E_{2222})(E_{1112} + E_{2212}) = 4d'_{11}d'_{12}, \nonumber \\
D &= 2(E_{1112} + E_{2212})^2 - \frac{1}{2}(E_{1111} - E_{2222})^2 = 2(d'_{12}^2 - d'_{11}^2). \tag{80}
\end{align}$$
Setting $\theta = \frac{1}{2} \arctan t$, we obtain the fourth-order polynomial in $t$,

$$ (A - C)t^4 + (2D - 4B)t^3 - 6At^2 + (2D + 4B)t + A + C = 0. \quad (81) $$

The roots $t_k$ of the polynomial are obtained by symbolic resolution of the polynomial using (Meurer et al., 2017), and then evaluated numerically. The root retained corresponds to either to the solution $\theta = \frac{1}{2} \arctan(t_k)$ or to the solution $\theta = \frac{1}{2} \arctan(t_k) - \frac{\pi}{2}$ that minimizes $d_{\theta r1}(E)^2$. Finally, the distance to orthotropy is obtained by injecting the minimizers $-d^*$ and $H^*$ both from Eq. (78), and $\theta$— and taking the square root of Eq. (77).

## B. Damaging loadings

Table 8 and Table 9 provide a description of the damaging loadings. Each loading is discretized into load steps. During KUBC and PBC loadings, the strain is imposed on the whole boundary $\partial\Omega$ of the Area Element. During EXPE loading, a displacement is imposed on the sub-parts of the boundary $\partial\Omega$. Each load step adds an increment of strain $\Delta\varepsilon(\mathbf{x})$, and displacement $\Delta\mathbf{u}(\mathbf{x})$, to the currently imposed strain, and displacement, respectively. The strain at load step $i$ is

$$ \begin{cases} \varepsilon^0 = \mathbf{0}, \\ \varepsilon^i = \varepsilon^{i-1} + \Delta\varepsilon^i. \end{cases} \quad (82) $$

The displacement at load step $i$, on a sub-part of the boundary $\partial\Omega_u$, is

$$ \forall \mathbf{x} \in \partial\Omega_u, \quad \begin{cases} \mathbf{u}^0(\mathbf{x}) = \mathbf{0} \\ \mathbf{u}^i(\mathbf{x}) = \mathbf{u}^{i-1}(\mathbf{x}) + \Delta\mathbf{u}^i(\mathbf{x}). \end{cases} \quad (83) $$

In practice, the boundary $\partial\Omega_u = \{\mathbf{x} = (x, y) | x = 0\}$ corresponds to the layer of particle having their center in the rectangle defined by $[0, \bar{t}_b] \times [0, L]$.

**Remark 9.** The rotation of particle is never imposed in those loadings.
<table><caption>**Table 8**</caption><thead><tr><th rowspan="2">Type</th><th rowspan="2">Name</th><th colspan="2">Load steps (i)</th><th colspan="3">Strain incs Δε<sup>i</sup> × 10<sup>6</sup></th></tr><tr><th>Start</th><th>End</th><th>ε<sub>xx</sub></th><th>ε<sub>yy</sub></th><th>ε<sub>xy</sub></th></tr></thead><tbody><tr><td rowspan="10">KUBC</td><td>bi-tension</td><td>1</td><td>100</td><td>6</td><td>6</td><td>0</td></tr><tr><td>shear</td><td>1</td><td>100</td><td>0</td><td>0</td><td>15</td></tr><tr><td>tension</td><td>1</td><td>100</td><td>10</td><td>0</td><td>0</td></tr><tr><td>shear → tension</td><td>1<br>50</td><td>50<br>100</td><td>0<br>20</td><td>0<br>0</td><td>10<br>0</td></tr><tr><td>tension → shear</td><td>1<br>50</td><td>50<br>100</td><td>10<br>0</td><td>0<br>0</td><td>0<br>20</td></tr><tr><td>tension → tension</td><td>1<br>50</td><td>50<br>100</td><td>10<br>0</td><td>0<br>10</td><td>0<br>0</td></tr><tr><td>(Willam et al., 1989)</td><td>1<br>100</td><td>100<br>200</td><td>4<br>2</td><td>0<br>1</td><td>0<br>3</td></tr><tr><td>bi-tension</td><td>1</td><td>100</td><td>5</td><td>5</td><td>0</td></tr><tr><td>shear</td><td>1</td><td>100</td><td>0</td><td>0</td><td>10</td></tr><tr><td>tension</td><td>1</td><td>100</td><td>3</td><td>0</td><td>0</td></tr><tr><td>shear → tension</td><td>1<br>50</td><td>50<br>100</td><td>0<br>40</td><td>0<br>0</td><td>10<br>0</td></tr><tr><td>tension → shear</td><td>1<br>50</td><td>50<br>100</td><td>10<br>0</td><td>0<br>0</td><td>0<br>20</td></tr><tr><td>tension → tension</td><td>1<br>50</td><td>50<br>100</td><td>6<br>0</td><td>0<br>6</td><td>0<br>0</td></tr><tr><td>(Willam et al., 1989)</td><td>1<br>100</td><td>100<br>200</td><td>8<br>4</td><td>0<br>6</td><td>0<br>2</td></tr></tbody></table>

KUBC and PBC damaging loadings. For PBC loadings, an particle in locked in displacement to prevent rigid body motions.

<table><caption>**Table 9**</caption><thead><tr><th rowspan="3">Name</th><th colspan="2">Load steps</th><th colspan="6">Displacement incs Δu'(x) × 10<sup>6</sup></th></tr><tr><th rowspan="2">Start</th><th rowspan="2">End</th><th colspan="2">x = 0</th><th colspan="2">x = L</th><th colspan="2">y = L</th></tr><tr><th>u<sub>x</sub></th><th>u<sub>y</sub></th><th>u<sub>x</sub></th><th>u<sub>y</sub></th><th>u<sub>x</sub></th><th>u<sub>y</sub></th></tr></thead><tbody><tr><td>bitension</td><td>1</td><td>100</td><td>0</td><td>-</td><td>1.2</td><td>-</td><td>-</td><td>0</td><td>-</td><td>1.2</td></tr><tr><td>tension</td><td>1</td><td>100</td><td>0</td><td>-</td><td>1</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>shear</td><td>1</td><td>100</td><td>-</td><td>0</td><td>-</td><td>0.8</td><td>0</td><td>-</td><td>0.8</td><td>0</td></tr><tr><td>simple shear</td><td>1</td><td>100</td><td>-</td><td>-</td><td>-</td><td>-</td><td>0</td><td>0</td><td>0</td><td>4</td></tr><tr><td>shear → tension</td><td>1<br>50</td><td>50<br>100</td><td>0<br>0</td><td>-</td><td>2</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>tension → shear</td><td>1<br>50</td><td>50<br>100</td><td>0<br>-</td><td>-</td><td>1</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>tension → tension</td><td>1<br>50</td><td>50<br>100</td><td>0<br>-</td><td>-</td><td>1.2</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr></tbody></table>

EXPE loadings. For most of those loadings, an additional constraint is applied to a particle to prevent rigid body motion.