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

: 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 behaviour, 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 disc shape. Intersections between randomly oriented cracks are accounted for, thereby highlighting a potential percolation phenomenon. The effects of micro-cracks on the material’s behaviour 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.


Introduction
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][2][3].Characterizing the behaviour 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 behaviour 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 performed using the finite element method with the Cast3M (http://www-cast3m.cea.fr/(version 2022, accessed on 4 April 2022)) software.This comparison between the two approaches offers key validation requirements when it comes to modelling the macroscopic behaviour 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 crack orientations as a function of distance from the excavated drifts.In the vicinity of the galleries, crack 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][8][9][10][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's macroscopic properties [14][15][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 modelling [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) behaviour, 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 behaviour in the Laplace-Carson space.This provides a straightforward framework to have recourse to for 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 behaviour of the COx claystone is represented here by non-ageing 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 behaviour, 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 behaviour [22,23].
Most classical approaches used to estimate the behaviour 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 behaviour 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 behaviour of the material can be represented by a set of Burgers models.This efficient method provides an analytical solution of the model 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 using 3D samples generated with disc-shaped voids of zero thickness, thus mimicking very accurately the behaviour of micro-cracked materials.In this context, the material's behaviour will be analysed by considering relevant cases of crack orientations including randomly oriented cracks as well as a parallel crack distribution in an isotropic medium.

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 behaviour of a rock mass must assume time-dependence.A rigorous approach is to use viscoelastic or viscoplastic behaviour laws to describe the evolution of the creep phenomenon for geomaterials [4,27,28].In analysing the behaviour 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-ageing 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 behaviour, while requiring a limited number of parameters.It is intended to use this model to represent the behaviour 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 two different Burgers models to describe independently the shear and bulk moduli denoted by µ and k, respectively.
Modelling 2024, 5, FOR PEER REVIEW 3 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 using 3D samples generated with disc-shaped voids of zero thickness, thus mimicking very accurately the behaviour of micro-cracked materials.In this context, the material's behaviour will be analysed by considering relevant cases of crack orientations including randomly oriented cracks as well as a parallel crack distribution in an isotropic medium.

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 behaviour of a rock mass must assume time-dependence.A rigorous approach is to use viscoelastic or viscoplastic behaviour laws to describe the evolution of the creep phenomenon for geomaterials [4,27,28].In analysing the behaviour 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-ageing 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 behaviour, while requiring a limited number of parameters.It is intended to use this model to represent the behaviour 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 two different Burgers models to describe independently the shear and bulk moduli denoted by µ and , respectively.
Equation ( 1) expresses the effective compressibility and shear moduli in Laplace-Carson space, denoted by  * and µ * , respectively, where  is the LC variable.The eight parameters of the Burgers model can be distinguished in Equation (1).The Maxwell part is represented by four parameters, including spherical components ( , ղ ) and deviatoric components (µ , ղ ).Additionally, the Kelvin part of the model is also composed of four parameters, with ( , ղ ) describing the spherical parameters and (µ , ղ ) 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 .Additionally, the Kelvin part of the model is also composed of four parameters, with (k K , η s K ) describing the spherical parameters and (µ K , η d K ) 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 behaviour.These tests include different series of loadings and have been designed to examine the deferred behaviour 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 conformity, 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 (ε nn and ε tt , respectively) in the time space are identified by calculating the inverse of the Carson-Laplace transform of: 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 σ 0 can be deduced in LC space from Equation (2): imum deviatoric stress that the material can withstand.The argillite is considered aniso-tropic, and therefore, its properties change as a function of the imposed loading with respect to the bedding plane.Here, for conformity, 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 (ԑ and ԑ , respectively) in the time space are identified by calculating the inverse of the Carson-Laplace transform of: 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  can be deduced in LC space from Equation (2): Inversion of the normal strain ԑ * from LC space yields: This expression involves the eight parameters of the Burgers model and the time.Inversion of the normal strain ε * nn from LC space yields: This expression involves the eight parameters of the Burgers model and the time.
Overall, the proposed model provides an acceptable representation of the argillite's behaviour in both the 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 Table 1.

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 is very close.The interested reader is invited to refer to the cited papers for more details that would be too lengthy to include here.

Isotropic Crack Distribution
Considering an isotropic distribution of cracks in a solid matrix (Figure 3), where cracks are characterized by the crack density parameter ϵ 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 for all cracks (Equation ( 5)).

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 is very close.The interested reader is invited to refer to the cited papers for more details that would be too lengthy to include here.

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

𝜖 = 𝑁𝑎
(5) The apparent effective stiffness tensor in the Laplace-Carson space defined by ℂ * is therefore isotropic, expressed in the (,) basis in Equation ( 6).The tensors  and  are, respectively, the projectors of the spherical and deviatoric parts of the fourth-order unit tensor.The homogenized moduli  * and µ * are calculated as a function of the crack density parameter , taking into account two distinct loading conditions: spherical loading, which is used to determine the compressibility modulus, and deviatoric loading, which is used to determine the shear modulus.A homogenized stiffness matrix is then constructed in the LC space, reflecting the material cracking conditions [26].The apparent effective stiffness tensor in the Laplace-Carson space defined by C hom * is therefore isotropic, expressed in the (J,K) basis in Equation ( 6).The tensors J and K are, respectively, the projectors of the spherical and deviatoric parts of the fourth-order unit tensor.
The homogenized moduli k hom * and µ hom * are calculated as a function of the crack density parameter ϵ, taking into account two distinct loading conditions: spherical loading, which is used to determine the compressibility modulus, and deviatoric loading, which is used to determine the shear modulus.A homogenized stiffness matrix is then constructed in the LC space, reflecting the material cracking conditions [26].
To simplify, here, the homogenized material behaviour in the LC space is studied by considering a general analogy with a homogenized macroscopic Burgers model, i.e., we assume that the behaviour of the homogenized material can be approached by Burgers models.This makes it possible to identify the homogenized modules k hom * and µ hom * in the LC space using a simplified form: This maintains the same structure of the Burgers model formulated in Equation ( 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 (Equation ( 5)) is established, as illustrated in Figure 4.In summary, this approach provides an alternative method for characterizing the behaviour of a micro-cracked material in the Laplace-Carson space [6].
To simplify, here, the homogenized material behaviour in the LC space is studied by considering a general analogy with a homogenized macroscopic Burgers model, i.e., we assume that the behaviour of the homogenized material can be approached by Burgers models.This makes it possible to identify the homogenized modules  * and µ * in the LC space using a simplified form: This maintains the same structure of the Burgers model formulated in Equation ( 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 (Equation ( 5)) is established, as illustrated in Figure 4.In summary, this approach provides an alternative method for characterizing the behaviour of a micro-cracked material in the Laplace-Carson space [6].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 (Equations ( 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].

Parallel Crack Distribution
In this section, the effect of a random distribution of parallel cracks in a linear viscoelastic material, as described in the previous sections, will be discussed.In this case, a transversely isotropic behaviour is obtained (Figure 5).In this context, it is convenient to work with the expression of the material compliance tensor in LC space  * () (Equation ( 11)), which is the inverse of the stiffness tensor using the Walpole basis [34].Using the 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 (Equations ( 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].

Parallel Crack Distribution
In this section, the effect of a random distribution of parallel cracks in a linear viscoelastic material, as described in the previous sections, will be discussed.In this case, a transversely isotropic behaviour 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) (Equation ( 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.
By decomposing the compliance tensor in the Walpole basis (Equation ( 11)), we obtain six scalars with indices  ranging from 1 to 6 denoted  * ().The tensors  are defined in accordance with the specifications given in the Walpole basis [34].Equation (12) describes the expressions of the scalars  * () in the Laplace-Carson space.Note that the directions (  ,   ,   ), respectively, correspond to the conventional directions in a Cartesian coordinate system (x, y, z). is the direction normal to the fractures along   (Figure 5).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 scalar indices in the Walpole's basis, the expression for the moduli  * and  * as a function of the cracking density yields: provided that  * and  * 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  * is obtained by imposing a uniaxial tensile load along the normal to the cracks, while that of  * requires a shear loading in the form  (  ⊗  +  ⊗   ), e.g., with  = 1 MPa.As the other modules will not be affected by the cracks, it is By decomposing the compliance tensor in the Walpole basis (Equation ( 11)), we obtain six scalars with indices 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].Equation ( 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, correspond 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).
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 scalar indices in the Walpole's basis, the expression for the moduli s hom * 2 and s hom * 4 as a function of the cracking density yields: provided that s hom * 2 and s hom * 4 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 hom * 2 is obtained by imposing a uniaxial tensile load along the normal to the cracks, while that of s hom * 4 requires a shear loading in the form σ (e 2 ⊗ n + n ⊗ e 2 ), e.g., with σ = 1 MPa.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.
The same identification technique as for the isotropic case was used, considering that the homogenized behaviour 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 eight coefficients in Equations ( 13) and ( 14) as a function of crack density (Equation ( 5)), as shown in Figure 6.The moduli affected by cracking in a Cartesian coordinate system, such as E hom 1 and µ hom 12 , can be determined from calculations in the Walpole basis as a function of the coefficients s hom 2 and s hom 4 [34].As mentioned above, investigating the transversely isotropic case is representative of the in situ behaviour 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.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 eight coefficients in Equations ( 13) and ( 14) as a function of crack density (Equation ( 5)), as shown in Figure 6.The moduli affected by cracking in a Cartesian coordinate system, such as  and , can be determined from calculations in the Walpole basis as a function of the coefficients and  [34].As mentioned above, investigating the transversely isotropic case is representative of the in situ behaviour 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.

Three-Dimensional Microstructures and Meshes
The Burgers model is implemented in the Mfront (https://tfel.sourceforge.net/gallery.html(accessed on 2 February 2022)) library, offering several advantages and convenient access for Cast3M users.Mfront is a code generator that enables convenient implementation of behavioural 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  = 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).

Numerical Simulations 4.1. Three-Dimensional Microstructures and Meshes
The Burgers model is implemented in the Mfront (https://tfel.sourceforge.net/gallery.html (accessed on 2 February 2022)) library, offering several advantages and convenient access for Cast3M users.Mfront is a code generator that enables convenient implementation of behavioural 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 Σ 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).The 3D sample generation procedure is based on a random distribution of inclusions of prescribed shape and size in a box [37,38].Disc-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 automatic software connected to Salome (https://www.salome-platform.org(accessed on 8 February 2022)) (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 The 3D sample generation procedure is based on a random distribution of inclusions of prescribed shape and size in a box [37,38].Disc-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 automatic software connected to Salome (https://www.salome-platform.org(accessed on 8 February 2022)) (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 crack mesh.Here, as the cracks are void, no interface elements are introduced, allowing us to drastically reduce the calculation time at the REV scale [7,40].
The 3D sample generation procedure is based on a random distribution of inclusions of prescribed shape and size in a box [37,38].Disc-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 automatic software connected to Salome (https://www.salome-platform.org(accessed on 8 February 2022)) (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 crack mesh.Here, as the cracks are void, no interface elements are introduced, allowing us to drastically reduce the calculation time at the REV scale [7,40].A series of test cases were examined featuring an increasing range of crack density parameters, and comparisons were subsequently conducted with analytical results.Start- 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 four 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 behaviour 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.
Modelling 2024, 5, FOR PEER REVIEW 10 ing with relatively low crack density parameters, the four 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 behaviour 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.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.
The analysis of Figure 10 demonstrates that the microstructures generated are fairly representative, and with minimized deviations, based on the calculation of the effective 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.
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 behaviour.The loading applied is a stress of σ = 1 MPa with Σ = σ e 1 ⊗ e 1 , constant with time.Simulations last around 2 h 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 are 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 number 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 (Tables 2 and 3), with periodic boundary conditions and the material subjected to a constant stress of σ = 1 MPa.Two simulations are required to calculate the effective properties: one with Σ = σ e 1 ⊗ e 1 for the calculation of the homogenized Young's modulus, and the second one with Σ = σ (e 1 ⊗ e 2 + e 2 ⊗ e 1 ) to determine the homogenized shear modulus.For all subsequent simulations, a mesh size of around 1.5 million elements was used.It is 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.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 (KUBCs), uniform stress boundary conditions (SUBCs) and finally periodic boundary conditions (PBCs).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 in 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 assumes that they provide the most accurate results, and they will be used as a reference for comparison with analytical results.

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  = 1 MPa.The time step 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

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 Σ 11 = 1 MPa.The time step 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 µ 12 = µ 13 , considering the case of parallel cracks where the Ox axis is the direction perpendicular to the corresponding fracture 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 Σ 11 = 1 MPa to calculate the effective Young's modulus, and the second requiring a deviatoric loading at Σ 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 transform, 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 (Equations ( 13) and ( 14)), which is subsequently used to determine the Young's modulus and the homogenized shear modulus.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 a 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.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 13c).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 param- 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 a 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.
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 13c).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.
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 a 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.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 13c).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.

Isotropic Case Results
We analyse 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 crack orientations.The first test is based on an REV with a defined

Isotropic Case Results
We analyse 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 crack orientations.The first test is based on an 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 (Equations ( 9) and ( 10)), enabling direct calculation of the homogenized Young's modulus and shear modulus in real time space.The normal strain ε nn is calculated as a function of the Young's modulus E hom and the applied constant stress σ = 1 MPa in LC space (Equation ( 3)).
Modelling 2024, 5, FOR PEER REVIEW 14 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 (Equations ( 9) and ( 10)), enabling direct calculation of the homogenized Young's modulus and shear modulus in real time space.The normal strain ԑ is calculated as a function of the Young's modulus  and the applied constant stress  = 1 MPa in LC space (Equation (3)).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 behaviour is extended to include more significantly damaged cases.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 behaviour is extended to include more significantly damaged cases.
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 parameters.Consequently, the validity of the analytical solution can be considered as reasonable up to a 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 behaviour 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.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 behaviour is extended to include more significantly damaged cases.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 parameters.Consequently, the validity of the analytical solution can be considered as reasonable up to a 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 an- alytical approach within the framework of investigating the behaviour 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.

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 behaviour 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 behaviour 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].

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 behaviour 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 behaviour 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 the crack density parameter to values greater than 1, which opens up possibilities for investigating phenomena 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 that is 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 a 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 behaviour of cracked materials, taking into account the respective limitations of each method and exploiting the results to obtain a more accurate and representative model.

Conclusions
In this study, the combined application of numerical and analytical approaches was explored in order to characterize the behaviour 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 behaviour 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 behaviour 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 behaviour, and the second featuring a parallel distribution of cracks, resulting in a transversely isotropic behaviour of the material.A simplified approach was used to calculate the homogenized moduli, assuming that the macroscopic behaviour 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 behaviour 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 behaviour.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.

Figure 1 .
Figure 1.Schematic representation of the Burgers model.

Figure 1 .
Figure 1.Schematic representation of the Burgers model.Equation (1) expresses the effective compressibility and shear moduli in Laplace-Carson space, denoted by k * and µ * , respectively, where p is the LC variable.The eight parameters of the Burgers model can be distinguished in Equation (1).The Maxwell part is represented by four parameters, including spherical components (k M , η s M ) and deviatoric components (µ M , η d M ).Additionally, the Kelvin part of the model is also composed of four parameters, with (k K , η s K ) describing the spherical parameters and (µ K , η d K ) 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 behaviour.These tests include different series of loadings and have been designed to examine the deferred behaviour of the material for different

Figure 2 .
Figure 2. Identification of the Burgers model parameters based on experimental creep test results.Figure 2. Identification of the Burgers model parameters based on experimental creep test results.

Figure 2 .
Figure 2. Identification of the Burgers model parameters based on experimental creep test results.Figure 2. Identification of the Burgers model parameters based on experimental creep test results.

Figure 3 .
Figure 3. Schematic representation of randomly oriented crack distribution in a solid matrix.

Figure 3 .
Figure 3. Schematic representation of randomly oriented crack distribution in a solid matrix.

Figure 4 .
Figure 4. Effective moduli calculation as a function of crack density, considering an isotropic crack distribution in the solid matrix.

Figure 4 .
Figure 4. Effective moduli calculation as a function of crack density, considering an isotropic crack distribution in the solid matrix.

Figure 5 .
Figure 5. Parallel crack distribution in the solid matrix.

Figure 5 .
Figure 5. Parallel crack distribution in the solid matrix.

Figure 6 .
Figure 6.Determination of the moduli affected by cracking for the transversely isotropic case in the Walpole basis as a function of crack density.

Figure 6 .
Figure 6.Determination of the moduli affected by cracking for the transversely isotropic case in the Walpole basis as a function of crack density.

Modelling 2024, 5 , 9 Figure 7 .
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.

Figure 7 .
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.

Figure 8 .
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.

Figure 8 .
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.

Figure 9 .
Figure 9.Even distribution of points across the sphere in order to define crack orientations in the microstructure and ensure an isotropic behaviour.

Figure 9 .
Figure 9.Even distribution of points across the sphere in order to define crack orientations in the microstructure and ensure an isotropic behaviour.

Modelling 2024, 5 , 11 Figure 10 .
Figure 10.Evaluation of REV representativity by varying the number of cracks distributed in the microstructure while keeping a constant crack density ϵ = 0.2.The effective Young's modulus in the direction normal to the cracking plane for the transverse isotropic case is calculated (left) as well as the effective Young's modulus for the isotropic case (right), and periodic stress boundary conditions were applied.

Figure 10 .
Figure 10.Evaluation of REV representativity by varying the number of cracks distributed in the microstructure while keeping a constant crack density ϵ = 0.2.The effective Young's modulus in the direction normal to the cracking plane for the transverse isotropic case is calculated (left) as well as the effective Young's modulus for the isotropic case (right), and periodic stress boundary conditions were applied.

Modelling 2024, 5 ,Figure 11 .
Figure 11.Effects of the applied boundary conditions KUBCs, SUBCs and PBCs 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.

Figure 11 .
Figure 11.Effects of the applied boundary conditions KUBCs, SUBCs and PBCs 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.

Modelling 2024, 5 ,Figure 12 .
Figure 12.Comparison between numerical and analytical results with  = 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 the transverse isotropic medium (b,c).

Figure 12 .
Figure 12.Comparison between numerical and analytical results with ϵ = 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 the transverse isotropic medium (b,c).

Figure 14 .
Figure 14.Comparison between numerical and analytical results with  = 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 an isotropically cracked medium.

Figure 15 .
Figure 15.Influence of increasing the crack density on the normal strain determined from a creep

Figure 14 .
Figure 14.Comparison between numerical and analytical results with ϵ = 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 an isotropically cracked medium.

Figure 14 .
Figure 14.Comparison between numerical and analytical results with  = 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 an isotropically cracked medium.

Figure 15 .
Figure 15.Influence of increasing the crack density on the normal strain determined from a creep test with  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.

Figure 15 .
Figure 15.Influence of increasing the crack density on the normal strain determined from a creep test with ϵ 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.

Figure 16 .
Figure 16.Results of a uniaxial creep test in an isotropic medium considering crack densities of 0.2, 0.3 and 0.4.

Figure 16 .
Figure 16.Results of a uniaxial creep test in an isotropic medium considering crack densities of 0.2, 0.3 and 0.4.

Table 1 .
Identified viscoelastic parameters of the COx claystone.

Table 2 .
Mesh convergence study considering an isotropic behaviour: computation of the homogenized Young's modulus and shear modulus.The relative error with respect to the last column is given in parentheses.

Table 3 .
Mesh convergence study considering a transversely isotropic behaviour: computation of the effective out-of-plane Young's modulus  and shear modulus  .The relative error with respect to the reference case (last column) is calculated.

Table 2 .
Mesh convergence study considering an isotropic behaviour: computation of the homogenized Young's modulus and shear modulus.The relative error with respect to the last column is given in parentheses.

Table 3 .
Mesh convergence study considering a transversely isotropic behaviour: computation of the effective out-of-plane Young's modulus E hom 1 and shear modulus µ hom 12 .The relative error with respect to the reference case (last column) is calculated.