Numerical implementation of a micromechanics-based viscoelastic model for geomaterials with isotropically distributed micro-fractures

The constitutive behavior of geomaterials is generally affected by the presence at different scales of discontinuity surfaces with different sizes and orientations. According to their mechanical behavior, such discontinuities can be distinguished as cracks or fractures. Fractures are interfaces that can transfer normal and tangential stresses, whereas cracks are discontinuities without stress transfer. Regarding the formulation of the behavior of materials with isotropic distribution of micro-cracks or fractures, previous works had essentially focused on their instantaneous response induced by structural loading. Few research works have addressed time-dependent (delayed) behavior of such materials. The present contribution describes the formulation and computational implementation of a micromechanics-based modeling for viscoelastic media with an isotropic distribution of micro-fractures. The homogenized viscoelastic properties are assessed by implementing a reasoning based on linear homogenization schemes (Mori-Tanaka) together with the correspondence principle for non-aging viscoelastic materials. It is shown that the homogenized viscoelastic behavior can be described by means of a generalized Maxwell rheological model. The computational implementation is developed within the finite element framework to analyze the delayed behavior of geomaterials with the presence of isotropically distributed micro-fractures under plane strain conditions. Several examples of applications are presented with the aim to illustrate the performance of the finite element modeling. The assessment of the approach accuracy and the corresponding code verification are performed by comparing the numerical predictions with analytical solutions for simple and complex geo-structures.

dependent respose addressing both fracture mechanics and viscoelasticity have been conducted recently. In order to study and characterize that behavior, several mathematical models have been analytically formulated. Among these models, there are those based on a micromechanics framework.
Some important micromechanics-based works that considered delayed behavior were developed by Nguyen et al. (2010Nguyen et al. ( , 2013 and Nguyen and Dormieux. (2014). These researchers formulated a micromechanics-based model for a viscoelastic medium where the heterogeneities are considered as cracks. However, as mentioned above, there are no stress transfer in this type of discontinuities. Subsequently, Nguyen´s analysis was extended by Aguiar and Maghous (2018), who added the fracture behavior studied by Maghous et al. (2011).
In this context, the main objective of the present paper is related to the finite element modeling and implementation within a computational code of the viscoelastic constitutive behavior originally formulated in Aguiar and Maghous (2018) for the fractured medium. It is emphasized that the main novelty of the approach lies in the ability of such a computational tool to allow for effective and accurate analysis of the delayed response of structures built in fractured media.
The paper content is organized as follows. Section 2 presents succinctly the theoretical formulation for the homogenized viscoelastic behavior of fractured mediums, considering an isotropic distribution of micro-fractures. Section 3 describes the algorithm used for time discretization of the viscoelastic model. Section 4 aims to verify the implementation code by presenting several examples of application where the accuracy of the numerical responses is assessed through comparison with the corresponding analytical solutions. Finally, section 5 illustrates the predictive capabilities of the computational tool by analyzing the delayed response of a complex geo-structure.

FORMULATION OF THE HOMOGENIZED VISCOELASTIC MODEL
As mentioned above, Aguiar and Maghous (2018) formulated the homogenized viscoelastic properties for a fractured medium. This approach is based on prior research works developed within the framework of elastic fractured mediums (Maghous et al. 2011(Maghous et al. , 2014 and the elastic-viscoelastic correspondence principle in the context of non-aging materials. There are two principal steps in order to formulate the macroscopic viscoelastic properties: in first place, it is necessary to formulate the overall elastic properties of the medium taking advantage of linear homogenization schemes. The second step is to associate the elastic-viscoelastic correspondence principle (Le, 2008) with a specific procedure to determinate the inverse Laplace-Carson transform of a function.
Since the paper focus is the study of the computational implementation of the viscoelastic model, only the formulation of the homogenized viscoelastic properties in a medium with isotropically distributed fractures will be shown below in a concise way. The formulation of the elastic properties in a fracture medium, the complete procedure of inverse Laplace-Carson transform and other basic concepts of the theoretical model can be found in Aguiar and Maghous (2018).

Homogenized viscoelastic properties in a medium with isotropic distribution of micro-fractures
Let  denote the representative elementary volume (REV) of a homogeneous material (the matrix phase) cut by a discrete distribution of randomly oriented short fractures (the embedded fractures phase). The REV is chosen so as to be statistically representative of the fractured medium and to comply with scale separation conditions (Zaoui, 2002). In particular, the characteristic size of fractures is assumed to be small with respect to the dimension of the REV, which in turn should be small when compared to a characteristic dimension of the structure body. For sake of simplicity and without loss of generality, we shall assume that the REV is a Lipschitz domain and its boundary  is located in the matrix phase (i.e.,  does not intersect the fractures phase). As classically considered in micromechanics-based approaches, the loading of the REV can be defined by homogeneous strain type boundary conditions (Zaoui, 2002): Latin American Journal of Solids and Structures, 2019, 16(9), e238 3/33 concentration problem stated on the REV under applied loading  (Aguiar and Maghous, 2018). In such conditions, the homogenized constitutive law for a fractured viscoelastic medium can be defined as where  is the macroscopic stress,  is the strain concentration tensor associated with the viscoelastic problem defined on the REV. The fourth order tensor  represents the link between the local strain at any point of the REV and applied macroscopic strain . In the above relationship, the symbol  denotes the Boltzmann operator and the symbol   represents the volume average over the matrix. It is pointed out that the time derivative ( ) should be understood in the sense of distributions in order to account for possible strain jumps in time. For practical determination of tensor  hom  , it is convenient to resort to the elastic-viscoelastic correspondence principle. This principle enables to formulate the time-domain viscoelastic problem in terms of an equivalent elastic problem in the Laplace-Carson domain.
In view of all the above mentioned, it is possible to take advantage of the elastic properties of a fractured medium. It should be taken into account that in the elastic problem, fractures are considered as randomly oriented oblate ellipsoids. This particular feature can be considered as an isotropic distribution of fractures in the REV as shown in Figure 1. It is also important to point out that for the determination of the homogenized stiffness tensor  hom in the elastic problem, it was adopted the Mori-Tanaka linear homogenization scheme. The same scheme will be used for determination of the relaxation tensor in the Laplace-Carson space, which can be denoted by  hom  .
The inverse Laplace transform of ( ) F p can be named as hom ( ) R t and the corresponding expression is where the subscripts k and  indicate that the parameters i R and i D are referring to the relaxation behavior in bulk and relaxation behavior in shear, respectively. The fourth-order tensors  and  are the spherical and deviatoric tensors, respectively. These tensors can be defined as  1 1 1 3   and      , where  is the fourth-order unitary tensor. Comparing Eq. (4) with the relaxation function of a Generalized Maxwell rheological model in Figure 2, it is noticed that exist a mathematical equivalence between them. From this equivalence, it can be stated that the mentioned rheological model represents the delayed behavior of non-aging linear viscoelastic materials and, therefore, it can be considered as the exact viscoelastic model. It should be emphasized that such a mechanical analogy results from direct application of the homogenization procedure (Aguiar and Maghous, 2018), without any a priori simplifying assumption that would be introduced to lead this analogy.
Latin American Journal of Solids and Structures, 2019, 16(9), e238 5/33 Figure 2 Generalized Maxwell model and its relaxation function (excerpt from Aguiar and Maghous, 2018) In the case of macroscopic isotropy, it is deduced that the bulk and shear viscoelastic behavior in (5) can also be exactly represented by a Generalized Maxwell rheological model, as displayed in Figure 3.  Aguiar and Maghous, 2018) As seen in Eq. (5), for an isotropic distribution of randomly oriented fractures, tensor  hom ( ) t depends on the moduli hom k and hom  . In this way, taking advantage of the elastic-viscoelastic correspondence principle, viscoelastic moduli hom k  and hom   in Laplace-Carson domain must be obtained in order to apply subsequently the inverse Laplace-Carson transform. The expressions for these viscoelastic moduli in Laplace-Carson domain are shown in (6).
where parameters ( , ) depend on the rheological model adopted for the viscoelastic behavior of the constituents (matrix and fractures). As mentioned above, having defined the moduli hom k  and hom   , it only remains to apply the inverse Laplace-Carson transform procedure in order to obtain the moduli in time domain: hom ( ) k t and hom ( ) t  .

ALGORITHM FOR TIME DISCRETIZATION OF THE VISCOELASTIC MODEL
In the previous section, it was defined the theoretical formulation of the homogenized viscoelastic model for geomaterials with isotropically distributed micro-fractures. Now, it is necessary to develop a numerical programming scheme which considers a time discretization in the finite element framework. Essentially, what needs to be implemented is the viscoelastic constitutive law that describes the material behavior when applied a homogeneous strain boundary condition.
Based on the conditions mentioned above, subroutines to be generated will be responsible for the calculation of homogenized stresses at Gauss points in finite element mesh of the analyzed structure. It is important to emphasize that the finite element approach is performed at the macroscopic scale in which the fractured medium is replaced by the homogenized viscoelastic material determined from previous micromechanical reasoning. In other words, the finite element discretization refers to the structure regarded as a homogenized structure. In this respect, the Greek letters used for denote stress and strain macroscopic tensors will be  and  , respectively.
The relaxation tensor  hom ( ) t for an isotropic distribution of microfractures, as well as the corresponding viscoelastic moduli hom ( ) k t and hom ( ) t  were defined in (5). It must be emphasized that numerical expressions for the viscoelastic moduli are previously obtained by any calculus software. Thus, these moduli are considered as input data for the program. In this way, the relaxation tensor can be expressed in matrix form as Once the relaxation tensor  hom ( ) t is computed, the macroscopic stress tensor can be determined, according to the constitutive law, as follows where ) (0  is the instantaneous strain tensor that appears right in the moment of load application. Latin American Journal of Solids and Structures, 2019, 16(9), e238 7/33 Regarding the code development, it will be considered n time steps with the same size t  and 0 0 t  . For the 1 i  time step, the corresponding time can be defined as An incremental procedure is applied in the integral equation in (9). Thus, the stress tensor at time 1 A linear variation is assumed for the strain in the interval where the variable of integration This assumption allows to express the derivative In order to update the stress tensor equation, (12) is replaced in (11). In a summarized form, the stress tensor where T K  and ( ) h t  are fourth-order tensors defined as At this point, it must be remembered that the relaxation function for a Generalized Maxwell rheological model have mathematical equivalence with the relaxation tensor  hom ( ) t components and it can be expressed in a generic way as where R  stands for the stiffness of spring model. R  and   are, respectively, the elastic spring stiffness and dashpot viscosity of branch  and p is the total number of branches. Since tensors T K  and 1 k A   depend on tensor  hom , their components can be conveniently expressed as Latin American Journal of Solids and Structures, 2019, 16(9), e238 8/33 The same equalities shown in (17) Similarly, for tensor where  corresponds to the bulk relaxation modulus and  to the shear relaxation modulus.
Latin American Journal of Solids and Structures, 2019, 16(9), e238 9/33 The precedent expressions make possible the stress calculation at Gauss point for all the considered step times. According to the way the stress tensor was formulated in (13), it is proposed to define the stress tensor components for three different time steps: TS = 1, TS = 2 and TS ≥ 3.
For TS = 1, where 1 1 i t t   and 0 i = , the stress tensor TS  can be defined as: : where and  1 is the homogenized relaxation tensor evaluated at TS = 1. It is known that 0  is the instantaneous strain tensor. The matrix form of tensor 1  is found in (25).
The matrix form of tensor TS  for TS ≥ 3 is presented in (29).  TS  TS   TS  TS   TS   T  T  TS  TS   T  T  TS  TS   T   TS   TS   TS   TS   TS   TS   T where SAE11, SAE22, SAE33 and SAE12 are defines as

VERIFICATION OF THE IMPLEMENTATION CODE
This section is focused on the verification of the implementation code through the analysis of two examples of application under plane strain conditions. Each one of them presents a time-dependent analytical solution that will be compared with the corresponding numerical responses. The correspondence between the delayed responses will serve to verify the developed code.
The same specimen will be used for the two examples. With respect to the specimen dimensions, it is a square with a size of 1.00 m. As regards the discretized geometry model, it is composed of four square finite elements, even though it is possible to use just one finite element since stress and strain fields are homogeneous for all these examples. Figure 4 displays the finite elements, nodes and dimensions of the geometric model.  Tables 1 and 2 show, respectively, the mechanical parameters used for matrix and fractures. These parameters were reported in Le (2008) for viscoelastic moduli of concrete, in Bart (2000) for elastic fracture behavior and the viscosity parameters of fractures were arbitrarily chosen according to Aguiar and Maghous (2018). The fracture density parameter is 0.20   and the number of fractures per unit volume    k . It is emphasized that the above rheological models used for the viscoelastic behaviors of the constituents (matrix, fractures) refer to the microscopic scale, while the homogenized viscoelastic behavior is still given by (4), which reflects the behavior at macroscopic scale. In particular, the parameters associated at this scale with the generalized Maxwell model (spring stiffness i E and dashpot viscosity i  ) are calculated in each case from the viscoelastic properties of the matrix and fractures.
Comparison between analytical solutions and numerical responses for the two examples of application under plane strain conditions is presented in the following subsections.

Example 1: Plane strain compression test
where '( ) k t stands for viscoelastic bulk modulus and '( ) t  for viscoelastic shear modulus. Variables y and z refer to the coordinates of position vector X . These moduli are defined as  (8) and (9) are showing the comparison between the analytical solution defined in (31) and the numerical prediction for vertical strain zz  and horizontal strain yy  .  (10) and (11).  The results observed from the above figures show a good agreement between the analytical and finite element solutions, exhibiting a relative error that remains lower than 6.5% for the time interval of analysis. Consequently, the implementation code has been successfully verified.
In the case of the specimen subjected to a vertical imposed displacement, the analytical solution tensorial fields are As shown in Figure 17, the analytical and numerical predictions show a good agreement for the early stages of analysis with increasing discrepancy with time due to the very small values of xx  , which leads to irrelevant relative error evaluations. Once again, the implementation code has been successfully verified.

Example 2: Compression under prescribed displacement rate
It will be considered an imposed instantaneous displacement 0.001 m.  Similarly, comparisons between analytical solution defined in (36) and numerical response for Burger-Maxwell case is illustrated in Figures (21) and (22). A direct error analysis is illustrated in Figure (23) by depicting the relative error between the analytical and numerical predictions obtained for the stress components along time.

Figure 23 Relative Error between the numerical and the analytical stress predictions
The results shown in the preceding figures indicate a satisfactory agreement with a maximum relative error less than 3% along the time interval of analysis. Consequently, the implementation code has been successfully verified for all the examples presented in this section.

ANALYSIS OF DELAYED BEHAVIOR OF DEEP UNDERGROUND GALLERIES
The implementation code was verified in section 4 through several examples using a simple quadrilateral specimen. In addition to that verification, this section is focused on the analysis of geo-structures of greater complexity under plane strain conditions. There exist some geo-structures that allow this type of analysis, such as the deep underground galleries or tunnels. The study of these geo-structures enables to determine a cross-section of interest for carrying out a two-dimensional analysis.
It will be considered two types of cross-section for deep tunnels: circular and horseshoe. Concerning to circular cross-section, two situations will be presented: unlined cross-section and lined cross-section. In both situations related to the circular tunnel, additional verifications of the implementation code performance will be performed by comparing the analytical solution and numerical response for delayed behavior in this type of geo-structures. Subsequently, for a horseshoe cross-section, its delayed response will be analyzed in order to highlight the predictive capabilities of this computational tool when it is not possible to have an analytical solution for viscoelastic behavior.

Unlined circular cross-section tunnel
The hypotheses considered in this problem are described as follows. The tunnel has radius R and a depth H, where H  R. The massif material is homogeneous and isotropic. The excavation process will be regarded as instantaneous.
The initial stress (geostatic) 0 1 p    , due to self-weight of the massif, will be uniform in the region around the gallery. The variable p is defined as p H   , where  represents the specific weight of the massif and 1 is the second-order unitary tensor. Due to symmetry determined by the tunnel geometry, only a quarter of the total structure will be used for respective analyses. The discretized geometry model consists of a 2D square domain of side length L H  . Figure 24 displays dimensions for the numerical analysis, as well as displacement and stress boundary conditions. It will be considered a depth H = 125.00 m , radius R = 2.50 m , specific weight 3 24000 N m   for calculations. The lateral and vertical pressure exerted by massif is 0.003 GPa p  . Concerning to viscoelastic characteristics of the material, the cases of delayed behavior adopted by matrix and fractures are the same as those used in section 4: Maxwell-Maxwell and Burger-Maxwell. Consequently, mechanical parameters correspond to those specified in Tables 1 and 2, as well as the fracture density parameter  and the number of fractures per unit volume  . Regarding to the discretized geometric model, Figure 25 shows the finite element mesh used for the analyses. It is composed of 192 finite element and 221 nodes. In relation to Maxwell-Maxwell case, the analytical solution, defined in (38) and (40), will be compared with the numerical response for radial convergence ( ) U t , radial displacement ( , ) rr r t  and radial stress ( , ) rr r t  . Regarding to radial displacements and stresses, the analysis will be performed for three different times: t1 corresponds to the instantaneous elastic response, t2 is an intermediate time and t3 corresponds to a very long time. Similarly, below is shown the comparison of analytical solution, defined in (38) and (40), with the numerical response for Burger-Maxwell case. As shown in Figures 26-31, analytical and numerical curves present a satisfactory correspondence. It can be considered as an additional verification of correct performance of the implementation code. Despite the consistent correspondence between analytical and numerical results, boundary effects related to displacements have been observed in Figures 27 and 30. Since boundary displacements are supposed to be zero in the discretized geometric model, it is necessary to perform future researches focused on discretized model dimensions and finite element mesh refinement in order to reduce these boundary effects.

Lined circular cross-section tunnel
The hypotheses considered in this section will be the same as those used in section 5.1. In this problem, there will be a perfectly rigid lining in the tunnel cross-section. This lining characteristic enables to obtain a particular analytical solution. Tunnel cross-section surface presents an instantaneous radial displacement 0 ( , ) u R t  after excavation that will be constrained with lining. Due to this displacement constraint, massif will exert a pressure ( ) q t on tunnel lining. Figure 32 shows the displacement conditions in tunnel cross-section surface after excavation and the pressure on the lining.  As shown in Figures 33-38, analytical and numerical curves present a satisfactory correspondence. This has been also an additional verification of correct performance of the implementation code. As in section 5.1, despite the consistent correspondence between analytical and numerical results, boundary effects related to displacements have been observed in Figures 34 and 37. Therefore, it is necessary to perform future researches focused on discretized model dimensions and finite element mesh refinement in order to reduce these boundary effects.

Horseshoe cross-section tunnel
This section is focused on showing the predictive capabilities or potentialities of the computational implementation when studied complex geo-structures that don't present an analytical viscoelastic solution due to their particular geometric characteristics. In this context, it was chosen a horseshoe cross-section deep tunnel to analyze its corresponding time-dependent displacement and stress fields.
Information concerning cross-section geometry, tunnel depth and specific weight of massif were obtained from Couto (2011) doctoral thesis, who performed some researches on an underground gallery located in Jagran, Pakistan. Regarding to the material viscoelastic properties, the cases of delayed behavior adopted by matrix and fractures are the same as those used in sections 5.1 and 5.2. In the same way, it will be used the fracture density parameter  , number of fractures per unit volume  and mechanical parameters from Tables 1 and 2 shown in section 4. Figure  39 As can be seen in Figure 39, there are two marked points A and B, whose corresponding axes (Y-axis and X-axis) will be used for analysis of displacements   Similarly, below is shown time-dependent displacements and stresses for Burger-Maxwell case. As shown in Figures 41-48, stresses and displacements present a coherent delayed behavior in the analyzed axes. As position vector moves away from cross-section surface, displacements values decrease and stresses values tend to the geostatic pressure p . In Figures 41, 43, 45 and 47, while boundary displacements tend to zero, residual values are still significant in the discretized geometric model. It is necessary to perform future researches focused on discretized model dimensions and finite element mesh refinement in order to reduce these boundary effects. Despite of these effects, the implementation code works in a satisfactory way for problems that address geo-structures with presence of isotropically distributed micro-fractures.

CONCLUSIONS
The study of delayed behavior in geomaterials becomes vitally important when geo-structures are subjected to long-term loading. This behavior has strong impact in geo-structures strains and becomes relevant when discontinuity surfaces appear in the matrix material. The homogenized viscoelastic model formulated by Aguiar and Maghous (2018) aims to represent, in an adequate way, this type of behavior. This constitutive law formulation was the theoretical basis for development of the implementation code in this paper. Based on the mentioned basis and the results produced by the computational tool, the following conclusions and suggesitions for future works were made: • Even though this model can use different rheological models for description of matrix and fractures delayed behavior, the general case adopts a Burger rheological model for both material constituents. The choice of rheological models for representation of delayed behavior depends on many factors, such as, for example, proposals studied in previous research works. It also depends on the possibility of using parameters obtained from laboratory tests and on the research objectives when studied viscosity effects on matrix and fractures.

•
The rheological model that represents, in an accurate way, the overall viscoelastic behavior of a homogenized medium with isotropic distribution of fractures is the Generalized Maxwell model. Due to its characteristics, this model is denominated as the exact model and was used for the present implementation. An isotropic distribution of micro-fractures enables to consider the matrix material as isotropic and, thus, the relaxation tensor can be represented by a combination of bulk and shear viscoelastic moduli.

•
The correct performance of the implementation code was verified through the examples of application shown in section 4 and through the analyses of deep circular cross-section tunnels in sections 5.1 and 5.2. This verification showed that numerical predictions have an optimal correspondence with the analytical solutions at short, intermediate and long times. It was observed that the greater the rheological model used for delayed behavior description, the greater the computational effort. Thus, Burger-Maxwell case required greater computational analysis time than Maxwell-Maxwell case.
• In the case of deep horseshoe cross-section tunnel analysis, it was noticed a coherent delayed behavior for both displacements and stresses. The greater the distance from cross-section surface, the less the displacement given by the numerical model. As position vector moves away from cross-section surface, stresses values tend to the geostatic pressure p related to the initial stress state. In this way, the predictive capabilities of the computational tool were shown in this paper.
• Many possibilities exist to extend this computational implementation work. For example, the study of mesh dimensions and refinement for deep tunnels with viscoelastic materials in order to reduce boundary effects related to displacements. It also can be studied the relation between geometric model size and the real tunnel depth with the aim of improving numerical results. It is suggested to add another type of fracture distribution to the implementation code and to analyze the delayed response when used another fracture density parameter or number of fractures per unit volume.
• Finally, one the main suggestions for extension of the present work is the validation of the viscoelastic model, using the parameters identification as a first step. It would be interesting to perform tridimensional analyses with the aim to assess the predictive capabilities of the model at short, intermediate and long times by comparing numerical predictions with results obtained from laboratory tests of a specific geomaterial subjected to long-term loadings. An example of these laboratory tests is the triaxial compression test.