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

Ricardo Haniel Moran Ramirez Cássio B. Aguiar Eduardo Bittencourt Samir Maghous About the authors

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, 1976Budiansky B, O’Connell RJ (1976). Elastic moduli of a cracked solid.International Jounal ofSolids Structures 12: 81-97.). 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. (2010Nguyen, 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., 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 (2018Aguiar, 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), 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 (2018Aguiar, 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) 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 (2018Aguiar, 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) 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, 2008Le 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.) 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 (2018Aguiar, 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).

2.1 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, 2002Zaoui A. (2002). Continuum micromechanics: survey. Journal of Engineering Mechanics ASCE 128(8): 808‐816.). 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): ξ_(x_)=__x_x_Ω, where __ represents the macroscopic strain (loading parameter), ξ_ is the displacement filed and 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 hom components. This determination relies upon the solution to a viscoelastic concentration problem stated on the REV under applied loading __ (Aguiar and Maghous, 2018Aguiar, 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). In such conditions, the homogenized constitutive law for a fractured viscoelastic medium can be defined as

Σ _ _ = hom _ _ = hom ( t τ ) : ˙ _ _ ( τ ) d τ w i t h hom = s A (1)

where Σ__ is the macroscopic stress, Ais the strain concentration tensor associated with the viscoelastic problem defined on the REV. The fourth order tensor Arepresents 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 tensorhom, 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.

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

After tensor hom 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 hom=Lc1(hom), which was developed by Aguiar and Maghous (2018Aguiar, 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). This procedure was conceived to work with the tensor hom 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 hom is generically referred to as Rhom. This function can be expressed as

R hom ( p ) = A ( p ) B ( p ) = i = 0 n a i p i 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 ai and bi 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 Rhom(p) changed to the laplace domain:

F ( p ) = R hom ( p ) p = a 0 p b 0 + C ( p ) B ( p ) w i t h C ( p ) = 1 b 0 i = 1 n [ ( b 0 a i a 0 b i ) p i 1 ] (3)

The inverse Laplace transform of F(p) can be named as Rhom(t) and the corresponding expression is

R hom ( t ) = [ a 0 b 0 + i = 1 n D i e R i t ] Y ( t ) w i t h D i = [ C ( p ) B ( p ) / p ] p = R i (4)

where Rhom(t) represents the non-zero component of tensor hom(t), Y(t) represent the Heaviside function at origin and the scalar Ri is the kth root of the polynomial B(p). In the context of macroscopic isotropy, the generic function Rhom(t) refers to either bulk khom(t) or shear μ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).

hom ( t ) = 3 k hom ( t ) J + 2 μ hom ( t ) K k hom ( t ) = [ a 0 k b 0 k + i = 1 n D i k e R i k t ] Y ( t ) , μ hom ( t ) = [ a 0 μ b 0 μ + i = 1 n D i μ e R i μ t ] Y ( t ) (5)

where the subscripts k and μ indicate that the parameters Riand Di are referring to the relaxation behavior in bulk and relaxation behavior in shear, respectively. The fourth-order tensors J and K are the spherical and deviatoric tensors, respectively. These tensors can be defined as J=131__1__ and K=IJ, where 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, 2018Aguiar, 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), 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, 2018Aguiar, 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)

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 khom and μhom(excerpt from Aguiar and Maghous, 2018Aguiar, 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)

As seen in Eq. (5), for an isotropic distribution of randomly oriented fractures, tensor hom(t) depends on the moduli khom and μhom. In this way, taking advantage of the elastic-viscoelastic correspondence principle, viscoelastic moduli khom 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).

k hom = k s 1 + ε M k , μ hom = μ s 1 + ε M μ (6)

where ε=Na3 is the fracture density parameter, N is the number of fractures per unit volume and a is the radius of the oblate fracture. The dimensionless functions Mk(ks,μs,akn,akt) and Mμ(ks,μs,akn,akt) are defined as

M k = 4 3 π k s / μ s π κ 1 + a k n / μ s , M μ = 16 π κ 4 15 6 κ 2 + 4 κ 3 + 9 π κ 4 ( 3 κ 1 + κ 4 ) ( 3 π κ 1 κ 4 + κ 2 ) ( 4 κ 3 + 9 π κ 4 ( κ 1 + κ 4 ) ) κ 1 = 3 k s + μ s 3 k s + 4 μ s , κ 2 = 3 k n a 3 k s + 4 μ s , κ 3 = 3 k t a 3 k s + 4 μ s , κ 4 = μ s 3 k s + 4 μ s (7)

where parameters (ks,μs) and (kn,kt) depend on the rheological model adopted for the viscoelastic behavior of the constituents (matrix and fractures). As mentioned above, having defined the moduli khom and μhom, it only remains to apply the inverse Laplace-Carson transform procedure in order to obtain the moduli in time domain: khom(t) and μ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 σ__ and ε__, respectively.

The relaxation tensor hom(t) for an isotropic distribution of microfractures, as well as the corresponding viscoelastic moduli khom(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

hom = [ ( k hom + 4 3 μ hom ) ( k hom 2 3 μ hom ) ( k hom 2 3 μ hom ) 0 0 0 ( k hom 2 3 μ hom ) ( k hom + 4 3 μ hom ) ( k hom 2 3 μ hom ) 0 0 0 ( k hom 2 3 μ hom ) ( k hom 2 3 μ hom ) ( k hom + 4 3 μ hom ) 0 0 0 0 0 0 2 μ hom 0 0 0 0 0 0 2 μ hom 0 0 0 0 0 0 2 μ hom ] (8)

Once the relaxation tensor hom(t) is computed, the macroscopic stress tensor can be determined, according to the constitutive law, as follows

σ _ _ ( t ) = hom ( t ) : ε _ _ ( 0 ) + 0 t hom ( s , t ) : ε _ _ s ( s ) d s (9)

where ε__(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 Δt and t0=0. For the i+1 time step, the corresponding time can be defined as

t i + 1 = t i + Δ t = ( i + 1 ) Δ t w i t h i = 0, ... , n 1 (10)

An incremental procedure is applied in the integral equation in (9). Thus, the stress tensor at time ti+1, is

σ _ _ ( t i + 1 ) = hom ( t i + 1 ) : ε _ _ ( 0 ) + k = 0 i [ t k t k + 1 hom ( s , t i + 1 ) : ε _ _ s ( s ) d s ] (11)

A linear variation is assumed for the strain in the interval where the variable of integration s[tk,tk+1]. This assumption allows to express the derivative ε__/s as

ε _ _ s = Δ ε _ _ k + 1 Δ t = ε _ _ ( t k + 1 ) ε _ _ ( t k ) Δ t (12)

In order to update the stress tensor equation, (12) is replaced in (11). In a summarized form, the stress tensor σ__(ti+1) can be expressed as

σ _ _ ( t i + 1 ) = K ˜ T : Δ ε _ _ i + 1 + σ _ _ h ( t ) (13)

where K˜T and σ__h(t) are fourth-order tensors defined as

K ˜ T = 1 Δ t ( t i t i + 1 hom ( s , t i + 1 ) d s ) , σ _ _ h ( t ) = σ _ _ 0 ( t i + 1 ) + k = 0 i 1 [ A ˜ k + 1 : Δ ε _ _ k + 1 ] σ _ _ 0 ( t i + 1 ) = hom ( t i + 1 ) : ε _ _ ( 0 ) , A ˜ k + 1 = 1 Δ t ( t k t k + 1 hom ( s , t i + 1 ) d s ) (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 hom(t) components and it can be expressed in a generic way as

R i j k l ( s , t ) = R + α = 1 p R α e ( t s τ α ) w i t h τ α = η α R α (15)

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 K˜T and A˜k+1 depend on tensor hom, their components can be conveniently expressed as

K i j k l T = R + α = 1 p [ R α ( Δ t τ α ) ( 1 e ( Δ t τ α ) ) ] , A i j k l k + 1 = R + α = 1 p [ R α ( Δ t τ α ) ( e ( i k ) ( Δ t τ α ) e ( i k + 1 ) ( Δ t τ α ) ) ] (16)

It is known that in an isotropic distribution of fractures, for the non-zero components of tensor hom, the following applies

R 1111 hom = R 2222 hom = R 3333 hom , R 1212 hom = R 1313 hom = R 2323 hom R 1122 hom = R 1133 hom = R 2211 hom = R 2233 hom = R 3311 hom = R 3322 hom (17)

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

K 1111 T = a 0 ω b 0 ω + α = 1 n ω ( D α ω ( R α ω Δ t ) ( 1 e R α ω Δ t ) ) + 4 3 [ a 0 μ b 0 μ + α = 1 n μ ( D α μ ( R α μ Δ t ) ( 1 e R α μ Δ t ) ) ] (18)

K 1122 T = a 0 ω b 0 ω + α = 1 n ω ( D α ω ( R α ω Δ t ) ( 1 e R α ω Δ t ) ) 2 3 [ a 0 μ b 0 μ + α = 1 n μ ( D α μ ( R α μ Δ t ) ( 1 e R α μ Δ t ) ) ] (19)

K 1111 T = a 0 ω b 0 ω + α = 1 n ω ( D α ω ( R α ω Δ t ) ( 1 e R α ω Δ t ) ) + 4 3 [ a 0 μ b 0 μ + α = 1 n μ ( D α μ ( R α μ Δ t ) ( 1 e R α μ Δ t ) ) ] (20)

Similarly, for tensor A˜k+1, the corresponding non-zero components are

A 1111 k + 1 = a 0 ω b 0 ω + α = 1 n ω ( D α ω ( R α ω Δ t ) ( e ( i k ) R α ω Δ t e ( i k + 1 ) R α ω Δ t ) ) + 4 3 [ a 0 μ b 0 μ + α = 1 n μ ( D α μ ( R α μ Δ t ) ( e ( i k ) R α μ Δ t e ( i k + 1 ) R α μ Δ t ) ) ] (21)

A 1122 k + 1 = a 0 ω b 0 ω + α = 1 n ω ( D α ω ( R α ω Δ t ) ( e ( i k ) R α ω Δ t e ( i k + 1 ) R α ω Δ t ) ) 2 3 [ a 0 μ b 0 μ + α = 1 n μ ( D α μ ( R α μ Δ t ) ( e ( i k ) R α μ Δ t e ( i k + 1 ) R α μ Δ t ) ) ] (22)

A 1212 k + 1 = 2 [ a 0 μ b 0 μ + α = 1 n μ ( D α μ ( R α μ Δ t ) ( e ( i k ) R α μ Δ t e ( i k + 1 ) R α μ Δ t ) ) ] (23)

where ω corresponds to the bulk relaxation modulus and μ 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 ti+1=t1 and i=0, the stress tensor σ__TS can be defined as:

σ _ _ 1 = K ˜ T : Δ ε _ _ 1 + 1 : ε _ _ 0 (24)

where Δε__1=ε__1ε__0 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).

[ σ 11 1 σ 22 1 σ 33 1 σ 12 1 σ 13 1 σ 23 1 ] = [ K 1111 T Δ ε 11 1 + K 1122 T Δ ε 22 1 + R 1111 1 ε 11 0 + R 1122 1 ε 22 0 K 1111 T Δ ε 22 1 + K 1122 T Δ ε 11 1 + R 1111 1 ε 22 0 + R 1122 1 ε 11 0 K 1122 T Δ ε 11 1 + K 1122 T Δ ε 22 1 + R 1122 1 ε 11 0 + R 1122 1 ε 22 0 K 1212 T Δ ε 12 1 + R 1212 1 ε 12 0 0 0 ] (25)

For TS = 2, where ti+1=t2 and i=1, the stress tensor σ__TS can be defined as

σ _ _ 2 = K ˜ T : Δ ε _ _ 2 + 2 : ε _ _ 0 + A ˜ 1 : Δ ε _ _ 1 (26)

where Δε__2=ε__2ε__1 and 2 is the homogenized relaxation tensor evaluated at TS = 2. A˜1 is the tensor A˜k+1 evaluated at TS=1 (k=0). The matrix form of tensor σ__2 is found in (27).

[ σ 11 2 σ 22 2 σ 33 2 σ 12 2 σ 13 2 σ 23 2 ] = [ K 1111 T Δ ε 11 2 + K 1122 T Δ ε 22 2 + R 1111 2 ε 11 0 + R 1122 2 ε 22 0 + A 1111 1 Δ ε 11 1 + A 1122 1 Δ ε 22 1 K 1111 T Δ ε 22 2 + K 1122 T Δ ε 11 2 + R 1111 2 ε 22 0 + R 1122 2 ε 11 0 + A 1111 1 Δ ε 22 1 + A 1122 1 Δ ε 11 1 K 1122 T Δ ε 11 2 + K 1122 T Δ ε 22 2 + R 1122 2 ε 11 0 + R 1122 2 ε 22 0 + A 1122 1 Δ ε 11 1 + A 1122 1 Δ ε 22 1 K 1212 T Δ ε 12 2 + R 1212 2 ε 12 0 + A 1212 1 Δ ε 12 1 0 0 ] (27)

For TS ≥ 3, where ti+1t3, i=TS1 and k[0,i1], the stress tensor σ__TS can be defined as:

σ _ _ T S = K ˜ T : Δ ε _ _ T S + T S : ε _ _ 0 + k = 0 i 1 [ A ˜ k + 1 : Δ ε _ _ k + 1 ] (28)

The matrix form of tensor σ__TS for TS ≥ 3 is presented in (29).

[ σ 11 T S σ 22 T S σ 33 T S σ 12 T S σ 13 T S σ 23 T S ] = [ K 1111 T Δ ε 11 T S + K 1122 T Δ ε 22 T S + R 1111 T S ε 11 0 + R 1122 T S ε 22 0 + S A E 11 K 1111 T Δ ε 22 T S + K 1122 T Δ ε 11 T S + R 1111 T S ε 22 0 + R 1122 T S ε 11 0 + S A E 22 K 1122 T Δ ε 11 T S + K 1122 T Δ ε 22 T S + R 1122 T S ε 11 0 + R 1122 T S ε 22 0 + S A E 33 K 1212 T Δ ε 12 T S + R 1212 T S ε 12 0 + S A E 12 0 0 ] (29)

where SAE11, SAE22, SAE33 and SAE12 are defines as

S A E 11 = k = 0 i 1 ( A 1111 k + 1 Δ ε 11 k + 1 + A 1122 k + 1 Δ ε 22 k + 1 ) , S A E 22 = k = 0 i 1 ( A 1111 k + 1 Δ ε 22 k + 1 + A 1122 k + 1 Δ ε 11 k + 1 ) S A E 33 = k = 0 i 1 ( A 1122 k + 1 Δ ε 11 k + 1 + A 1122 k + 1 Δ ε 22 k + 1 ) , S A E 12 = k = 0 i 1 ( A 1212 k + 1 Δ ε 12 k + 1 ) (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 (2008Le 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.) for viscoelastic moduli of concrete, in Bart (2000Bart, M. (2000). Contribution à la modélisation du comportement hydromécanique des massifs rocheux avec fractures, Ph.D. Thesis (in French), Université Lille 1, France.) for elastic fracture behavior and the viscosity parameters of fractures were arbitrarily chosen according to Aguiar and Maghous (2018Aguiar, 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). The fracture density parameter is ε=0.20 and the number of fractures per unit volume 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

Table 2
Fractures mechanical parameters for Maxwell model

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 ke,Ms and ke,Ks, whereas that associated in shear are μe,Ms and μe,Ks. Viscosity parameters associated to a dashpot in bulk are denoted by kv,Ms and kv,Ks, whereas that associated in shear are μv,Ms and μv,Ks. As regards the fractures behavior, the elastic parameter associated with a spring under normal stress is denoted by ke,Mn, whereas that associated under tangential stress is ke,Mt. The viscosity parameter associated to a dashpot under normal stress is denoted by kv,Mn, whereas that associated under tangential stress is kv,Mt. 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 Ei 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.

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 σY(t), whereas Figure 7B shows a vertical imposed displacement δzY(t)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 σ__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 (σ__(t),ξ_(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

ξ _ ( X _ , t ) = [ σ ( 2 μ ' ( t ) 3 k ' ( t ) ) 4 μ ' ( t ) ( 3 k ' ( t ) + μ ' ( t ) ) ] y e _ y + [ σ ( 3 k ' ( t ) + 4 μ ' ( t ) ) 4 μ ' ( t ) ( 3 k ' ( t ) + μ ' ( t ) ) ] z e _ z ε _ _ ( t ) = [ σ ( 2 μ ' ( t ) 3 k ' ( t ) ) 4 μ ' ( t ) ( 3 k ' ( t ) + μ ' ( t ) ) ] e _ y e _ y + [ σ ( 3 k ' ( t ) + 4 μ ' ( t ) ) 4 μ ' ( t ) ( 3 k ' ( t ) + μ ' ( t ) ) ] e _ z e _ z σ _ _ ( t ) = [ σ ( 3 k ' ( t ) 2 μ ' ( t ) ) 2 ( 3 k ' ( t ) + μ ' ( t ) ) ] e _ x e _ x σ e _ z e _ z (31)

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

k ' ( t ) = 1 L c 1 ( 1 / k ( p ) hom ) , μ ' ( t ) = 1 L c 1 ( 1 / μ ( p ) hom ) (32)

where k(p)hom and μ(p)hom are the homogenized bulk and shear moduli in Laplace-Carson domain, respectively. Lc1 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 σ=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 εzz and horizontal strain εyy.

Figure 8
Comparison of delayed responses for vertical strain εzz (Maxwell-Maxwell)

Figure 9
Comparison of delayed responses for horizontal strain ε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 εzz (Burger-Maxwell)

Figure 11
Comparison of delayed responses for horizontal strain ε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

ξ _ ( X _ , t ) = [ ( 3 k ( t ) 2 μ ( t ) 3 k ( t ) + 4 μ ( t ) ) ε z y ] e _ y + ε z z e _ z ε _ _ ( t ) = [ ( 3 k ( t ) 2 μ ( t ) 3 k ( t ) + 4 μ ( t ) ) ε z ] e _ y e _ y + ε z e _ z e _ z σ _ _ ( t ) = [ 2 μ ( t ) ( 3 k ( t ) 2 μ ( t ) ) 3 k ( t ) + 4 μ ( t ) ε z ] e _ x e _ x + [ 4 μ ( t ) ( 3 k ( t ) + μ ( t ) ) 3 k ( t ) + 4 μ ( t ) ε z ] e _ z e _ z (33)

where variables y and z are coordinates of the position vector X_. The strain caused by imposed displacement is denoted by εz=δz/H. The bulk and shear moduli in time domain are denoted by k(t) and μ(t), respectively. The numerical expressions of k(t) and μ(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 δ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 σzz and stress σxx perpendicular to YZ-plane.

Figure 13
Comparison of delayed responses for stress σzz (Maxwell-Maxwell)

Figure 14
Comparison of delayed responses for stress σ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 σzz (Burger-Maxwell)

Figure 16
Comparison of delayed responses for stress σ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 σ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 δ˙z under plane strain conditions. The applied displacement δz(t) is time-dependent and can be defined as

δ z ( t ) = ( δ ˙ z t + δ z 0 ) Y ( t ) (34)

where δ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 σ__0=0 and the strain εz(t) related to displacement δz(t) can be expressed as

ε z ( t ) = δ z ( t ) H = ( δ ˙ z H t + δ z 0 H ) Y ( t ) = ( ε ˙ z t + ε z 0 ) Y ( t ) (35)

where ε˙z can be considered as the vertical strain rate and ε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 (σ__(t),ξ_(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

ξ _ ( X _ , t ) = [ ( ε ˙ y t + ε y 0 ) y ] e _ y + [ ( ε ˙ z t + ε z 0 ) z ] e _ z ε _ _ ( t ) = [ ε ˙ y t + ε y 0 ] e _ y e _ y + [ ε ˙ z t + ε z 0 ] e _ z e _ z σ _ _ ( t ) = [ ( ε ˙ y + ε ˙ z ) ( 0 t R 1122 ( s , t ) d s ) + ( ε y 0 + ε z 0 ) ( R 1122 ( t ) ) ] e _ x e _ x + [ ε ˙ y ( 0 t R 1111 ( s , t ) d s ) + ε ˙ z ( 0 t R 1122 ( s , t ) d s ) + ε y 0 R 1111 ( t ) + ε z 0 R 1122 ( t ) ] e _ y e _ y + [ ε ˙ y ( 0 t R 1122 ( s , t ) d s ) + ε ˙ z ( 0 t R 1111 ( s , t ) d s ) + ε y 0 R 1122 ( t ) + ε z 0 R 1111 ( t ) ] e _ z e _ z (36)

where R1111(s,t) and R1122(s,t) are relaxation tensor (t) components. The instantaneous horizontal strain εy0 can be obtained from the viscoelastic solution shown in section 4.1 for imposed displacement problem. The horizontal strain rate ε˙y is:

ε ˙ y = ε ˙ z ( 0 t R 1122 ( s , t ) d s ) ε y 0 R 1111 ( t ) ε z 0 R 1122 ( t ) ( 0 t R 1111 ( s , t ) d s ) (37)

It will be considered an imposed instantaneous displacement δz0=0.001m. and a displacement rate δ˙z=108m/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 σzz and stress σxx perpendicular to YZ-plane.

Figure 19
Comparison of delayed responses for stress σzz (Maxwell-Maxwell)

Figure 20
Comparison of delayed responses for stress σ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 σzz (Burger-Maxwell)

Figure 22
Comparison of delayed responses for stress σ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 HR. The massif material is homogeneous and isotropic. The excavation process will be regarded as instantaneous. The initial stress (geostatic) σ__0=p1__, 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.

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 X_ does not depend on the angular coordinate and can be defined as X_=re_r, where r is the radial distance. The problem solution involves calculation of stress and displacement fields (σ__(r,t),ξ_(r,t)) by using basic concepts of continuum mechanics and viscoelasticity for t>t0+, where t0+ is considered as the instant right after excavation. The temporal evolution of solution fields is shown as follows

ξ _ ( r , t ) = p R 2 r M ( t , t 0 + ) e _ r σ _ _ ( r , t ) = [ p ( R 2 r 2 Y t 0 + ( t ) 1 ) ] e _ r e _ r + [ p ( R 2 r 2 Y t 0 + ( t ) + 1 ) ] e _ θ e _ θ + ( p ) e _ z e _ z (38)

where Yt0+(t) is the Heaviside function at t0+ and the modulus M(t,t0+) depends on the rheological models adopted by the matrix and fractures. This modulus can be defined as

M ( t , t 0 + ) = L c 1 ( 1 2 μ * ( p ) ) (39)

where μ*(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 ) = ξ _ ( R , t ) R = p M ( t , t 0 + ) (40)

It will be considered a depth H=125.00m, radius R=2.50m, specific weight γ=24000N/m3 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 ε and the number of fractures per unit volume 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 ξ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.

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,t0+) 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 X_=re_r in polar coordinates. The problem solution involves calculation of stress and displacement fields (σ__(r,t),ξ_(r,t)) for t>t0+. The temporal evolution of solution fields is

ξ _ ( r , t ) = p R 2 2 μ ( t 0 + , t 0 + ) r Y t 0 + ( t ) e _ r σ _ _ ( r , t ) = [ p ( R 2 r 2 μ ( t , t 0 + ) μ ( t 0 + , t 0 + ) 1 ) ] e _ r e _ r + [ p ( R 2 r 2 μ ( t , t 0 + ) μ ( t 0 + , t 0 + ) + 1 ) ] e _ θ e _ θ + ( p ) e _ z e _ z (41)

where μ(t,t0+) and μ(t0+,t0+) represent viscoelastic shear modulus evaluated at times t and t0+, respectively. Pressure q(t) exerted by massif on tunnel lining is defined as follows

q ( t ) = σ r r ( R , t ) = p ( 1 μ ( t , t 0 + ) μ ( t 0 + , t 0 + ) ) (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 ξrr(r,t) and radial stress σ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 (2011Couto, 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.) 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 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=γH for deep galleries, where massif specific weight is γ=25kN/m3 and initial stress state before excavation is σ__0=p1__. 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 (ξxx,ξyy) and stresses (σxx,σyy) as a function of distances rx and ry 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 ξyy(ry) in Y-axis (Maxwell-Maxwell)

Figure 42
Vertical stresses σyy(ry) in Y-axis (Maxwell-Maxwell)

Figure 43
Horizontal displacements ξxx(rx) in X-axis (Maxwell-Maxwell)

Figure 44
Horizontal stresses σxx(rx) in X-axis (Maxwell-Maxwell)

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

Figure 45
Vertical displacements ξyy(ry) in Y-axis (Burger-Maxwell)

Figure 46
Vertical stresses σyy(ry) in Y-axis (Burger-Maxwell)

Figure 47
Horizontal displacements ξxx(rx) in X-axis (Burger-Maxwell)

Figure 48
Horizontal stresses σxx(rx) 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 (2018Aguiar, 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) 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.

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
Associação Brasileira de Ciências Mecânicas Av. Rio Branco, 124/14º andar, 20040-001 Rio de Janeiro RJ Brasil, Tel.: (55 21) 2221 0438 - Rio de Janeiro - RJ - Brazil
E-mail: abcm@abcm.org.br