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

Abstract

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.

Keywords:
Fracture; micromechanics; viscoelasticity; finite element

1 INTRODUCTION

One of the most common causes of material degradation in civil engineering structures is the presence of discontinuity surfaces, which can appear with different sizes and orientations. The existence of these small thickness regions increases the risk of mechanical properties deterioration during the service life of structures because of a significative reduction of stiffness, strength, ductility and a serious increase of permeability. From a mechanical behavior perspective, there are two principal types of discontinuities: cracks and fractures. Essentially, cracks have no stress transfer between faces, whereas fractures can transfer normal and tangential stresses. Consequently, a medium can be considered as cracked or fractured depending on the characteristics of its discontinuities.

Geomaterials under long-term loading present two types of response: instantaneous and delayed. The instantaneous response of fractured geomaterials has been object of research over the years (see for instance Budiansky and O’Connell, 1976). However, less attention has been given to delayed response in fractured mediums. Given the importance of delayed components on stress/strain response of geomaterials, researches focused on time-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. (2010, 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.

2 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, 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).

2.1 Homogenized viscoelastic properties in a medium with isotropic distribution of micro-fractures

Let $\Omega$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 $\partial \Omega$ is located in the matrix phase (i.e., $\partial \Omega$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): $\underset{\_}{\xi }(\underset{\_}{x})=\underset{\_}{\underset{\_}{\in }}\cdot \underset{\_}{x}\forall \underset{\_}{x}\in \partial \Omega$, where $\underset{\_}{\underset{\_}{\in }}$ represents the macroscopic strain (loading parameter), $\underset{\_}{\xi }$ is the displacement filed and $\underset{\_}{x}$ is the position vector. Figure 1 provides a schematized representation of the REV and associated loading mode.

From a constitutive viewpoint, the matrix material is assumed to exhibit a non-aging linear viscoelastic behavior defined by the fourth-order relaxation tensor ${ℂ}^S$. The fractures are discontinuities that are able to transfer stress and can therefore be regarded as interfaces with a specific non-aging linear viscoelastic behavior.

The formulation of non-aging viscoelastic properties of a homogenized medium involves the determination of the fourth-order relaxation tensor $ℝ{\mathrm{hom}}^{}$ components. This determination relies upon the solution to a viscoelastic concentration problem stated on the REV under applied loading $\underset{\_}{\underset{\_}{\in }}$ (Aguiar and Maghous, 2018). In such conditions, the homogenized constitutive law for a fractured viscoelastic medium can be defined as

$\underset{\_}{\underset{\_}{\Sigma }}=ℝ{\mathrm{hom}}^{}\circ \underset{\_}{\underset{\_}{\in }}={\displaystyle {\int }_{-\infty }^{\infty }{ℝ}^{\mathrm{hom}}(t-\tau ):\underset{\_}{\underset{\_}{\dot{\in }}}(\tau )d\tau }withℝ{\mathrm{hom}}^{}=\langle ℂs^{}\circ \mathbb{A}\rangle$ (1)

where $\underset{\_}{\underset{\_}{\Sigma }}$ is the macroscopic stress, $\mathbb{A}$is the strain concentration tensor associated with the viscoelastic problem defined on the REV. The fourth order tensor $\mathbb{A}$represents the link between the local strain at any point of the REV and applied macroscopic strain $\underset{\_}{\underset{\_}{\in }}$. In the above relationship, the symbol $\circ$ denotes the Boltzmann operator and the symbol $\langle \rangle$ represents the volume average over the matrix. It is pointed out that the time derivative $\underset{\_}{\underset{\_}{\dot{\in }}}(\tau )=\partial \underset{\_}{\underset{\_}{\in }}/\partial \tau$ should be understood in the sense of distributions in order to account for possible strain jumps in time. For practical determination of tensor$ℝ{\mathrm{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 $ℂ{\mathrm{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 $ℝ{\text{hom}\ast }^{}$.

Figure 1
Representative elementary volume of the fractured medium and loading mode

After tensor $ℝ{\mathrm{hom}\ast }^{}$ determination in Laplace-Carson domain, it becomes essential to express the relaxation tensor components in the real time domain. This raises the need for a specific inverse Laplace-Carson transform procedure $ℝ{\mathrm{hom}}^{}={\mathcal{L}}_c^{-1}(ℝ{\mathrm{hom}\ast }^{})$, which was developed by Aguiar and Maghous (2018). This procedure was conceived to work with the tensor $ℝ{\mathrm{hom}\ast }^{}$ obtained from a Mori-Tanaka homogenization scheme when used classical rheological models for description of matrix and fractures behavior. What follows is a summary of the main considerations of the procedure. Any non-zero component of relaxation tensor $ℝ{\mathrm{hom}\ast }^{}$ is generically referred to as $R^{\mathrm{hom}\ast }$. This function can be expressed as

$R^{\mathrm{hom}\ast }(p)=\frac{A(p)}{B(p)}=\frac{{\displaystyle \sum _{i=0}^n{a}_i{p}^i}}{{\displaystyle \sum _{i=0}^n{b}_i{p}^i}}$ (2)

where $n$ is the common degree of polynomials $A(p)$ and $B(p)$. The variable $p$ belongs to the Laplace-Carson domain. Scalar $a_i$ and $b_i$ are real coefficients that depend on the rheological viscoelastic models adopted for matrix and fractures. By sake of simnplicity, the following steps deal with $F(p)$, which corresponds to function $R^{\mathrm{hom}\ast }(p)$ changed to the laplace domain:

$F(p)=\frac{R^{\mathrm{hom}\ast }(p)}{p}=\frac{a_0}{p{b}_0}+\frac{C(p)}{B(p)}withC(p)=\frac{1}{b_0}{\displaystyle \sum _{i=1}^n\left[(b_0{a}_i-a_0{b}_i)p^{i-1}\right]}$ (3)

The inverse Laplace transform of $F(p)$ can be named as $R^{\mathrm{hom}}(t)$ and the corresponding expression is

$R^{\mathrm{hom}}(t)=\left[\frac{a_0}{b_0}+{\displaystyle \sum _{i=1}^n{D}_i}{e}^{R_{i}t}\right]Y(t)with{D}_i={\left[\frac{C(p)}{\partial B(p)/\partial p}\right]}_{p=R_i}$ (4)

where $R^{\mathrm{hom}}(t)$ represents the non-zero component of tensor $ℝ{\mathrm{hom}}^{}(t)$, $Y(t)$ represent the Heaviside function at origin and the scalar $R_i$ is the kth root of the polynomial $B(p)$. In the context of macroscopic isotropy, the generic function $R^{\mathrm{hom}}(t)$ refers to either bulk $k^{\mathrm{hom}}(t)$ or shear ${\mu }^{\mathrm{hom}}(t)$ relaxation moduli. The definition of the relaxation tensor in the case of randomly oriented fractures within an isotropic viscoelastic matrix and its corresponding moduli are presented in (5).

$\begin{array}{l}ℝ{\mathrm{hom}}^{}(t)=3k^{\mathrm{hom}}(t)\mathbb{J}+2{\mu }^{\mathrm{hom}}(t)\mathbb{K}\\ \\ k^{\mathrm{hom}}(t)=\left[\frac{a_0^k}{b_0^k}+{\displaystyle \sum _{i=1}^n{D}_i^k}{e}^{R_i^{k}t}\right]Y(t),{\mu }^{\mathrm{hom}}(t)=\left[\frac{a_0^{\mu }}{b_0^{\mu }}+{\displaystyle \sum _{i=1}^n{D}_i^{\mu }}{e}^{R_i^{\mu }t}\right]Y(t)\end{array}$ (5)

where the subscripts $k$ and $\mu$ indicate that the parameters $R_i^{}$and $D_i$ are referring to the relaxation behavior in bulk and relaxation behavior in shear, respectively. The fourth-order tensors $\mathbb{J}$ and $\mathbb{K}$ are the spherical and deviatoric tensors, respectively. These tensors can be defined as $\mathbb{J}=\frac{1}{3}\underset{\_}{\underset{\_}{1}}\otimes \underset{\_}{\underset{\_}{1}}$ and $\mathbb{K}=\mathbb{I}-\mathbb{J}$, where $\mathbb{I}$ 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.

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.

Figure 3
Generalized Maxwell model for viscoelastic moduli $k^{\mathrm{hom}}$ and ${\mu }^{\mathrm{hom}}$(excerpt from Aguiar and Maghous, 2018)

As seen in Eq. (5), for an isotropic distribution of randomly oriented fractures, tensor $ℝ{\mathrm{hom}}^{}(t)$ depends on the moduli $k^{\mathrm{hom}}$ and ${\mu }^{\mathrm{hom}}$. In this way, taking advantage of the elastic-viscoelastic correspondence principle, viscoelastic moduli $k^{\mathrm{hom}\ast }$ and ${\mu }^{\mathrm{hom}\ast }$ 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).

$k^{\mathrm{hom}\ast }=\frac{k^{s\ast }}{1+\epsilon M_k},{\mu }^{\mathrm{hom}\ast }=\frac{{\mu }^{s\ast }}{1+\epsilon M_{\mu }}$ (6)

where $\epsilon =\mathcal{N}{a}^3$ is the fracture density parameter, $\mathcal{N}$ is the number of fractures per unit volume and $a$ is the radius of the oblate fracture. The dimensionless functions $M_k(k^{s\ast },{\mu }^{s\ast },a{k}_n^{\ast },a{k}_t^{\ast })$ and $M_{\mu }(k^{s\ast },{\mu }^{s\ast },a{k}_n^{\ast },a{k}_t^{\ast })$ are defined as

$\begin{array}{l}{M}_k=\frac{\frac{4}{3}\pi k^{s\ast }/{\mu }^{s\ast }}{\pi {\kappa }_1+a{k}_n^{\ast }/{\mu }^{s\ast }},M_{\mu }=\frac{16\pi {\kappa }_4}{15}\frac{6{\kappa }_2+4{\kappa }_3+9\pi {\kappa }_4(3{\kappa }_1+{\kappa }_4)}{(3\pi {\kappa }_1{\kappa }_4+{\kappa }_2)(4{\kappa }_3+9\pi {\kappa }_4({\kappa }_1+{\kappa }_4))}\\ \\ {\kappa }_1=\frac{3k^{s\ast }+{\mu }^{s\ast }}{3k^{s\ast }+4{\mu }^{s\ast }},{\kappa }_2=\frac{3k_n^{\ast }a}{3k^{s\ast }+4{\mu }^{s\ast }},{\kappa }_3=\frac{3k_t^{\ast }a}{3k^{s\ast }+4{\mu }^{s\ast }},{\kappa }_4=\frac{{\mu }^{s\ast }}{3k^{s\ast }+4{\mu }^{s\ast }}\end{array}$ (7)

where parameters $(k^{s\ast },{\mu }^{s\ast })$ and $(k_n^{\ast },k_t^{\ast })$ depend on the rheological model adopted for the viscoelastic behavior of the constituents (matrix and fractures). As mentioned above, having defined the moduli $k^{\mathrm{hom}\ast }$ and ${\mu }^{\mathrm{hom}\ast }$, it only remains to apply the inverse Laplace-Carson transform procedure in order to obtain the moduli in time domain: $k^{\mathrm{hom}}(t)$ and ${\mu }^{\mathrm{hom}}(t)$.

3 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 $\underset{\_}{\underset{\_}{\sigma }}$ and $\underset{\_}{\underset{\_}{\epsilon }}$, respectively.

The relaxation tensor $ℝ{\mathrm{hom}}^{}(t)$ for an isotropic distribution of microfractures, as well as the corresponding viscoelastic moduli $k^{\mathrm{hom}}(t)$ and ${\mu }^{\mathrm{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

${ℝ}^{\mathrm{hom}}=\left[\begin{array}{cccccc}\left(k^{\mathrm{hom}}+\frac{4}{3}{\mu }^{\mathrm{hom}}\right)& \left(k^{\mathrm{hom}}-\frac{2}{3}{\mu }^{\mathrm{hom}}\right)& \left(k^{\mathrm{hom}}-\frac{2}{3}{\mu }^{\mathrm{hom}}\right)& 0& 0& 0\\ \left(k^{\mathrm{hom}}-\frac{2}{3}{\mu }^{\mathrm{hom}}\right)& \left(k^{\mathrm{hom}}+\frac{4}{3}{\mu }^{\mathrm{hom}}\right)& \left(k^{\mathrm{hom}}-\frac{2}{3}{\mu }^{\mathrm{hom}}\right)& 0& 0& 0\\ \left(k^{\mathrm{hom}}-\frac{2}{3}{\mu }^{\mathrm{hom}}\right)& \left(k^{\mathrm{hom}}-\frac{2}{3}{\mu }^{\mathrm{hom}}\right)& \left(k^{\mathrm{hom}}+\frac{4}{3}{\mu }^{\mathrm{hom}}\right)& 0& 0& 0\\ 0& 0& 0& 2{\mu }^{\mathrm{hom}}& 0& 0\\ 0& 0& 0& 0& 2{\mu }^{\mathrm{hom}}& 0\\ 0& 0& 0& 0& 0& 2{\mu }^{\mathrm{hom}}\end{array}\right]$ (8)

Once the relaxation tensor $ℝ{\mathrm{hom}}^{}(t)$ is computed, the macroscopic stress tensor can be determined, according to the constitutive law, as follows

$\underset{\_}{\underset{\_}{\sigma }}(t)={ℝ}^{\mathrm{hom}}(t):\underset{\_}{\underset{\_}{\epsilon }}(0)+{\displaystyle {\int }_0^t{ℝ}^{\mathrm{hom}}(s,t):\frac{\partial \underset{\_}{\underset{\_}{\epsilon }}}{\partial s}}(s)ds$ (9)

where $\underset{\_}{\underset{\_}{\epsilon }}(0)$ is the instantaneous strain tensor that appears right in the moment of load application.

Regarding the code development, it will be considered $n$ time steps with the same size $\Delta t$ and $t_0=0$. For the $i+\text{1}$ time step, the corresponding time can be defined as

$t_{i+1}=t_i+\Delta t=(i+1)\cdot \Delta twithi=\mathrm{0,}\mathrm{...},n-1$ (10)

An incremental procedure is applied in the integral equation in (9). Thus, the stress tensor at time $t_{i+1}$, is

$\underset{\_}{\underset{\_}{\sigma }}(t_{i+1})={ℝ}^{\mathrm{hom}}(t_{i+1}):\underset{\_}{\underset{\_}{\epsilon }}(0)+{\displaystyle \sum _{k=0}^i\left[{\displaystyle {\int }_{t_k}^{t_{k+1}}{ℝ}^{\mathrm{hom}}(s,t_{i+1}):\frac{\partial \underset{\_}{\underset{\_}{\epsilon }}}{\partial s}}(s)ds\right]}$ (11)

A linear variation is assumed for the strain in the interval where the variable of integration $s\in [t_k,t_{k+1}]$. This assumption allows to express the derivative $\partial \underset{\_}{\underset{\_}{\epsilon }}/\partial s$ as

$\frac{\partial \underset{\_}{\underset{\_}{\epsilon }}}{\partial s}=\frac{\Delta {\underset{\_}{\underset{\_}{\epsilon }}}_{k+1}}{\Delta t}=\frac{\underset{\_}{\underset{\_}{\epsilon }}(t_{k+1})-\underset{\_}{\underset{\_}{\epsilon }}(t_k)}{\Delta t}$ (12)

In order to update the stress tensor equation, (12) is replaced in (11). In a summarized form, the stress tensor $\underset{\_}{\underset{\_}{\sigma }}(t_{i+1})$ can be expressed as

$\underset{\_}{\underset{\_}{\sigma }}(t_{i+1})={\underset{˜}{K}}^T:\Delta {\underset{\_}{\underset{\_}{\epsilon }}}_{i+1}+{\underset{\_}{\underset{\_}{\sigma }}}^h(t)$ (13)

where ${\underset{˜}{K}}^T$ and ${\underset{\_}{\underset{\_}{\sigma }}}^h(t)$ are fourth-order tensors defined as

$\begin{array}{l}{\underset{˜}{K}}^T=\frac{1}{\Delta t}\left({\displaystyle {\int }_{t_i}^{t_{i+1}}{ℝ}^{\mathrm{hom}}(s,t_{i+1})}ds\right),{\underset{\_}{\underset{\_}{\sigma }}}^h(t)={\underset{\_}{\underset{\_}{\sigma }}}^0(t_{i+1})+{\displaystyle \sum _{k=0}^{i-1}\left[{\underset{˜}{A}}^{k+1}:\Delta {\underset{\_}{\underset{\_}{\epsilon }}}_{k+1}\right]}\\ \\ {\underset{\_}{\underset{\_}{\sigma }}}^0(t_{i+1})={ℝ}^{\mathrm{hom}}(t_{i+1}):\underset{\_}{\underset{\_}{\epsilon }}(0),{\underset{˜}{A}}^{k+1}=\frac{1}{\Delta t}\left({\displaystyle {\int }_{t_k}^{t_{k+1}}{ℝ}^{\mathrm{hom}}(s,t_{i+1})}ds\right)\end{array}$ (14)

At this point, it must be remembered that the relaxation function for a Generalized Maxwell rheological model have mathematical equivalence with the relaxation tensor $ℝ{\mathrm{hom}}^{}(t)$ components and it can be expressed in a generic way as

$R_{ijkl}(s,t)=R_{\infty }+{\displaystyle \sum _{\alpha =1}^p{R}_{\alpha }}{e}^{-\left(\frac{t-s}{{\tau }_{\alpha }}\right)}with{\tau }_{\alpha }=\frac{{\eta }_{\alpha }}{R_{\alpha }}$ (15)

where $R_{\infty }$ stands for the stiffness of spring model. $R_{\alpha }$ and ${\eta }_{\alpha }$ are, respectively, the elastic spring stiffness and dashpot viscosity of branch $\alpha$ and $p$ is the total number of branches. Since tensors ${\underset{˜}{K}}^T$ and ${\underset{˜}{A}}^{k+1}$ depend on tensor ${ℝ}^{\mathrm{hom}}$, their components can be conveniently expressed as

$K_{ijkl}^T=R_{\infty }+{\displaystyle \sum _{\alpha =1}^p\left[\frac{R_{\alpha }}{\left(\raisebox{1ex}{$\Delta t$}\!\left/ \!\raisebox{-1ex}{${\tau }_{\alpha }$}\right.\right)}\left(1-e^{-\left(\raisebox{1ex}{$\Delta t$}\!\left/ \!\raisebox{-1ex}{${\tau }_{\alpha }$}\right.\right)}\right)\right]},A_{ijkl}^{k+1}=R_{\infty }+{\displaystyle \sum _{\alpha =1}^p\left[\frac{R_{\alpha }}{\left(\raisebox{1ex}{$\Delta t$}\!\left/ \!\raisebox{-1ex}{${\tau }_{\alpha }$}\right.\right)}\left(e^{-(i-k)\left(\frac{\Delta t}{{\tau }_{\alpha }}\right)}-e^{-(i-k+1)\left(\frac{\Delta t}{{\tau }_{\alpha }}\right)}\right)\right]}$ (16)

It is known that in an isotropic distribution of fractures, for the non-zero components of tensor ${ℝ}^{\mathrm{hom}}$, the following applies

$\begin{array}{l}{R}_{1111}^{\mathrm{hom}}=R_{2222}^{\mathrm{hom}}=R_{3333}^{\mathrm{hom}},R_{1212}^{\mathrm{hom}}=R_{1313}^{\mathrm{hom}}=R_{2323}^{\mathrm{hom}}\\ \\ R_{1122}^{\mathrm{hom}}=R_{1133}^{\mathrm{hom}}=R_{2211}^{\mathrm{hom}}=R_{2233}^{\mathrm{hom}}=R_{3311}^{\mathrm{hom}}=R_{3322}^{\mathrm{hom}}\end{array}$ (17)

The same equalities shown in (17) applies also to tensors ${\underset{˜}{K}}^T$ and ${\underset{˜}{A}}^{k+1}$. Now, considering (5), it can be determined the expressions for non-zero components of tensor ${\underset{˜}{K}}^T$

$K_{1111}^T=\frac{a_0^{\omega }}{b_0^{\omega }}+{\displaystyle \sum _{\alpha =1}^{n\omega }\left(\frac{D_{\alpha }^{\omega }}{\left(-R_{\alpha }^{\omega }\Delta t\right)}(1-e^{R_{\alpha }^{\omega }\Delta t})\right)}+\frac{4}{3}\left[\frac{a_0^{\mu }}{b_0^{\mu }}+{\displaystyle \sum _{\alpha =1}^{n\mu }\left(\frac{D_{\alpha }^{\mu }}{\left(-R_{\alpha }^{\mu }\Delta t\right)}(1-e^{R_{\alpha }^{\mu }\Delta t})\right)}\right]$ (18)

$K_{1122}^T=\frac{a_0^{\omega }}{b_0^{\omega }}+{\displaystyle \sum _{\alpha =1}^{n\omega }\left(\frac{D_{\alpha }^{\omega }}{\left(-R_{\alpha }^{\omega }\Delta t\right)}(1-e^{R_{\alpha }^{\omega }\Delta t})\right)}-\frac{2}{3}\left[\frac{a_0^{\mu }}{b_0^{\mu }}+{\displaystyle \sum _{\alpha =1}^{n\mu }\left(\frac{D_{\alpha }^{\mu }}{\left(-R_{\alpha }^{\mu }\Delta t\right)}(1-e^{R_{\alpha }^{\mu }\Delta t})\right)}\right]$ (19)

$K_{1111}^T=\frac{a_0^{\omega }}{b_0^{\omega }}+{\displaystyle \sum _{\alpha =1}^{n\omega }\left(\frac{D_{\alpha }^{\omega }}{\left(-R_{\alpha }^{\omega }\Delta t\right)}(1-e^{R_{\alpha }^{\omega }\Delta t})\right)}+\frac{4}{3}\left[\frac{a_0^{\mu }}{b_0^{\mu }}+{\displaystyle \sum _{\alpha =1}^{n\mu }\left(\frac{D_{\alpha }^{\mu }}{\left(-R_{\alpha }^{\mu }\Delta t\right)}(1-e^{R_{\alpha }^{\mu }\Delta t})\right)}\right]$ (20)

Similarly, for tensor ${\underset{˜}{A}}^{k+1}$, the corresponding non-zero components are

$\begin{array}{l}{A}_{1111}^{k+1}=\frac{a_0^{\omega }}{b_0^{\omega }}+{\displaystyle \sum _{\alpha =1}^{n\omega }\left(\frac{D_{\alpha }^{\omega }}{\left(-R_{\alpha }^{\omega }\Delta t\right)}(e^{(i-k)R_{\alpha }^{\omega }\Delta t}-e^{(i-k+1)R_{\alpha }^{\omega }\Delta t})\right)}+\\ \frac{4}{3}\left[\frac{a_0^{\mu }}{b_0^{\mu }}+{\displaystyle \sum _{\alpha =1}^{n\mu }\left(\frac{D_{\alpha }^{\mu }}{\left(-R_{\alpha }^{\mu }\Delta t\right)}(e^{(i-k)R_{\alpha }^{\mu }\Delta t}-e^{(i-k+1)R_{\alpha }^{\mu }\Delta t})\right)}\right]\end{array}$ (21)

$\begin{array}{l}{A}_{1122}^{k+1}=\frac{a_0^{\omega }}{b_0^{\omega }}+{\displaystyle \sum _{\alpha =1}^{n\omega }\left(\frac{D_{\alpha }^{\omega }}{\left(-R_{\alpha }^{\omega }\Delta t\right)}(e^{(i-k)R_{\alpha }^{\omega }\Delta t}-e^{(i-k+1)R_{\alpha }^{\omega }\Delta t})\right)}-\\ \frac{2}{3}\left[\frac{a_0^{\mu }}{b_0^{\mu }}+{\displaystyle \sum _{\alpha =1}^{n\mu }\left(\frac{D_{\alpha }^{\mu }}{\left(-R_{\alpha }^{\mu }\Delta t\right)}(e^{(i-k)R_{\alpha }^{\mu }\Delta t}-e^{(i-k+1)R_{\alpha }^{\mu }\Delta t})\right)}\right]\end{array}$ (22)

$A_{1212}^{k+1}=2\left[\frac{a_0^{\mu }}{b_0^{\mu }}+{\displaystyle \sum _{\alpha =1}^{n\mu }\left(\frac{D_{\alpha }^{\mu }}{\left(-R_{\alpha }^{\mu }\Delta t\right)}(e^{(i-k)R_{\alpha }^{\mu }\Delta t}-e^{(i-k+1)R_{\alpha }^{\mu }\Delta t})\right)}\right]$ (23)

where $\omega$ corresponds to the bulk relaxation modulus and $\mu$ to the shear relaxation modulus.

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 $t_{i+1}=t_1$ and $i=0$, the stress tensor ${\underset{\_}{\underset{\_}{\sigma }}}^{TS}$ can be defined as:

${\underset{\_}{\underset{\_}{\sigma }}}^1={\underset{˜}{K}}^T:\Delta {\underset{\_}{\underset{\_}{\epsilon }}}^1+{ℝ}^1:{\underset{\_}{\underset{\_}{\epsilon }}}^0$ (24)

where $\Delta {\underset{\_}{\underset{\_}{\epsilon }}}^1={\underset{\_}{\underset{\_}{\epsilon }}}^1-{\underset{\_}{\underset{\_}{\epsilon }}}^0$ and ${ℝ}^1$ is the homogenized relaxation tensor evaluated at TS = 1. It is known that ${\underset{\_}{\underset{\_}{\epsilon }}}^0$ is the instantaneous strain tensor. The matrix form of tensor ${\underset{\_}{\underset{\_}{\sigma }}}^1$ is found in (25).

$\left[\begin{array}{c}{\sigma }_{11}^1\\ {\sigma }_{22}^1\\ {\sigma }_{33}^1\\ \begin{array}{l}{\sigma }_{12}^1\\ {\sigma }_{13}^1\\ {\sigma }_{23}^1\end{array}\end{array}\right]=\left[\begin{array}{c}\begin{array}{l}{K}_{1111}^T\Delta {\epsilon }_{11}^1+K_{1122}^T\Delta {\epsilon }_{22}^1+R_{1111}^1{\epsilon }_{11}^0+R_{1122}^1{\epsilon }_{22}^0\\ K_{1111}^T\Delta {\epsilon }_{22}^1+K_{1122}^T\Delta {\epsilon }_{11}^1+R_{1111}^1{\epsilon }_{22}^0+R_{1122}^1{\epsilon }_{11}^0\\ K_{1122}^T\Delta {\epsilon }_{11}^1+K_{1122}^T\Delta {\epsilon }_{22}^1+R_{1122}^1{\epsilon }_{11}^0+R_{1122}^1{\epsilon }_{22}^0\end{array}\\ K_{1212}^T\Delta {\epsilon }_{12}^1+R_{1212}^1{\epsilon }_{12}^0\\ 0\\ 0\end{array}\right]$ (25)

For TS = 2, where $t_{i+1}=t_2$ and $i=1$, the stress tensor ${\underset{\_}{\underset{\_}{\sigma }}}^{TS}$ can be defined as

${\underset{\_}{\underset{\_}{\sigma }}}^2={\underset{˜}{K}}^T:\Delta {\underset{\_}{\underset{\_}{\epsilon }}}^2+{ℝ}^2:{\underset{\_}{\underset{\_}{\epsilon }}}^0+{\underset{˜}{A}}^1:\Delta {\underset{\_}{\underset{\_}{\epsilon }}}^1$ (26)

where $\Delta {\underset{\_}{\underset{\_}{\epsilon }}}^2={\underset{\_}{\underset{\_}{\epsilon }}}^2-{\underset{\_}{\underset{\_}{\epsilon }}}^1$ and ${ℝ}^2$ is the homogenized relaxation tensor evaluated at TS = 2. ${\underset{˜}{A}}^1$ is the tensor ${\underset{˜}{A}}^{k+1}$ evaluated at TS=1 ($k=0$). The matrix form of tensor ${\underset{\_}{\underset{\_}{\sigma }}}^2$ is found in (27).

$\left[\begin{array}{c}{\sigma }_{11}^2\\ {\sigma }_{22}^2\\ {\sigma }_{33}^2\\ \begin{array}{l}{\sigma }_{12}^2\\ {\sigma }_{13}^2\\ {\sigma }_{23}^2\end{array}\end{array}\right]=\left[\begin{array}{c}\begin{array}{l}{K}_{1111}^T\Delta {\epsilon }_{11}^2+K_{1122}^T\Delta {\epsilon }_{22}^2+R_{1111}^2{\epsilon }_{11}^0+R_{1122}^2{\epsilon }_{22}^0+A_{1111}^1\Delta {\epsilon }_{11}^1+A_{1122}^1\Delta {\epsilon }_{22}^1\\ K_{1111}^T\Delta {\epsilon }_{22}^2+K_{1122}^T\Delta {\epsilon }_{11}^2+R_{1111}^2{\epsilon }_{22}^0+R_{1122}^2{\epsilon }_{11}^0+A_{1111}^1\Delta {\epsilon }_{22}^1+A_{1122}^1\Delta {\epsilon }_{11}^1\\ K_{1122}^T\Delta {\epsilon }_{11}^2+K_{1122}^T\Delta {\epsilon }_{22}^2+R_{1122}^2{\epsilon }_{11}^0+R_{1122}^2{\epsilon }_{22}^0+A_{1122}^1\Delta {\epsilon }_{11}^1+A_{1122}^1\Delta {\epsilon }_{22}^1\end{array}\\ K_{1212}^T\Delta {\epsilon }_{12}^2+R_{1212}^2{\epsilon }_{12}^0+A_{1212}^1\Delta {\epsilon }_{12}^1\\ 0\\ 0\end{array}\right]$ (27)

For TS ≥ 3, where $t_{i+1}\ge t_3$, $i=TS-1$ and $k\in \left[\mathrm{0,}i-1\right]$, the stress tensor ${\underset{\_}{\underset{\_}{\sigma }}}^{TS}$ can be defined as:

${\underset{\_}{\underset{\_}{\sigma }}}^{TS}={\underset{˜}{K}}^T:\Delta {\underset{\_}{\underset{\_}{\epsilon }}}^{TS}+{ℝ}^{TS}:{\underset{\_}{\underset{\_}{\epsilon }}}^0+{\displaystyle \sum _{k=0}^{i-1}\left[{\underset{˜}{A}}^{k+1}:\Delta {\underset{\_}{\underset{\_}{\epsilon }}}^{k+1}\right]}$ (28)

The matrix form of tensor ${\underset{\_}{\underset{\_}{\sigma }}}^{TS}$ for TS ≥ 3 is presented in (29).

$\left[\begin{array}{c}{\sigma }_{11}^{TS}\\ {\sigma }_{22}^{TS}\\ {\sigma }_{33}^{TS}\\ \begin{array}{l}{\sigma }_{12}^{TS}\\ {\sigma }_{13}^{TS}\\ {\sigma }_{23}^{TS}\end{array}\end{array}\right]=\left[\begin{array}{c}\begin{array}{l}{K}_{1111}^T\Delta {\epsilon }_{11}^{TS}+K_{1122}^T\Delta {\epsilon }_{22}^{TS}+R_{1111}^{TS}{\epsilon }_{11}^0+R_{1122}^{TS}{\epsilon }_{22}^0+SAE11\\ K_{1111}^T\Delta {\epsilon }_{22}^{TS}+K_{1122}^T\Delta {\epsilon }_{11}^{TS}+R_{1111}^{TS}{\epsilon }_{22}^0+R_{1122}^{TS}{\epsilon }_{11}^0+SAE22\\ K_{1122}^T\Delta {\epsilon }_{11}^{TS}+K_{1122}^T\Delta {\epsilon }_{22}^{TS}+R_{1122}^{TS}{\epsilon }_{11}^0+R_{1122}^{TS}{\epsilon }_{22}^0+SAE33\end{array}\\ K_{1212}^T\Delta {\epsilon }_{12}^{TS}+R_{1212}^{TS}{\epsilon }_{12}^0+SAE12\\ 0\\ 0\end{array}\right]$ (29)

where SAE11, SAE22, SAE33 and SAE12 are defines as

$\begin{array}{l}SAE11={\displaystyle \sum _{k=0}^{i-1}(A_{1111}^{k+1}\Delta {\epsilon }_{11}^{k+1}+A_{1122}^{k+1}\Delta {\epsilon }_{22}^{k+1})},SAE22={\displaystyle \sum _{k=0}^{i-1}(A_{1111}^{k+1}\Delta {\epsilon }_{22}^{k+1}+A_{1122}^{k+1}\Delta {\epsilon }_{11}^{k+1})}\\ \\ SAE33={\displaystyle \sum _{k=0}^{i-1}(A_{1122}^{k+1}\Delta {\epsilon }_{11}^{k+1}+A_{1122}^{k+1}\Delta {\epsilon }_{22}^{k+1})},SAE12={\displaystyle \sum _{k=0}^{i-1}(A_{1212}^{k+1}\Delta {\epsilon }_{12}^{k+1})}\end{array}$ (30)

Finally, expressions (25), (27) and (29) can be considered as the basis for the implementation code.

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

Figure 4
Elements of the discretized geometric model

With regard to viscoelastic characteristics of the material, the related parameters depend on the rheological model chosen to describe matrix and fractures behaviors at microscopic scale. The general case considers a Burger rheological model for representation of the matrix and fractures delayed behavior. However, for the examples in this paper, it will be analyzed two particular cases of viscoelastic behavior. These two cases combine rheological models adopted by the material constituents. In the first case, matrix and fractures adopt a Maxwell rheological model, while in the second case, matrix adopt a Burger rheological model and fractures adopt a Maxwell rheological model.

Figures 5 and 6 display, respectively, the rheological model elements and their corresponding parameters for case 1 and 2. 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 $\epsilon =0.20$ and the number of fractures per unit volume $\mathcal{N}=1.0$.

Figure 5
Rheological model elements for case 1 (Maxwell-Maxwell)

Figure 6
Rheological model elements for case 2 (Burger-Maxwell)

Table 1
Matrix mechanical parameters for Burger model

MATRIX
Parameter	Value	Unit	Parameter	Value	Unit
$k_{e,M}^s$	24.42	GPa	$k_{e,K}^s$	39.27	GPa
$k_{v,M}^s$	7.33x108	GPa-s	$k_{v,K}^s$	5.07x107	GPa-s
${\mu }_{e,M}^s$	13.27	GPa	${\mu }_{e,K}^s$	14.07	GPa
${\mu }_{v,M}^s$	3.88x108	GPa-s	${\mu }_{v,K}^s$	1.27x107	GPa-s

Table 2
Fractures mechanical parameters for Maxwell model

FRACTURES
Parameter	Value	Unit
$k_{e,M}^n$	42.22	GPa / m
$k_{v,M}^n$	7x107	GPa-s / m
$k_{e,M}^t$	16.88	GPa / m
$k_{v,M}^t$	4x107	GPa-s / m

Burger model is composed by a Maxwell component (subscript M) connected in series with a Kelvin component (subscript K). With respect to matrix material, elastic parameters associated with a spring in bulk are denoted by $k_{e,M}^s$ and $k_{e,K}^s$, whereas that associated in shear are ${\mu }_{e,M}^s$ and ${\mu }_{e,K}^s$. Viscosity parameters associated to a dashpot in bulk are denoted by $k_{v,M}^s$ and $k_{v,K}^s$, whereas that associated in shear are ${\mu }_{v,M}^s$ and ${\mu }_{v,K}^s$. As regards the fractures behavior, the elastic parameter associated with a spring under normal stress is denoted by $k_{e,M}^n$, whereas that associated under tangential stress is $k_{e,M}^t$. The viscosity parameter associated to a dashpot under normal stress is denoted by $k_{v,M}^n$, whereas that associated under tangential stress is $k_{v,M}^t$. 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 $E_i$ and dashpot viscosity ${\eta }_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.

4.1 Example 1: Plane strain compression test

Figure 7 shows a specimen with height H subjected to two different loading cases under plane strain conditions. Figure 7A shows an imposed compressive stress $-\sigma \cdot Y(t)$, whereas Figure 7B shows a vertical imposed displacement $-{\delta }_{z}Y(t){\underset{\_}{e}}_z$, where $Y(t)$ denotes the Heaviside step function at origin. Material will be considered homogeneous, isotropic and linear viscoelastic. The initial stress state is ${\underset{\_}{\underset{\_}{\sigma }}}^0=0$.

Figure 7
YZ-plane for specimen subjected to compressive stress (a) and imposed displacement (b)

The problem solution in both situations involves the calculation of time-dependent stress and displacement fields $(\underset{\_}{\underset{\_}{\sigma }}(t),\underset{\_}{\xi }(\underset{\_}{X},t))$ by using basic concepts of continuum mechanics and viscoelasticity. Once the static and kinematic admissibility of the respective solution fields for the specimen subjected to compressive stress has been determined, the analytical solution tensorial fields are

$\begin{array}{l}{\underset{\_}{\xi }}_{}(\underset{\_}{X},t)=\left[\frac{-\sigma (2\mu \text{'}(t)-3k\text{'}(t))}{4\mu \text{'}(t)(3k\text{'}(t)+\mu \text{'}(t))}\right]y{\underset{\_}{e}}_y+\left[\frac{-\sigma (3k\text{'}(t)+4\mu \text{'}(t))}{4\mu \text{'}(t)(3k\text{'}(t)+\mu \text{'}(t))}\right]z{\underset{\_}{e}}_z\\ \\ \underset{\_}{\underset{\_}{\epsilon }}(t)=\left[\frac{-\sigma (2\mu \text{'}(t)-3k\text{'}(t))}{4\mu \text{'}(t)(3k\text{'}(t)+\mu \text{'}(t))}\right]{\underset{\_}{e}}_y\otimes {\underset{\_}{e}}_y+\left[\frac{-\sigma (3k\text{'}(t)+4\mu \text{'}(t))}{4\mu \text{'}(t)(3k\text{'}(t)+\mu \text{'}(t))}\right]{\underset{\_}{e}}_z\otimes {\underset{\_}{e}}_z\\ \\ \underset{\_}{\underset{\_}{\sigma }}(t)=\left[\frac{-\sigma (3k\text{'}(t)-2\mu \text{'}(t))}{2(3k\text{'}(t)+\mu \text{'}(t))}\right]{\underset{\_}{e}}_x\otimes {\underset{\_}{e}}_x-\sigma {\underset{\_}{e}}_z\otimes {\underset{\_}{e}}_z\end{array}$ (31)

where $k\text{'}(t)$ stands for viscoelastic bulk modulus and $\mu \text{'}(t)$ for viscoelastic shear modulus. Variables $y$ and $z$ refer to the coordinates of position vector $\underset{\_}{X}$. These moduli are defined as

$k\text{'}(t)=\frac{1}{{\mathcal{L}}_c^{-1}(1/k_{(p)}^{\mathrm{hom}\ast })},\mu \text{'}(t)=\frac{1}{{\mathcal{L}}_c^{-1}(1/{\mu }_{(p)}^{\mathrm{hom}\ast })}$ (32)

where $k_{(p)}^{\mathrm{hom}\ast }$ and ${\mu }_{(p)}^{\mathrm{hom}\ast }$ are the homogenized bulk and shear moduli in Laplace-Carson domain, respectively. ${\mathcal{L}}_c^{-1}$ denotes the inverse Laplace-Carson transform. These moduli depend on the rheological models adopted for representation of matrix and fractures behavior. It will be considered a compressive stress $\sigma =0.02GPa$ for calculations.

Referring to Maxwell-Maxwell case, Figures (8) and (9) are showing the comparison between the analytical solution defined in (31) and the numerical prediction for vertical strain ${\epsilon }_{zz}$ and horizontal strain ${\epsilon }_{yy}$.

Figure 8
Comparison of delayed responses for vertical strain ${\epsilon }_{zz}$ (Maxwell-Maxwell)

Figure 9
Comparison of delayed responses for horizontal strain ${\epsilon }_{yy}$ (Maxwell-Maxwell)

Similarly, the comparison between analytical solution defined in (31) and numerical response for Burger-Maxwell case are depicted in Figures (10) and (11).

Figure 10
Comparison of delayed responses for vertical strain ${\epsilon }_{zz}$ (Burger-Maxwell)

Figure 11
Comparison of delayed responses for horizontal strain ${\epsilon }_{yy}$ (Burger-Maxwell)

As regard the error analysis, Figure (12a) and Figure (12b) plot the relative error between the analytical and numerical predictions as function of time.

Figure 12
Relative Error between the numerical and the analytical strain predictions

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

$\begin{array}{l}\underset{\_}{\xi }(\underset{\_}{X},t)=\left[-\left(\frac{3k(t)-2\mu (t)}{3k(t)+4\mu (t)}\right){\epsilon }_{z}y\right]{\underset{\_}{e}}_y+{\epsilon }_{z}z{\underset{\_}{e}}_z\\ \\ \underset{\_}{\underset{\_}{\epsilon }}(t)=\left[-\left(\frac{3k(t)-2\mu (t)}{3k(t)+4\mu (t)}\right){\epsilon }_z\right]{\underset{\_}{e}}_y\otimes {\underset{\_}{e}}_y+{\epsilon }_z{\underset{\_}{e}}_z\otimes {\underset{\_}{e}}_z\\ \\ \underset{\_}{\underset{\_}{\sigma }}(t)=\left[\frac{2\mu (t)(3k(t)-2\mu (t))}{3k(t)+4\mu (t)}{\epsilon }_z\right]{\underset{\_}{e}}_x\otimes {\underset{\_}{e}}_x+\left[\frac{4\mu (t)(3k(t)+\mu (t))}{3k(t)+4\mu (t)}{\epsilon }_z\right]{\underset{\_}{e}}_z\otimes {\underset{\_}{e}}_z\end{array}$ (33)

where variables $y$ and $z$ are coordinates of the position vector $\underset{\_}{X}$. The strain caused by imposed displacement is denoted by ${\epsilon }_z=-{\delta }_z/H$. The bulk and shear moduli in time domain are denoted by $k(t)$ and $\mu (t)$, respectively. The numerical expressions of $k(t)$ and $\mu (t)$ can be obtained by means of calculus softwares and depend on the rheological models adopted for matrix and fractures. It will be considered an imposed displacement ${\delta }_z=0.01m.$ for calculations.

In the Maxwell-Maxwell case, Figures (13) and (14) present the comparison between analytical solution defined in (33) and numerical response for vertical stress ${\sigma }_{zz}$ and stress ${\sigma }_{xx}$ perpendicular to YZ-plane.

Figure 13
Comparison of delayed responses for stress ${\sigma }_{zz}$ (Maxwell-Maxwell)

Figure 14
Comparison of delayed responses for stress ${\sigma }_{xx}$ (Maxwell-Maxwell)

Similarly, comparisons between analytical solution defined in (33) and numerical response for Burger-Maxwell case are shown in Figures (15) and (16).

Figure 15
Comparison of delayed responses for stress ${\sigma }_{zz}$ (Burger-Maxwell)

Figure 16
Comparison of delayed responses for stress ${\sigma }_{xx}$ (Burger-Maxwell)

The variations in time of the relative error between the analytical and numerical predictions are plotted in Figures (17a) and (17b).

Figure 17
Relative Error between the numerical and the analytical stress predictions

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 ${\sigma }_{xx}$, which leads to irrelevant relative error evaluations. Once again, the implementation code has been successfully verified.

4.2 Example 2: Compression under prescribed displacement rate

Figure 18 shows a specimen with height H subjected to a constant imposed displacement rate ${\dot{\delta }}_z$ under plane strain conditions. The applied displacement $-{\delta }_z(t)$ is time-dependent and can be defined as

$-{\delta }_z(t)=-({\dot{\delta }}_{z}t+{\delta }_{z0})Y(t)$ (34)

where ${\delta }_{z0}$ stands for the vertical displacement imposed at the time $t=0^{+}$ and $Y(t)$ denotes the Heaviside step function at origin.

Figure 18
YZ-plane for specimen subjected to constant vertical displacement rate

Material will be considered homogeneous, isotropic and linear viscoelastic. The initial stress state will be considered as ${\underset{\_}{\underset{\_}{\sigma }}}^0=0$ and the strain ${\epsilon }_z(t)$ related to displacement $-{\delta }_z(t)$ can be expressed as

${\epsilon }_z(t)=-\frac{{\delta }_z(t)}{H}=-\left(\frac{{\dot{\delta }}_z}{H}t+\frac{{\delta }_{z0}}{H}\right)Y(t)=-({\dot{\epsilon }}_{z}t+{\epsilon }_{z0})Y(t)$ (35)

where ${\dot{\epsilon }}_z$ can be considered as the vertical strain rate and ${\epsilon }_{z0}$ as the instantaneous vertical strain.

As seen in sections 4.1 and 4.2, the problem solution involves the calculation of time-dependent stress and displacement fields $(\underset{\_}{\underset{\_}{\sigma }}(t),\underset{\_}{\xi }(\underset{\_}{X},t))$ by using basic concepts of continuum mechanics and viscoelasticity. Once the static and kinematic admissibility of the respective solution fields has been determined, the analytical solution tensorial fields are

$\begin{array}{l}\underset{\_}{\xi }(\underset{\_}{X},t)=\left[({\dot{\epsilon }}_{y}t+{\epsilon }_{y0})y\right]{\underset{\_}{e}}_y+\left[({\dot{\epsilon }}_{z}t+{\epsilon }_{z0})z\right]{\underset{\_}{e}}_z\\ \\ \underset{\_}{\underset{\_}{\epsilon }}(t)=\left[{\dot{\epsilon }}_{y}t+{\epsilon }_{y0}\right]{\underset{\_}{e}}_y\otimes {\underset{\_}{e}}_y+\left[{\dot{\epsilon }}_{z}t+{\epsilon }_{z0}\right]{\underset{\_}{e}}_z\otimes {\underset{\_}{e}}_z\\ \\ \underset{\_}{\underset{\_}{\sigma }}(t)=\left[({\dot{\epsilon }}_y+{\dot{\epsilon }}_z)\left({\displaystyle {\int }_0^t{R}_{1122}^{(s,t)}}ds\right)+({\epsilon }_{y0}+{\epsilon }_{z0})\left(R_{1122}^{(t)}\right)\right]{\underset{\_}{e}}_x\otimes {\underset{\_}{e}}_x+\\ \left[{\dot{\epsilon }}_y\left({\displaystyle {\int }_0^t{R}_{1111}^{(s,t)}}ds\right)+{\dot{\epsilon }}_z\left({\displaystyle {\int }_0^t{R}_{1122}^{(s,t)}}ds\right)+{\epsilon }_{y0}{R}_{1111}^{(t)}+{\epsilon }_{z0}{R}_{1122}^{(t)}\right]{\underset{\_}{e}}_y\otimes {\underset{\_}{e}}_y+\\ \left[{\dot{\epsilon }}_y\left({\displaystyle {\int }_0^t{R}_{1122}^{(s,t)}}ds\right)+{\dot{\epsilon }}_z\left({\displaystyle {\int }_0^t{R}_{1111}^{(s,t)}}ds\right)+{\epsilon }_{y0}{R}_{1122}^{(t)}+{\epsilon }_{z0}{R}_{1111}^{(t)}\right]{\underset{\_}{e}}_z\otimes {\underset{\_}{e}}_z\end{array}$ (36)

where $R_{1111}^{(s,t)}$ and $R_{1122}^{(s,t)}$ are relaxation tensor $ℝ(t)$ components. The instantaneous horizontal strain ${\epsilon }_{y0}$ can be obtained from the viscoelastic solution shown in section 4.1 for imposed displacement problem. The horizontal strain rate ${\dot{\epsilon }}_y$ is:

${\dot{\epsilon }}_y=\frac{-{\dot{\epsilon }}_z\left({\displaystyle {\int }_0^t{R}_{1122}^{(s,t)}}ds\right)-{\epsilon }_{y0}{R}_{1111}^{(t)}-{\epsilon }_{z0}{R}_{1122}^{(t)}}{\left({\displaystyle {\int }_0^t{R}_{1111}^{(s,t)}}ds\right)}$ (37)

It will be considered an imposed instantaneous displacement ${\delta }_{z0}=0.001m.$ and a displacement rate ${\dot{\delta }}_z={10}^{-8}m/s$ for calculations. In relation to Maxwell-Maxwell case, Figures (19) and (20) show the comparison between analytical solution defined in (36) and numerical response for vertical stress ${\sigma }_{zz}$ and stress ${\sigma }_{xx}$ perpendicular to YZ-plane.

Figure 19
Comparison of delayed responses for stress ${\sigma }_{zz}$ (Maxwell-Maxwell)

Figure 20
Comparison of delayed responses for stress ${\sigma }_{xx}$ (Maxwell-Maxwell)

Similarly, comparisons between analytical solution defined in (36) and numerical response for Burger-Maxwell case is illustrated in Figures (21) and (22).

Figure 21
Comparison of delayed responses for stress ${\sigma }_{zz}$ (Burger-Maxwell)

Figure 22
Comparison of delayed responses for stress ${\sigma }_{xx}$ (Burger-Maxwell)

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.

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

5.1 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$\gg$R. The massif material is homogeneous and isotropic. The excavation process will be regarded as instantaneous. The initial stress (geostatic) ${\underset{\_}{\underset{\_}{\sigma }}}^0=-p\underset{\_}{\underset{\_}{1}}$, due to self-weight of the massif, will be uniform in the region around the gallery. The variable $p$ is defined as $p=\gamma H$, where $\gamma$ represents the specific weight of the massif and $\underset{\_}{\underset{\_}{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.

Figure 24
Dimensions and boundary conditions for the geometric model

The analytical solution involving delayed response will be shown in polar coordinates due to the geometric characteristics of the tunnel cross-section. Position vector $\underset{\_}{X}$ does not depend on the angular coordinate and can be defined as $\underset{\_}{X}=r\cdot {\underset{\_}{e}}_r$, where $r$ is the radial distance. The problem solution involves calculation of stress and displacement fields $(\underset{\_}{\underset{\_}{\sigma }}(r,t),\underset{\_}{\xi }(r,t))$ by using basic concepts of continuum mechanics and viscoelasticity for $t>t_0{}^{+}$, where $t_0{}^{+}$ is considered as the instant right after excavation. The temporal evolution of solution fields is shown as follows

$\begin{array}{l}\underset{\_}{\xi }(r,t)=-\frac{p{R}^2}{r}M(t,t_0{}^{+}){\underset{\_}{e}}_r\\ \\ \underset{\_}{\underset{\_}{\sigma }}(r,t)=\left[p\left(\frac{R^2}{r^2}{Y}_{t_0^{+}}(t)-1\right)\right]{\underset{\_}{e}}_r\otimes {\underset{\_}{e}}_r+\left[-p\left(\frac{R^2}{r^2}{Y}_{t_0^{+}}(t)+1\right)\right]{\underset{\_}{e}}_{\theta }\otimes {\underset{\_}{e}}_{\theta }+\left(-p\right){\underset{\_}{e}}_z\otimes {\underset{\_}{e}}_z\end{array}$ (38)

where $Y_{t_0^{+}}(t)$ is the Heaviside function at $t_0{}^{+}$ and the modulus $M(t,t_0{}^{+})$ depends on the rheological models adopted by the matrix and fractures. This modulus can be defined as

$M(t,t_0{}^{+})={\mathcal{L}}_c^{-1}\left(\frac{1}{2{\mu }^{*}(p)}\right)$ (39)

where ${\mu }^{*}(p)$ is the shear relaxation modulus in Laplace-Carson domain.

The radial convergence is an important parameter related to the displacements in unlined cross-sections and can be expressed as

$U(t)=-\frac{\underset{\_}{\xi }(R,t)}{R}=pM(t,t_0{}^{+})$ (40)

It will be considered a depth $H=125.00m$, radius $R=2.50m$, specific weight $\gamma =24000N/m^3$ for calculations. The lateral and vertical pressure exerted by massif is $p=0.003GPa$. 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 $\epsilon$ and the number of fractures per unit volume $\mathcal{N}$. 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.

Figure 25
Finite element mesh for circular cross-section tunnel

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 ${\xi }_{rr}(r,t)$ and radial stress ${\sigma }_{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.

Figure 26
Comparison of radial convergence $U(t)$ (Maxwell-Maxwell)

Figure 27
Comparison of radial displacements as a function of $r$ (Maxwell-Maxwell)

Figure 28
Comparison of radial stresses as a function of $r$ (Maxwell-Maxwell)

Similarly, below is shown the comparison of analytical solution, defined in (38) and (40), with the numerical response for Burger-Maxwell case.

Figure 29
Comparison of radial convergence $U(t)$ (Burger-Maxwell)

Figure 30
Comparison of radial displacements as a function of $r$ (Burger-Maxwell)

Figure 31
Comparison of radial stresses as a function of $r$ (Burger-Maxwell)

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.

5.2 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 $u(R,t_0{}^{+})$ 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.

Figure 32
Displacement conditions and pressure on tunnel lining

As in section 5.1, position vector will be considered as $\underset{\_}{X}=r{\underset{\_}{e}}_r$ in polar coordinates. The problem solution involves calculation of stress and displacement fields $(\underset{\_}{\underset{\_}{\sigma }}(r,t),\underset{\_}{\xi }(r,t))$ for $t>t_0{}^{+}$. The temporal evolution of solution fields is

$\begin{array}{l}\underset{\_}{\xi }(r,t)=-\frac{p{R}^2}{2\mu (t_0{}^{+},t_0{}^{+})\cdot r}{Y}_{t_0^{+}}(t){\underset{\_}{e}}_r\\ \\ \underset{\_}{\underset{\_}{\sigma }}(r,t)=\left[p\left(\frac{R^2}{r^2}\frac{\mu (t,t_0{}^{+})}{\mu (t_0{}^{+},t_0{}^{+})}-1\right)\right]{\underset{\_}{e}}_r\otimes {\underset{\_}{e}}_r+\left[-p\left(\frac{R^2}{r^2}\frac{\mu (t,t_0{}^{+})}{\mu (t_0{}^{+},t_0{}^{+})}+1\right)\right]{\underset{\_}{e}}_{\theta }\otimes {\underset{\_}{e}}_{\theta }+(-p){\underset{\_}{e}}_z\otimes {\underset{\_}{e}}_z\end{array}$ (41)

where $\mu (t,t_0{}^{+})$ and $\mu (t_0{}^{+},t_0{}^{+})$ represent viscoelastic shear modulus evaluated at times $t$ and $t_0{}^{+}$, respectively. Pressure $q(t)$ exerted by massif on tunnel lining is defined as follows

$q(t)=-{\sigma }_{rr}(R,t)=p\left(1-\frac{\mu (t,t_0{}^{+})}{\mu (t_0{}^{+},t_0{}^{+})}\right)$ (42)

In relation to Maxwell-Maxwell case, the analytical solution, defined in (41) and (42), will be compared with the numerical response for pressure $q(t)$, radial displacement ${\xi }_{rr}(r,t)$ and radial stress ${\sigma }_{rr}(r,t)$.

Figure 33
Comparison of pressure $q(t)$ exerted on tunnel lining (Maxwell-Maxwell)

Figure 34
Comparison of radial displacements as a function of $r$ (Maxwell-Maxwell)

Figure 35
Comparison of radial stresses as a function of $r$ (Maxwell-Maxwell)

Similarly, below is shown the comparison of solutions for Burger-Maxwell case.

Figure 36
Comparison of pressure $q(t)$ exerted on tunnel lining (Burger-Maxwell)

Figure 37
Comparison of radial displacements as a function of $r$ (Burger-Maxwell)

Figure 38
Comparison of radial stresses as a function of $r$ (Burger-Maxwell)

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.

5.3 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 $\epsilon$, number of fractures per unit volume $\mathcal{N}$ and mechanical parameters from Tables 1 and 2 shown in section 4. Figure 39 displays a schematic representation of the geometric model and the horseshoe cross-section.

Figure 39
Geometric model dimensions and horseshoe cross-section analysis axes

Values adopted by variables shown in Figure 39 are: $H=270m$, $L1=45.68m$, $L2=16.00m$, $L3=61.68m$ and $L4=48.00m$. For this paper, it was arbitrarily chosen a semielliptical geometry for tunnel vault. The values for variables of horseshoe cross-section are: $a=6.10m$, $b=2.32m$ and $c=6.84m$. Figure 40 illustrates the stress and strain boundary conditions, as well as the finite element mesh generated for respective analyses.

It is known that $p=\gamma H$ for deep galleries, where massif specific weight is $\gamma =25kN/m^3$ and initial stress state before excavation is ${\underset{\_}{\underset{\_}{\sigma }}}^0=-p\underset{\_}{\underset{\_}{1}}$. 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 $({\xi }_{xx},{\xi }_{yy})$ and stresses $({\sigma }_{xx},{\sigma }_{yy})$ as a function of distances $r_x$ and $r_y$ for three different times. These times were defined in section 5.1.

Figure 40
Boundary conditions and finite element mesh for horseshoe cross-section tunnel

In relation to Maxwell-Maxwell case, time-dependent displacements and stresses in X-axis (point B) and Y-axis (point A) are shown below.

Figure 41
Vertical displacements ${\xi }_{yy}(r_y)$ in Y-axis (Maxwell-Maxwell)

Figure 42
Vertical stresses ${\sigma }_{yy}(r_y)$ in Y-axis (Maxwell-Maxwell)

Figure 43
Horizontal displacements ${\xi }_{xx}(r_x)$ in X-axis (Maxwell-Maxwell)

Figure 44
Horizontal stresses ${\sigma }_{xx}(r_x)$ in X-axis (Maxwell-Maxwell)

Similarly, below is shown time-dependent displacements and stresses for Burger-Maxwell case.

Figure 45
Vertical displacements ${\xi }_{yy}(r_y)$ in Y-axis (Burger-Maxwell)

Figure 46
Vertical stresses ${\sigma }_{yy}(r_y)$ in Y-axis (Burger-Maxwell)

Figure 47
Horizontal displacements ${\xi }_{xx}(r_x)$ in X-axis (Burger-Maxwell)

Figure 48
Horizontal stresses ${\sigma }_{xx}(r_x)$ in X-axis (Burger-Maxwell)

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.

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

References

• Aguiar, C.B., Maghous, S., (2018). Micromechanical approach to effective viscoelastic properties of micro‐fractured geomaterials. International Journal for Numerical and Analytical Methods in Geomechanics 42(16): 2018-2046
• Bart, M. (2000). Contribution à la modélisation du comportement hydromécanique des massifs rocheux avec fractures, Ph.D. Thesis (in French), Université Lille 1, France.
• Budiansky B, O’Connell RJ (1976). Elastic moduli of a cracked solid.International Jounal ofSolids Structures 12: 81-97.
• Couto, E.C. (2011). Um modelo tridimensional para túneis escavados em rocha reforçada por tirantes passivos, Ph.D. Thesis (in Portuguese), Federal University of Rio Grande do Sul, Brazil.
• Le QV. (2008). Modélisation multi‐échelle des matériaux viscoélastiques hétérogènes. Application à l'identification et à l'estimation du fluage propre de béton d'enceintes de centrales nucléaires. Ph. D. Thesis (in French), Université Paris‐Est, France.
• Maghous, S., Dormieux, L., Kondo, D., Shao, J. F. (2011). Micromechanics approach to poroelastic behavior ofa jointed rock. International Journal for Numerical and Analytical Methods in Geomechanics 37: 111-129.
• Maghous, S., Lorenci, G., Bittencourt, E. (2014). Effective poroelastic behavior of a jointed rock. Mechanics Research Communications 59: 64-69.
• Nguyen, S. T., Dormieux, L., Le Pape, Y., Sanhuja, J. (2010). Crack propagation in viscoelastic structures:Theoretical and numerical analyses. Computational Materials Sciences 50: 83-91.
• Nguyen, S. T., Jeannin, L., Dormieux, L., Renard, F. (2013). Fracturing of viscoelastic geomaterials and application to sedimentary layered rocks. Mechanics Research Communications 49: 50-56.
• Nguyen, S. T., Dormieux, L. (2014). Propagation of micro-cracks in viscoelastic materials: Analytical andnumerical methods. International Journal of Damage Mechanics 24(4): 562-581
• Zaoui A. (2002). Continuum micromechanics: survey. Journal of Engineering Mechanics ASCE 128(8): 808‐816.

Edited by

Editor:

Marco L. Bittencourt.

Publication Dates

• Publication in this collection
  25 Nov 2019
• Date of issue
  2019

History

• Received
  01 July 2019
• Reviewed
  14 Oct 2019
• Accepted
  04 Nov 2019
• Published
  06 Nov 2019