# Micromechanical Estimates Compared to FE-Based Methods for Modelling the Behaviour of Micro-Cracked Viscoelastic Materials

Sarah Abou Chakra, Benoît Bary, Eric Lemarchand, Christophe Bourcier, Sylvie Granet, Jean Talandier

► To cite this version:

Sarah Abou Chakra, Benoît Bary, Eric Lemarchand, Christophe Bourcier, Sylvie Granet, et al. Micromechanical Estimates Compared to FE-Based Methods for Modelling the Behaviour of Micro-Cracked Viscoelastic Materials. Modelling, 2024, 5 (2), pp.625-641. 10.3390/modelling5020033 hal-04680472

HAL Id: hal-04680472

https://cnrs.hal.science/hal-04680472v1

Submitted on 20 Nov 2024

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 Attribution 4.0 International License
# Micromechanical Estimates compared to FE-based Methods for Modelling the Behaviour of Micro-cracked Viscoelastic Materials

Sarah Abou Chakra ¹*, Benoît Bary ², Eric Lemarchand ³, Christophe Bourcier ⁴, Sylvie Granet ⁵ and Jean Talandier ⁶

¹ Service de Recherche en Corrosion et Comportement des Matériaux, Université Paris-Saclay, CEA, France ; sarah.abouchakra@cea.fr 6

² Service de Recherche en Corrosion et Comportement des Matériaux, Université Paris-Saclay, CEA, France ; benoit.bary@cea.fr 8

³ Laboratoire Navier, Ecole des Ponts ParisTech, France ; eric.lemarchand@enpc.fr 10

⁴ Service de Génie Logiciel pour la Simulation, Université Paris-Saclay, CEA, France ; christophe.bourcier@cea.fr 11

5 Département ERMES, EDF - R&D, France, sylvie.granet@edf.fr 13

⁶ ANDRA, France, jean.talandier@andra.fr 14

* Correspondence: sarah.abouchakra@cea.fr 15

**Abstract:** The purpose of this study is to investigate the effective behaviour of a micro-cracked material whose matrix bulk and shear moduli are ruled by a linear viscoelastic Burgers model. The analysis includes a detailed study of randomly oriented and distributed cracks displaying an overall isotropic behavior, as well as aligned cracks resulting in a transversely isotropic medium. Effective material properties are approximated with the assumption that the homogenized equivalent medium exhibits the characteristics of a Burgers model, leading to the identification of short-term and long-term homogenized modules in the Laplace-Carson space through simplified formulations. The crucial advantage of this analytical technique consists in avoiding calculations of the inverse Laplace-Carson transform. The micromechanical estimates are validated through comparisons with FE numerical simulations on 3D microstructures generated with zero-thickness void cracks of disk shape. Intersections between randomly oriented cracks are accounted for, thereby highlighting potential percolation phenomenon. The effects of micro-cracks on the material behavior are then studied with the aim of providing high-performance creep models for macrostructure calculations at a moderate computation cost through the application of analytical homogenization techniques.

**Keywords:** Cracked media; Viscoelasticity; Burgers model; 3D numerical simulation 30

**Citation:** To be added by editorial staff during production.

Academic Editor: Firstname Last-name

Received: date

Revised: date

Accepted: date

Published: date

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

**Copyright:** © 2024 by the authors.
Submitted for possible open access
publication under the terms and
conditions of the Creative Commons
Attribution (CC BY) license
(https://creativecommons.org/licenses/by/4.0/).

## 1. Introduction 32

In the context of long-term deep disposal of nuclear waste in the Callovo-Oxfordian (COx) claystone, the excavation of tunnels in the host rock inevitably leads to macroscopic fracturing, raising potentially significant challenges for nuclear waste management [1–3]. Characterizing the behavior of the COx claystone is complex, as various factors are involved in determining the appropriate behaviour of the claystone. In particular, the creep properties of the argillite play a crucial role in the crack sealing process [4], an anticipated key phase expected to take place following the closure of the underground waste disposal [5]. This study aims to provide a description for the creep phenomenon within a simplified linear viscoelastic framework. Furthermore, it is intended to model the behavior of a fractured viscoelastic material by approximating its effective properties and identifying the parameters by means of a macroscopic Burgers model [6]. For the purpose of validating the analytical approach, the results are compared with numerical simulations per-
formed using the finite element method with the Cast3M<fnref n="1" /> software. This comparison between the two approaches offers key validation requirements when it comes to modelling the macroscopic behavior of damaged materials while exploiting the use of multi-scale approaches. The excavation of tunnels for the Cigéo project causes fracturing of the COx claystone. Characterizing the cracking state through in situ observations provides means for dimensioning the damaged zone and determining cracks orientations as a function of distance from the excavated drifts. In the vicinity of the galleries, cracks orientations appear widely variable and random. With increasing distance from the drifts, cracks show preferential orientations mostly parallel to the galleries [2].

Several techniques are employed to accurately describe the complex phenomenon of fracture in materials [7–11]. Alongside damage models, which incorporate the effect of cracks by introducing a damage variable [12], methods based on continuum mechanics can also capture the presence of cracking in a solid material [13]. In particular, micromechanical techniques based on the calculation of homogenized moduli can be used to estimate the effect of fractures on the material macroscopic properties [14–16]. Cracks are considered as heterogeneities distributed in the solid matrix, and solutions based on micromechanical estimates are presented to deduce homogenized properties. The approach used by Budiansky and O'Connell [17] involves the utilization of a crack density parameter to characterize the fracturing effect, assuming that all cracks are circular with an identical shape and size. The crack density parameter can be associated macroscopically with a damage variable, whose effects are calculated from analytical homogenization schemes [18]. In addition, other techniques, notably finite element modeling [19], can be used to create microstructures representative of the fracturing condition, with considerably fewer limitations in terms of design simplifications. This numerical approach can be directly compared to analytical techniques obtained by approximation schemes, which is a key objective of this study.

Materials subjected to constant loads can exhibit a delayed (creep) behavior, characterized by continuous deformation over time. Linear viscoelastic models can depict such effects in a simplified context [20]. Moreover, it is possible to use the Laplace-Carson transform (LC) and proceed to operate with an equivalent linear elastic behavior in the Laplace-Carson space. This provides a straightforward framework to recourse on analytical solutions obtained with classical homogenization procedures developed for elastic materials affected by the presence of cracks and/or heterogeneities [21]. For this purpose, the behavior of the COx claystone is represented here by non-aging linear viscoelastic Burgers models combining Maxwell and Kelvin-Voigt models. Using a minimum number of parameters, the Burgers model is indeed able to capture the main features of a viscoelastic material behavior, particularly providing an adequate representation of the long-term characteristics of the COx claystone. In linear viscoelastic applications, dealing with heterogeneous media considerably complexifies the formulated equations in LC space. This difficulty renders model equation inversion into real time space challenging. Commonly, numerical solutions, such as the collocation method, are used to provide effective resolutions to the equations describing the material's behavior [22,23].

The most classical approaches used to estimate the behavior of heterogeneous materials with homogenization techniques are based on the Eshelby solution [24]. These approaches have proved to be very efficient for approximating the effective behavior of such materials in various contexts, including for micro-cracked linear elastic materials [19,25]. In this case, an alternative approach has been proposed in [26], based on the linear relationship between macroscopic stress and the discontinuity of local displacements across the crack. This approach will be explored in the present study, together with an approximation method assuming that the effective behavior of the material can be represented by a set of Burgers models. This efficient method provides an analytical solution of the model

<footnote n="1">¹ http://www-cast3m.cea.fr/</footnote>
equations in simplified forms, making inversion calculations analytically achievable. One key advantage is that their parameters can be directly expressed as a function of the parameters of the viscoelastic models describing the sound material, and the crack density parameter characterizing the cracking state. This differs from more classical approaches based on the application of homogenization schemes in the LC space, which generally requires a numerical inversion to characterize the model in the time space. The analytical results will be compared to finite element simulations obtained on 3D samples generated with disk-shaped voids of zero thickness, thus mimicking very accurately the behavior of micro-cracked materials. In this context, the material's behavior will be analyzed by considering relevant cases of crack orientations including randomly oriented cracks as well as a parallel crack distribution in an isotropic medium.

## 2. Viscoelastic Burgers model

The majority of bedrock masses, when subjected to a permanent stress, exhibit to some degree a deferred behaviour, which is expressed by an increase in strain under constant loading. This phenomenon indicates that a complete description of the constitutive behavior of a rock mass must assume time-dependency. A rigorous approach is to use viscoelastic or viscoplastic behavior laws to describe the evolution of creep phenomenon for geomaterials [4,27,28]. In analyzing the behavior of the Callovo-Oxfordian claystone, the effect of delayed long-term response of the material must be taken into account to ensure an accurate representation of the COx claystone [29]. In this context, we focus mainly on non-aging linear viscoelastic models, which allow establishing a linear relation between the stress and strain in the Laplace-Carson space using the correspondence principle. The development of the Burgers model equations provides an efficient representation of the material's short- and long-term behavior, while requiring a limited number of parameters. It is intended to use this model to represent the behavior of the COx claystone. The indices M and K designate respectively the Maxwell and Kelvin parts of the model, while exponents s and d are attributed to the spherical and deviatoric parts of the model (Figure 1). Moreover, we propose to use 2 different Burgers models to describe independently the shear and bulk moduli denoted by $\mu$ and $k$, respectively.

<figure><img src="image_3.png" /><figcaption>Figure 1. Schematic representation of the Burgers model</figcaption></figure>

$$ \frac{1}{k^*} = \frac{1}{k_M} + \frac{1}{p \frac{\eta_M^s}{3}} + \frac{1}{k_K + p \frac{\eta_K^s}{3}}; \quad \frac{1}{\mu^*} = \frac{1}{\mu_M} + \frac{1}{p \frac{\eta_M^d}{2}} + \frac{1}{\mu_K + p \frac{\eta_K^d}{2}} \quad (1) $$

Eq. 1 expresses the effective compressibility and shear moduli in Laplace-Carson space, denoted by $k^*$ and $\mu^*$ respectively, where p is the LC variable. The eight parameters of the Burgers Model can be distinguished in Eq. 1. The Maxwell part is represented by four parameters, including spherical components ($k_M, \eta_M^s$) and deviatoric components ($\mu_M, \eta_M^d$). Additionally, the Kelvin part of the model is also composed of four parameters, with ($k_K, \eta_K^s$) describing the spherical parameters and ($\mu_K, \eta_K^d$) representing the deviatoric
components. The parameters of the Burgers models are identified through an analytical resolution of the strain under constant stress loading, followed by comparisons with experimental creep testing results. Armand et al. [30] have published a series of creep tests widely used in the literature to characterize this phenomenon. Cores of COx were extracted from the Meuse/Haute-Marne underground research laboratory and the samples were drilled perpendicular to the bedding plane of the rock, as the COx claystone naturally exhibits a transversely isotropic behavior. These tests include different series of loadings and have been designed to examine the deferred behavior of the material for different loading conditions. The applied loadings are a function of the maximum deviatoric stress that the material can withstand. The argillite is considered anisotropic and therefore its properties change as a function of the imposed loading with respect to the bedding plane. Here, for commodity the solid material will be considered isotropic, meaning that the elastic moduli do not vary with respect to the bedding plane. The model parameters have been identified under loading conditions corresponding to 75% of the maximum deviatoric stress.

The use of numerical models with a number of parameters can render the precise identification of each parameter from a single experimental test complex. Although some simplifications can be achieved at the limits as time tends to zero or infinity, it remains difficult to identify each parameter individually. Consequently, optimization methods were employed to refine the model parameters. The least-squares method is used to optimize the Burgers model parameters using standard Python libraries (Figure 2). The normal and tangential components of the strain tensor ($\epsilon_{nn}$ and $\epsilon_{tt}$ respectively) in the time space are identified by calculating the inverse of the Carson Laplace transform of:

$$ \epsilon^* = \frac{1 + \nu^*}{E^*} \boldsymbol{\sigma}^* - \frac{\nu^*}{E^*} \mathrm{tr}(\boldsymbol{\sigma}^*) \boldsymbol{I} \quad (2) $$

Where the superscript * stands for a quantity expressed in the LC space. The normal and tangential strain components due to the application of a uniaxial constant stress $\sigma_0$ can be deduced in LC space from Eq. 2:

$$ \epsilon_{nn}^* = \frac{\sigma_0}{E^*} \quad \text{and} \quad \epsilon_{tt}^* = \frac{-\sigma_0 \nu^*}{E^*} \quad (3) $$

Inversion of the normal strain $\epsilon_{nn}^*$ from LC space reads:

$$ \epsilon_{nn} = \frac{\sigma_0}{E_M} + \frac{\sigma_0}{9k_M} \frac{t}{\left(\frac{\eta_M^s}{3k_M}\right)} + \frac{\sigma_0}{3\mu_M} \frac{t}{\left(\frac{\eta_M^d}{2\mu_M}\right)} + \frac{\sigma_0}{9k_K} \left(1 - e^{\frac{t}{\eta_K^s}}\right) + \frac{\sigma_0}{3\mu_K} \left(1 - e^{\frac{t}{\eta_K^d}}\right) \quad (4) $$
with $v_M = \frac{3k_M - 2\mu_M}{6k_M + 2\mu_M}$ and $\frac{1}{E_M} = \frac{1}{9k_M} + \frac{1}{3\mu_M}$

This expression involves the eight parameters of the Burgers model and the time.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Identification of the Burgers model parameters based on experimental creep test results</figcaption></figure>

Overall, the proposed model provides an acceptable representation of the argillite's behavior in both short and long term, see figure 2 for the comparison between experimental and numerical results of a creep test. The characteristic parameters of the Burgers model are summarized in the Table 1.

<table><caption>**Table 1.** Identified viscoelastic parameter of the COx claystone</caption><thead><tr><th>Parts</th><th>k(GPa)</th><th>μ(GPa)</th><th>η<sup>s</sup>(GPa.s)</th><th>η<sup>d</sup>(GPa.s)</th></tr></thead><tbody><tr><td>Maxwell</td><td>6.872</td><td>3.870</td><td>35.381 · 10<sup>7</sup></td><td>12.463 · 10<sup>7</sup></td></tr><tr><td>Kelvin</td><td>7.308</td><td>3.099</td><td>26.364 · 10<sup>6</sup></td><td>4.020 · 10<sup>6</sup></td></tr></tbody></table>

## 3. Representation of cracking in a viscoelastic material

In this work, we make use of the approach proposed in [6,26,31,32] for considering the presence of cracks, which is not directly based on Eshelby's classical approach [24], although being very close. The interested reader is invited to refer to the cited papers for more details that would be too lengthy to include here.

### 3.1. Isotropic crack distribution

Considering an isotropic distribution of cracks in a solid matrix (Figure 3), where cracks are characterized by the crack density parameter $\epsilon$ in which $N$ is the number of cracks per unit volume of REV (Representative Elementary Volume) and $a$ is the crack radius, assuming an identical circular shape and a uniform size to all cracks (Eq. 5).

$$ \epsilon = Na^3 \qquad (5) $$

The apparent effective stiffness tensor in the Laplace-Carson space defined by $C^{hom*}$ is therefore isotropic, expressed in the $(J, K)$ basis in (Eq. 6). The tensors $J$ and $K$ are respectively the projectors of the spherical and deviatoric parts of the 4th-order unit tensor.

$$ C^{hom*} = 3k^{hom*}J + 2\mu^{hom*}K \qquad (6) $$
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Schematic representation of randomly oriented crack distribution in a solid matrix</figcaption></figure>

The homogenized moduli $k^{hom*}$ and $\mu^{hom*}$ are calculated as a function of the crack density parameter $\epsilon$, taking into account two distinct loading conditions: spherical loading which is used to determine the compressibility modulus, while deviatoric loading is used to determine the shear modulus. A homogenized stiffness matrix is then constructed in the LC space, reflecting the material cracking conditions [26].

$$ \frac{1}{k^{hom*}} = \frac{1}{k^*} + \frac{8\epsilon}{3} \frac{(1 - \nu^*)}{\mu^*} \qquad (7) $$

$$ \frac{1}{\mu^{hom*}} = \frac{1}{\mu^*} + \frac{32\epsilon}{\mu^*} \frac{(1 - \nu^*)(5 - \nu^*)}{45(2 - \nu^*)} \qquad (8) $$

To simplify, here the homogenized material behavior in the LC space is studied by considering a general analogy with a homogenized macroscopic Burgers model, i.e. we assume that the behavior of the homogenized material can be approached by Burgers models. This makes it possible to identify the homogenized modules $k^{hom*}$ and $\mu^{hom*}$ in the LC space using a simplified form:

$$ \frac{1}{k^{hom*}} = \frac{1}{k_M(\epsilon)} + \frac{1}{p \eta_M^s(\epsilon)/3} + \frac{1}{k_K(\epsilon) + p \eta_K^s(\epsilon)/3} \qquad (9) $$

$$ \frac{1}{\mu^{hom*}} = \frac{1}{\mu_M(\epsilon)} + \frac{1}{p \eta_M^d(\epsilon)/2} + \frac{1}{\mu_K(\epsilon) + p \eta_K^d(\epsilon)/2} \qquad (10) $$

This maintains the same structure of the Burgers model formulated in Eq. 1. Proceeding with a Taylor expansion of the expressions of the moduli with respect to the variable $p$ in the Carson space in the vicinity of 0 and infinity, the eight parameters of the Burgers models can be identified as a function of the crack density parameter (Figure 4). A correlation of the parameters as a function of crack density (Eq. 5) is established, as illustrated in Figure 4. In summary, this approach provides an alternative method for characterizing the behavior of a micro-cracked material in the Laplace-Carson space [6].

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Effective moduli calculation as a function of crack density, considering an isotropic crack distribution in the solid matrix</figcaption></figure>

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204
Determining the parameter equations as a function of crack density for a macroscopically cracked material of the Burgers type enables a straightforward formulation of the homogenized moduli that define this problem. This approach can provide an explicit formulation of the effective Burgers model equations in the time space as a function of the parameters of the undisturbed material and the crack density parameter (Eq. 9 and 10). This significantly reduces the complexity of the calculation of the modules in the time-space, as it avoids the direct application of the inverse of the Laplace-Carson transform on such equations. This in turn leads to the development of straightforward analytical solutions, offering results comparable to those obtained by alternative numerical approaches such as the finite element method [19] or the FFT approach, which also enable an analysis of the effective moduli on microstructures [7,33].

## 3.2. Parallel crack distribution

In this section, the effect of a random distribution of parallel cracks in a linear visco-elastic material as described in the previous sections, will be discussed. In this case, a transversely isotropic behavior is obtained (Figure 5). In this context, it is convenient to work with the expression of the material compliance tensor in LC space $S^*(p)$ (Eq. 11), which is the inverse of the stiffness tensor using the Walpole basis [34]. Using the Walpole basis simplifies the expression of the viscoelastic compliance tensor in the case of anisotropy, making the analysis of modules affected by the fracturing more convenient.

$$S^*(p) = \sum_{i=1}^{6} s_i^*(p)E_i \qquad (11)$$

with:

$$ \begin{aligned} s_1^*(p) &= \frac{1 + v^*}{2\mu^*(1 + v^*)} ; & s_2^*(p) &= \frac{1}{2\mu^*(1 + v^*)} ; & s_3^*(p) &= \frac{1}{2\mu^*} \\ s_4^*(p) &= \frac{1}{2\mu^*} ; & s_5^*(p) &= s_6^*(p) = -\frac{v^*}{2\mu^*(1 + v^*)} \end{aligned} \qquad (12) $$

By decomposing the compliance tensor in the Walpole basis (Eq. 11), we obtain six scalars with indexes i ranging from 1 to 6 denoted $s_i^*(p)$. The tensors $E_i$ are defined in accordance with the specifications given in the Walpole basis [34]. Eq. 12 describes the expressions of the scalars $s_i^*(p)$ in the Laplace-Carson space. Note that the directions $(e_1, e_2, e_3)$ respectively corresponds to the conventional directions in a Cartesian coordinate system $(x, y, z)$. $n$ is the direction normal to the fractures along $e_1$ (Figure 5).

<figure><img src="image_7.png" /><figcaption>Figure 5. Parallel crack distribution in the solid matrix</figcaption></figure>
As demonstrated e.g. in [26], the calculations highlight that for a transversely isotropic medium, some of the modules are not influenced by the cracking state of the material. This observation has been confirmed in an analysis of a linear elastic behaviour case, where only the Young's modulus perpendicular to the bedding plane and the shear coefficients were affected by the cracking [35]. Following the notation for the scalars indexes in the Walpole's basis, the expression for the moduli $s_2^{hom*}$ and $s_4^{hom*}$ as a function of the cracking density yields:

$$s_2^{hom*} = s_{M2}^e(\epsilon) + \frac{1}{p} s_{M2}^v(\epsilon) + (s_{K2}^e(\epsilon)^{-1} + p s_{K2}^v(\epsilon)^{-1})^{-1} \quad (13)$$

$$s_4^{hom*} = s_{M4}^e(\epsilon) + \frac{1}{p} s_{M4}^v(\epsilon) + (s_{K4}^e(\epsilon)^{-1} + p s_{K4}^v(\epsilon)^{-1})^{-1} \quad (14)$$

provided that $s_2^{hom*}$ and $s_4^{hom*}$ are the two moduli affected by the cracking [26]. The exponents e and v are attributed to the elastic and viscous parts, respectively. The model parameters are determined based on the assumption that the macroscopic model is essentially of the Burgers type, allowing the parameters to be identified with simplified analytical expressions, analogous to the isotropic case discussed in the previous section. The determination of $s_2^{hom*}$ is obtained by imposing a uniaxial tensile load along the normal to the cracks, while that of $s_4^{hom*}$ requires a shear loading in the form $\sigma (e_2 \otimes n + n \otimes e_2)$, e.g. with $\sigma = 1MPa$. As the other modules will not be affected by the cracks, it is possible to attribute to them the modules of the undisturbed material without requiring further calculations.

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Determination of the moduli affected by cracking for the transversely isotropic case in the Walpole basis as a function of crack density</figcaption></figure>

The same identification technique as for the isotropic case was used, considering that the homogenized behavior of the cracked material corresponds to a macroscopic Burgers model, while accounting for the crack density parameter. Conserving the same structure of the Burgers model formulation, it is possible to define a correlation of the 8 coefficients in Eq. 13 and 14 as a function of crack density (Eq. 5) as shown in figure 6. The moduli affected by cracking in a Cartesian coordinate system, such as $E_1^{hom}$ and $\mu_{12}^{hom}$, can be determined from calculations in the Walpole basis as a function of the coefficients $s_2^{hom}$ and $s_4^{hom}$ [34]. As mentioned above, investigating the transversely isotropic case is representative of the in situ behavior of the fractured COx claystone, since in reality, fracture formation in parts of the damaged zone is produced with a preferential orientation parallel to the excavated drifts [2]. A more in-depth analysis of the analytical results is carried out by comparing the results with the numerical simulations obtained with the Cast3M finite element code.

# 4. Numerical simulations

## 4.1. 3D microstructures and meshes
The Burgers model is implemented in the Mfront<fnref n="2" /> library, offering several advantages and convenient access for Cast3M users. Mfront is a code generator that enables convenient implementation of behavioral laws, which can be efficiently linked to Cast3M to perform finite element calculations [36]. The parameters of the Burgers model have been calibrated in the previous section. A simple uniaxial creep test with a loading of $\Sigma_{11} = 1$ MPa is performed to verify the perfect agreement between analytical and numerical results for the case where a homogeneous isotropic material is portrayed (Figure 7).

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Cast3M verification of the conformity between the implemented Burgers model and the analytical results, obtained following a creep test using the undisturbed material</figcaption></figure>

The 3D sample generation procedure is based on a random distribution of inclusions of prescribed shape and size in a box [37,38]. Disk-shaped cracks of zero thickness were introduced as they are suitable for direct comparisons with analytical results [7,39,40]. From this viewpoint, this generation procedure constitutes an improvement with respect to the previous method in which the inclusions are constituted of volume inclusions of non-zero thickness [19]. The geometry of the REV is periodic, which implies that periodic boundary conditions can be applied. The 3D meshes with tetrahedral elements were generated according to the selected periodic geometries using an automatic software connected to Salome<fnref n="3" /> (Figure 8). In this study, we use the MG-CADSurf algorithm since it allows us to generate periodic meshes. A mesh refinement is imposed systematically on the surfaces of the cracks by increasing the density of elements, to improve the quality of the simulation results. The technique employed to create the cracks is intended to be used with specific interface elements of zero thickness, i.e. double-nodes are created at each point of the cracks mesh. Here as the cracks are void, no interface elements are introduced, allowing to drastically reduce the calculation time at the REV scale [7,40].

<footnote n="2">² https://tfel.sourceforge.net/gallery.html</footnote>

<footnote n="3">³ https://www.salome-platform.org</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Microstructure meshes generation in Salome with a crack density equal to 0.1 with 100 cracks distributed in the solid matrix (ratio of crack radius to cube side length is 20/200), considering a parallel crack distribution (left) and an isotropic crack distribution (right) in the REV</figcaption></figure>

A series of test cases were examined featuring an increasing range of crack density parameters, and comparisons were subsequently conducted with analytical results. Starting with relatively low crack density parameters, the 4 crack density parameters investigated are 0.05, 0.1, 0.15 and 0.2. A progressively increasing number of cracks has been used, in order to keep a constant size of the micro-cracks in the REV. Thus, for the corresponding testing cases, the number of cracks used is 50, 100, 150 and 200, respectively. This is applicable to both isotropic and transversely isotropic cases. In order to achieve a consistent isotropic behavior of the microstructure, the "Fibonacci sphere" method is used, defining a uniform distribution of points on a sphere [41]. The number of these points is equivalent to the number of cracks considered, and the uniformly distributed points on the sphere serve to define the crack orientations in the microstructure (Figure 9). The placement of the cracks however remains completely random in all considered cases.

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Even distribution of points across the sphere in order to define crack orientations in the microstructure and ensure an isotropic behavior</figcaption></figure>

Sufficient numbers of micro-cracks must be included in the solid matrix to ensure a good representativeness, which is why a comprehensive study must be carried out prior to calculations to minimize errors attributable to the REV generation.
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Evaluation of REV representativity by varying the number of cracks distributed in the microstructure while keeping a constant crack density $\epsilon = 0.2$. The effective Young's modulus in the direction normal to the cracking plane for the transverse isotropic case are calculated (left) as well as the effective Young's modulus for the isotropic case (right) and periodic stress boundary conditions were applied</figcaption></figure>

The analysis of Figure 10 demonstrates that the microstructures generated are fairly representative, and with minimized deviations, based on the calculation of the effective Young's modulus as a function of time, whereas the number of cracks involved is varied from 100 to 200 cracks while the crack density parameter is constant at 0.2. An exponential evolution of the time steps was considered, using 100 time steps in total, ensuring a high concentration of time steps at the early stages of the calculation where there is the greatest variation in the material's behavior. The loading applied is a stress of $\sigma = 1MPa$ with $\Sigma = \sigma \cdot e_1 \otimes e_1$, constant with time. Simulations last around 2 hours on average in these cases. The number of elements in the meshes in each case majorly affects the simulation time. A mesh convergence study must be carried out to ensure that a sufficient number of mesh elements is used to guarantee optimized results with minimized deviations. In this case, convergence is achieved with a relatively low error of no more than 1% in cases where the number of tetrahedral mesh elements averages 1.5 million. This is shown by calculating the homogenized Young and shear modulus considering microstructures with an increasing numbers of elements. Consequently, increasing the number of elements beyond this level has no significant effect on the investigated mechanical properties. The calculations are illustrated for the last time step of 100 days (Table 2 and Table 3), with periodic boundary conditions and material subjected to a constant stress of $\sigma = 1MPa$. Two simulations are required to calculate the effective properties: one with $\Sigma = \sigma \cdot e_1 \otimes e_1$ for the calculation of the homogenized Young's modulus, and the second one with $\Sigma = \sigma \cdot (e_1 \otimes e_2 + e_2 \otimes e_1)$ to determine the homogenized shear modulus. For all subsequent simulations, a mesh size of around 1.5 million elements was used. It's important to note that the number of elements is strongly related to the number of cracks, and technically to the crack density, given that the mesh refinement is imposed to accurately describe the cracks.

<table><caption>**Table 2.** Mesh convergence study considering an isotropic behavior: computation of the homogenized Young's modulus and shear modulus. The relative error with respect to the last column is given in parentheses</caption><thead><tr><th></th><th>Mesh 1</th><th>Mesh 2</th><th>Mesh 3</th><th>Mesh 4</th></tr></thead><tbody><tr><td>Number of mesh elements</td><td>0.614 M</td><td>1.115 M</td><td><strong>1.602 M</strong></td><td>2.55 M</td></tr><tr><td>E<sup>hom</sup> (GPa)</td><td>2.9877 (2.45%)</td><td>2.8543 (1.78%)</td><td><strong>2.8941 (0.41%)</strong></td><td>2.9063</td></tr><tr><td>μ<sup>hom</sup> (GPa)</td><td>1.2299 (4.39%)</td><td>1.2066 (2.41%)</td><td><strong>1.1578 (1.71%)</strong></td><td>1.1781</td></tr></tbody></table>
<table><caption>**Table 3.** Mesh convergence study considering a transversely isotropic behavior: computation of the effective out-of-plane Young's modulus $E_1^{hom}$ and shear modulus $\mu_{12}^{hom}$. The relative error with respect to the reference case (last column) is calculated.</caption><thead><tr><th></th><th>Mesh 1</th><th>Mesh 2</th><th><strong>Mesh 3</strong></th><th>Mesh 4</th></tr></thead><tbody><tr><th>Number of mesh elements</th><td>0.675 M</td><td>0.919 M</td><td><strong>1.88 M</strong></td><td>2.88 M</td></tr><tr><th><i>E</i><sub>1</sub><sup>hom</sup> (GPa)</th><td>1.3627 (6.71%)</td><td>1.4053 (3.79%)</td><td><strong>1.4493 (0.78%)</strong></td><td>1.4607</td></tr><tr><th><i>μ</i><sub>12</sub><sup>hom</sup> (GPa)</th><td>0.8918 (2.27%)</td><td>0.9054 (0.77%)</td><td><strong>0.9089 (0.38%)</strong></td><td>0.9125</td></tr></tbody></table>

The specific nature of imposed boundary conditions can also influence the numerical results, and in this case three boundary conditions were tested: uniform displacement boundary conditions (KUBC), uniform stress boundary conditions (SUBC) and finally periodic boundary conditions (PBC). Modules are calculated by applying a constant load of stress equal to 1 MPa and the full load is applied from the very first time step. As expected, the calculation results reported on Figure 11 show that the most rigid response is obtained when imposing uniform displacement boundary conditions, while the softest response is obtained when imposing homogeneous stress boundary conditions. Considering that periodic boundary conditions yield an intermediate response both when testing the isotropic and transversely isotropic material, one assume that they provide the most accurate results and they will be used as a reference for comparison with analytical results.

<figure><img src="image_13.png" /><figcaption>**Figure 11.** Effects of the applied boundary conditions KUBC, SUBC and PBC on the calculated homogenized Young's modulus (a) and shear modulus (b) for the isotropic medium, and on the Young's modulus in the direction normal to the crack (c) and the effective shear modulus (d) for the transverse isotropic case, as a function of time</figcaption></figure>

348

349

350

351

352

353

354

355

356

357

358

359

360

361

362

363

364

365

366

367

368
## 4.2. *Transversely isotropic case results*

Starting with the transversely isotropic medium and applying periodic stress boundary conditions, a uniaxial creep test is carried out under constant normal stress $\Sigma_{11} = 1$ MPa. The time steps intervals are identical to the ones defined in section 4.1. At $t = 0$, the initial conditions correspond to the REV assumed at rest and without any loading. At $t = 0 +$, the full loading is applied. Again, the size of REV meshes is around 1.5 million elements, consistent with the mesh convergence study presented in section 4.1. Comparisons are subsequently conducted between numerical and analytical methods through the calculation of homogenized moduli and the determination of the normal strain as a function of time. In the present case, the moduli affected by cracking are Young's modulus $E_{11}$ and shear modulus $\mu_{12} = \mu_{13}$, considering the case of parallel cracks where the Ox axis is the direction perpendicular to the corresponding fractures plane in a Cartesian reference coordinate system. The normal strain is calculated as a function of time for the considered crack density parameter (Figure 12). The effective moduli are determined from two separate simulations: the first involving an axial loading at constant stress $\Sigma_{11} = 1$ MPa to calculate the effective Young's modulus, and the second requiring a deviatoric loading at $\Sigma_{12} = 1$ MPa to assess the effective shear modulus (Figure 12). The analytical approach however involves calculating the inverse of the Laplace transform, and thus returning to temporal space. The inverse of the Laplace-Carson transform can be derived directly from the simplified expressions outlined in the previous sections, since it follows a classical form of the Laplace transforms, which can be analytically computed. The analytical results are deduced from the equations presented in section 3.2. For the transverse isotropic case, homogenized moduli are calculated using Walpole's basis (Eq. 13 and 14), which is used to determine subsequently the Young's modulus and the homogenized shear modulus.

<figure><img src="image_14.png" /><figcaption>**Figure 12.** Comparison between numerical and analytical results with $\epsilon = 0.05$ and a number of cracks equal to 50, the strain under constant stress loading is calculated (a) as well as the homogenized moduli affected by the cracking for transverse isotropic medium (b) and (c)</figcaption></figure>

Initial results show that for the relatively low crack density parameter equal to 0.05, whereby the material still retains a large proportion of its mechanical strength, the numerical and analytical approaches produce very similar results for this case (Figure 12). However, upon investigating the cases with greater crack density parameter, the numerical results gradually start to deviate from the analytical estimations for the transverse isotropic case, as shown in Figure 13 which depicts the evolution of the creep strain as a function of time for crack density parameters of 0.1, 0.15 and 0.2.
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Creep strain calculation featuring numerical-analytical comparisons with an increasing crack density of 0.1 (a), 0.15 (b) and 0.2 (c)</figcaption></figure>

In the final case where the crack density parameter is equal to 0.2, it becomes clear that the numerical results are no longer comparable with the analytical ones (Figure 13 c)). In this situation, a different and more appropriate analytical technique is required since the previously described analytical method is only suitable for low crack density parameters. In this regard, previous studies have demonstrated that, in the case of parallel cracks within a linear elastic framework, the differential scheme gives very good results when compared to numerical calculations using the finite element method [35]. This aspect will be investigated in a future contribution.

### 4.3. *Isotropic case results*

We analyze in this section the creep response of the material when the cracks are randomly oriented, assuming periodic stress boundary conditions. Identical simulation conditions in this section are used to determine the properties of an isotropic medium where the only difference relative to the transverse isotropic case presented in section 4.1 lies in the definition of cracks orientations. The first test is based on a REV with a defined crack density parameter of 0.05 (Figure 14), and the numerical results in terms of normal strain and effective moduli evolutions show good agreement with the analytical ones, similarly to the transverse isotropic results. The analytical solution was obtained using the simplified equations for homogenized moduli in Laplace-Carson space (Eq. 9 and 10), enabling direct calculation of the homogenized Young's modulus and shear modulus in real time space. The normal strain $\epsilon_{nn}$ is calculated as a function of the Young's modulus $E^{hom}$ and the applied constant stress $\sigma = 1MPa$ in LC space (Eq. 3).

<figure><img src="image_16.png" /><figcaption>**Figure 14.** Comparison between numerical and analytical results with $\epsilon = 0.05$, the normal strain under constant stress loading is calculated (a) as well as the effective Young's modulus (b) and the shear modulus (c) for isotropically cracked medium</figcaption></figure>

405

406

407

408

409

410

411

412

413

414

415

416

417

418

419

420

421

422

423

424

425

426

427

428

429

430

431
For the remaining analyses performed with crack density parameters equal to 0.1, 0.15 and 0.2, the results show that there is a better agreement between the analytical results and the numerical ones than for the transverse isotropic case (Figure 15). Therefore, it can be concluded that the validity domain of the analytical approach for this particular case of isotropic behavior is extended to include more significantly damaged cases.

432

433

434

435

436

<figure><img src="image_17.png" /><figcaption>**Figure 15.** Influence of increasing the crack density on the normal strain determined from a creep test with $\epsilon$ equal to 0.1 (a), 0.15 (b) and 0.2 (c). Numerical results are then compared to the analytical ones for the tested cases</figcaption></figure>

437

438

439

Additional investigations were carried out with crack density parameters of 0.3 and 0.4 (Figure 16) for confirming the conclusion of the previous analysis. As shown in Figure 16, the results were consistently similar, maintaining a relatively good agreement between the analytical and numerical approaches even for higher crack density parameter. Consequently, the validity of the analytical solution can be considered as reasonable up to crack density parameter reaching 0.3-0.4. This is an important conclusion since the reliable analytical solution can be used with confidence to perform structure simulations, thus benefiting from a significant reduction in computational time. The application of this analytical approach within the framework of investigating the behavior of damaged claystone can be advantageous for a simplified isotropic representation. As mentioned, this applies primarily to situations where cracks are assumed to have a random orientation in some areas of the damaged zone.

440

441

442

443

444

445

446

447

448

449

450

451

<figure><img src="image_18.png" /><figcaption>**Figure 16.** Results of a uniaxial creep test in an isotropic medium considering crack densities of 0.2, 0.3 and 0.4</figcaption></figure>

452

453

454

455
## 5. Discussion

This research suggests an analysis of micro-cracked viscoelastic materials by comparing a numerical approach based on the finite element method and an analytical simplified approach. This analytical approach consists in approximating the effective moduli of the material in the Laplace Carson space by Burgers models, thereby providing analytical solutions that can be directly compared with the numerical calculations. One key advantage is that the parameters of these approximated Burgers models can be calculated as a function of the parameters of the linear viscoelastic models defining the behavior of the sound material, and the crack density parameter. This makes the approach adapted for structure simulations as the models can be implemented as a classical behavior law. The presented results include uniaxial creep tests, involving the calculation of normal strain as well as the evaluation of homogenized moduli affected by the cracking. The analysis for the transverse isotropic medium shows that the analytical approach only seems satisfactory for crack density parameters below around 0.1-0.15. For higher values, the numerical calculations diverge significantly from the analytical results. It is envisaged to investigate the use of more adapted analytical approaches for this case of parallel cracks such as methods deriving from the differential scheme, as it has previously proven to be valid in a linear elastic context [35].

In the case of an isotropic crack distribution, the analytical results were consistently closer to the numerical ones, including in the case of a crack density parameter up to 0.4. In that case, although the material is considered severely damaged, the analytical calculations were still in line with the numerical results, proving the validity of this method in the specific case of randomly oriented cracks in an isotropic medium.

Numerical approaches can be used to increase crack density parameter to values greater than 1, which opens up possibilities for investigating phenomenon occurring in severely damaged materials. As crack density increases, interactions between cracks inevitably intensify, creating microstructures highly dense in fractures. This can lead to the formation of preferential cracking paths, particularly in the case of randomly oriented cracks, generating a percolation phenomenon particularly interesting to explore. It would also be pertinent to consider irregular crack shapes to achieve an accurate representation of cracks in reality. Numerical approaches provide means of studying more complex crack geometries, which can be compared with the conventional example of penny shaped cracks. Overall, this study highlights the importance of combining analytical and numerical approaches to fully understand the behavior of cracked materials, taking into account the respective limitations of each method and exploiting the results to obtain a more accurate and representative model.

## 6. Conclusions

In this study, the combined application of numerical and analytical approaches was explored in order to characterize the behavior of the cracked viscoelastic COx claystone. This approach features a concrete application for representing the argillite in the context of the Cigéo project devoted to the storage of nuclear waste. Indeed, cracking of the host rock results from the excavation of the tunnels at 500m below ground level. Given the significant creep properties of the claystone, it is essential to account for both cracking and creep of the material. For simplification purposes, a linear viscoelastic behavior has been adopted, which allows the application of the Laplace Carson transform. This reduces the problem to a linear elastic behaviour in the Carson space, greatly simplifying the equations to provide straightforward analytical solutions. The Burgers model is selected as the viscoelastic model as it is capable, with a minimum number of parameters, of providing an appropriate representation of the behavior of the undamaged COx claystone.

In addition, two different cases were investigated, one involving an isotropic distribution of cracks in the solid matrix, resulting in an isotropic behavior, and the second featuring a parallel distribution of cracks, resulting in a transversely isotropic behavior of
the material. A simplified approach was used to calculate the homogenized moduli, assuming that the macroscopic behavior of the damaged material conforms to a Burgers model. This approach makes the identification of the model parameters feasible as a function of the parameters of the sound material and the crack density parameter, using simplified analytical expressions. While examining the transversely isotropic case, increasing the crack density parameter induces a progressive divergence in the approaches, and the analytical and numerical solutions are no longer fully comparable for parameters greater than 0.1-0.15. Advanced numerical studies can be pursued by investigating the material's behavior at crack densities greater than 0.2. The objective is to represent the behaviour of a severely damaged material, particularly in the context of describing the damaged zone behaviour for applications associated with the Cigéo project. With this background, two different cases of crack orientations are examined, consistent with in situ observations of the damaged zone. The first case involves a parallel crack distribution with a preferential orientation, while the second features randomly oriented cracks in an isotropic matrix.

Ultimately, further testing with an isotropic medium showed that the analytical solution gives fairly good results even for relatively high crack density parameters, as shown when comparing it with the numerical results. Potentially, time-consuming calculations of complex microstructures using the finite element method can be replaced by analytical solutions, as they have been proven to provide very similar results in certain conditions. Furthermore, this approach would provide basic solutions for multi-scale applications, in order to obtain the macroscopic response of a cracked material. Further investigations would involve examining different crack geometries, given that the crack shape can influence the calculations of the effective properties of the material. Alternatively, different analytical approaches can be explored, particularly in the case of aligned cracks exhibiting a transverse isotropic behavior. This could involve the use of homogenization schemes that could produce results more consistent with the numerical ones. The application of the differential scheme would be relevant, given its proven validity in the case of a transversely isotropic medium in a linear elastic framework [19]. The application of this approach with the differential scheme in a linear viscoelastic context is currently being investigated, and will be discussed in a forthcoming contribution.

**Funding:** Special recognition is expressed to our external partners at Andra and EDF for their financial support, which was essential to the realization of this project.

**Conflicts of Interest:** The authors declare no conflicts of interest in this work.

## References

1. Andra, "Dossier d'options de sûreté – Partie exploitation (DOS-EXPL)." 2016.

2. R. De La Vaissière, G. Armand, and J. Talandier, "Gas and water flow in an excavation-induced fracture network around an underground drift: A case study for a radioactive waste repository in clay rock," *Journal of Hydrology*, vol. 521, pp. 141–56, Feb. 2015, doi: 10.1016/j.jhydrol.2014.11.067.

3. G. Armand, F. Leveau, C. Nussbaum, R. de La Vaissiere, A. Noiret, D. Jaeggi, P. Landrein, and C. Righini, "Geometry and Properties of the Excavation-Induced Fractures at the Meuse/Haute-Marne URL Drifts," *Rock Mechanics and Rock Engineering*, vol. 47, pp. 21–41, Jan. 2014, doi: 10.1007/s00603-012-0339-6.

4. S. Rahal, A. Sellier, and G. Casaux-Ginestet, "Poromechanical consolidation and basic creep interactions around tunnel excavation," *International Journal of Rock Mechanics and Mining Sciences*, vol. 94, pp. 55–63, Apr. 2017.

5. Andra, "Dossier d'options de sûreté – Partie après fermeture (DOS-AF)." 2016.

6. S. T. Nguyen, L. Dormieux, Y. L. Pape, and J. Sanahuja, "A Burger Model for the Effective Behavior of a Microcracked Viscoelastic Solid," *International Journal of Damage Mechanics*, vol. 20, no. 8, pp. 1116–29, Nov. 2011, doi: 10.1177/1056789510395554.
7. J. Bluthé, B. Bary, and E. Lemarchand, "Contribution of FE and FFT-based methods to the determination of the effective elastic and conduction properties of composite media with flat inclusions and infinite contrast," *International Journal of Solids and Structures*, vol. 216, Jan. 2021, doi: 10.1016/j.ijsolstr.2021.01.007.

8. J. J. Marigo, C. Maurini, and K. Pham, "An overview of the modelling of fracture by gradient damage models," *Meccanica*, vol. 51, no. 12, pp. 3107–28, Dec. 2016, doi: 10.1007/s11012-016-0538-4.

9. Y. J. Cao, W. Q. Shen, J. F. Shao, and W. Wang, "A novel FFT-based phase field model for damage and cracking behavior of heterogeneous materials," *International Journal of Plasticity*, vol. 133, p. 102786, Oct. 2020, doi: 10.1016/j.ijplas.2020.102786.

10. L. Dormieux and D. Kondo, "Stress-based estimates and bounds of effective elastic properties: The case of cracked media with unilateral effects," *Computational Materials Science*, vol. 46, no. 1, pp. 173–9, Jul. 2009, doi: 10.1016/j.commatsci.2009.02.027.

11. A. Tsitova, F. Bernachy-Barbe, B. Bary, and F. Hild, "Experimental and numerical analyses of the interaction of creep with mesoscale damage in cementitious materials," *Mechanics of Materials*, vol. 184, p. 104715, Sep. 2023, doi: 10.1016/j.mechmat.2023.104715.

12. J. Mazars, "A description of micro- and macroscale damage of concrete structures," *Engineering Fracture Mechanics*, vol. 25, no. 5, pp. 729–37, Jan. 1986, doi: 10.1016/0013-7944(86)90036-6.

13. L. Dormieux, D. Kondo, and F. J. Ulm, "A micromechanical analysis of damage propagation in fluid-saturated cracked media," *Comptes Rendus. Mécanique*, vol. 334, no. 7, pp. 440–6, 2006, doi: 10.1016/j.crme.2006.05.007.

14. P. P. Castañeda and J. R. Willis, "The effect of spatial distribution on the effective behavior of composite materials and cracked media," *Journal of the Mechanics and Physics of Solids*, vol. 43, no. 12, pp. 1919–51, Dec. 1995, doi: 10.1016/0022-5096(95)00058-Q.

15. A. V. Hershey, "The Elasticity of an Isotropic Aggregate of Anisotropic Cubic Crystals," *Journal of Applied Mechanics*, vol. 21, no. 3, pp. 236–40, Jun. 2021, doi: 10.1115/1.4010899.

16. Z. Hashin, "The differential scheme and its application to cracked materials," *Journal of the Mechanics and Physics of Solids*, vol. 36, no. 6, pp. 719–34, Jan. 1988, doi: 10.1016/0022-5096(88)90005-1.

17. B. Budiansky and R. J. O'connell, "Elastic moduli of a cracked solid," *International Journal of Solids and Structures*, vol. 12, no. 2, pp. 81–97, Jan. 1976.

18. J. F. Shao and J. W. Rudnicki, "A microcrack-based continuous damage model for brittle geomaterials," *Mechanics of Materials*, vol. 32, no. 10, pp. 607–19, Oct. 2000, doi: 10.1016/S0167-6636(00)00024-7.

19. J. Bluthé, B. Bary, and E. Lemarchand, "Closure of parallel cracks: Micromechanical estimates versus finite element computations," *European Journal of Mechanics - A/Solids*, vol. 81, p. 103952, May 2020.

20. J. M. Ricaud and R. Masson, "Effective properties of linear viscoelastic heterogeneous media: Internal variables formulation and extension to ageing behaviours," *International Journal of Solids and Structures*, vol. 46, no. 7, pp. 1599–606, Apr. 2009, doi: 10.1016/j.ijsolstr.2008.12.007.

21. S. Beurthey and A. Zaoui, "Structural morphology and relaxation spectra of viscoelastic heterogeneous materials," *European Journal of Mechanics - A/Solids*, vol. 19, no. 1, pp. 1–16, Jan. 2000, doi: 10.1016/S0997-7538(00)00157-1.

22. R. Schapery, "A Simple Collocation Method for Fitting Viscoelastic Models to Experimental Data," Feb. 1962.

23. L. P. Ferreira, P. C. F. Lopes, R. Leiderman, F. T. S. Aragão, and A. M. B. Pereira, "An image-based numerical homogenization strategy for the characterization of viscoelastic composites," *International Journal of Solids and Structures*, vol. 267, p. 112142, Apr. 2023, doi: 10.1016/j.ijsolstr.2023.112142.

24. J. D. Eshelby, "The Determination of the Elastic Field of an Ellipsoidal Inclusion, and Related Problems," *Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences*, vol. 241, no. 1226, pp. 376–96, 1957.

25. L. Dormieux and D. Kondo, "Approche micromécanique du couplage perméabilité-endommagement," *Comptes Rendus Mécanique*, vol. 332, no. 2, Art. no. 2, Feb. 2004, doi: 10.1016/j.crme.2003.11.003.
26. S. T. Nguyen, "Propagation de fissures et endommagement par microfissures des matériaux viscoélastiques linéaires non vieillissants," phdthesis, Université Paris-Est, 2010. [Online]. Available: https://pastel.hal.science/tel-00598511 597

27. S. Cuvilliez, I. Djouadi, S. Raude, and R. Fernandes, "An elastoviscoplastic constitutive model for geomaterials: Application to hydromechanical modelling of claystone response to drift excavation," *Computers and Geotechnics*, vol. 85, pp. 321–40, May 2017. 599

28. M. Mánica, A. Gens, J. Vaunat, and D. F. Ruiz, "A time-dependent anisotropic model for argillaceous rocks. Application to an underground excavation in Callovo-Oxfordian claystone," *Computers and Geotechnics*, vol. 85, pp. 341–50, May 2017. 601

29. C. Zhang and T. Rothfuchs, "Experimental study of the hydro-mechanical behaviour of the Callovo-Oxfordian argillite," *Applied Clay Science*, vol. 26, no. 1, Art. no. 1, Aug. 2004, doi: 10.1016/j.clay.2003.12.025. 603

30. G. Armand, N. Conil, J. Talandier, and D. M. Seyedi, "Fundamental aspects of the hydromechanical behaviour of Callovo-Oxfordian claystone: From experimental studies to model calibration and validation," *Computers and Geotechnics*, vol. 85, pp. 277–86, May 2017. 605

31. S. T. Nguyen, "Generalized Kelvin model for micro-cracked viscoelastic materials," *Engineering Fracture Mechanics*, vol. 127, pp. 226–34, Sep. 2014, doi: 10.1016/j.engfracmech.2014.06.010. 608

32. S. T. Nguyen, M. Vu, N. Vu, N. Thi Thu Nga, A. Gasser, and A. Rekik, "Generalized Maxwell model for micro-cracked viscoelastic materials," *International Journal of Damage Mechanics*, vol. 26, Apr. 2015, doi: 10.1177/1056789515608231. 610

33. N. G. March, D. R. Gunasegaram, and A. B. Murphy, "Evaluation of computational homogenization methods for the prediction of mechanical properties of additively manufactured metal parts," *Additive Manufacturing*, vol. 64, p. 103415, Feb. 2023, doi: 10.1016/j.addma.2023.103415. 612

34. L. J. Walpole, "Elastic Behavior of Composite Materials: Theoretical Foundations," in *Advances in Applied Mechanics*, vol. 21, C. S. Yih, Ed., Elsevier, 1981, pp. 169–242. doi: 10.1016/S0065-2156(08)70332-6. 615

35. J. Bluthe, "Modélisation couplée et simulation numérique des phénomènes de colmatage des argilites," phdthesis, Paris Est, 2019. [Online]. Available: http://www.theses.fr/2019PESC1027 617

36. T. Helfer, B. Michel, J. M. Proix, M. Salvo, J. Sercombe, and M. Casella, "Introducing the open-source mfront code generator: Application to mechanical behaviours and material knowledge management within the PLEIADES fuel element modelling platform," *Computers & Mathematics with Applications*, vol. 70, no. 5, p. 994, 2015, doi: 10.1016/j.camwa.2015.06.027. 619

37. B. Bary, C. Bourcier, and T. Helfer, "Analytical and 3D numerical analysis of the thermoviscoelastic behavior of concrete-like materials including interfaces," *Advances in Engineering Software*, vol. 112, pp. 16–30, Oct. 2017, doi: 10.1016/j.advengsoft.2017.06.006. 622

38. T. Honório, B. Bary, and F. Benboudjema, "Multiscale estimation of ageing viscoelastic properties of cement-based materials: A combined analytical and numerical approach to estimate the behaviour at early age," *Cement and Concrete Research*, vol. 85, pp. 137–55, 2016, doi: 10.1016/j.cemconres.2016.03.010. 625

39. L. Charpin and A. Ehrlacher, "Estimating the poroelastic properties of cracked materials," *Acta Mech*, vol. 225, no. 9, pp. 2501–19, Sep. 2014, doi: 10.1007/s00707-013-1082-0. 628

40. K. Vasylevskyi, B. Drach, and I. Tsukrov, "On micromechanical modeling of orthotropic solids with parallel cracks," *International Journal of Solids and Structures*, vol. 144–145, pp. 46–58, Jul. 2018, doi: 10.1016/j.ijsolstr.2018.02.038. 630

41. J. H. Hannay and J. F. Nye, "Fibonacci numerical integration on a sphere," *J. Phys. A: Math. Gen.*, vol. 37, no. 48, p. 11591, Nov. 2004, doi: 10.1088/0305-4470/37/48/005. 632