Non linear elasto-plastic analysis of cylindrical cavity in rock mass using a Hoek-Brown criterion

Análise não linear elastoplástica de cavidades cilíndricas em maciços rochosos usando o critério de Hoek-Brown

Abstract

This paper aims to present an elastic, perfectly plastic, constitutive model based on the Hoek-Brown failure criterion and with non-associative plasticity. The objective is to apply the model to the non-linear analysis of geotechnical problems like excavations in rock mass. The computational implementation was carried out with a computational program called ANLOG (Non-Linear Analysis of Geotechnical Problem) system based on a displacement formulation of the finite element method. Due to the non-linear nature of the constitutive model, the study adopts an incremental iterative Newton-Raphson procedure with automatic load increments to guarantee the global level equilibrium. In addition, to guarantee the consistency condition at the local level, the study adopts, for the stress integration, an explicit algorithm with automatic sub-increments of strain. To validate the computational implementation and applicability of the numerical model, the study uses theoretical results to compare with ones obtained with the numerical simulation of cylindrical cavity in rock mass.

Keywords:
Hoek-Brown failure criterion; finite element method; elastoplasticity; stress integration algorithm; cylindrical cavity; tunnel; rock mass

Resumo

Esse artigo apresenta um modelo constitutivo elástico perfeitamente plástico com plasticidade associada e com base no critério de resistência de Hoek-Brown. O objetivo é aplicar esse modelo para análise não linear de problemas geotécnicos como escavações em maciços rochosos. As implementações computacionais foram realizadas no sistema ANLOG (Análise não linear de obras geotécnicas) com base na formulação em deslocamento do método dos elementos finitos. Devido à natureza não linear do modelo constitutivo, o estudo adota um procedimento incremental iterativo do tipo Newton-Raphson com incrementos automáticos de modo a garantir o equilíbrio em nível global. Além disto, para garantir a condição de consistência em nível local, o estudo adota um esquema de explicito de integração de tensão com subincrementos automáticos de deformação. Para validar as implementações computacionais e a aplicabilidade do modelo numérico gerado, o estudo usa os resultados da simulação numérica de uma cavidade cilíndrica em maciços rochosos.

Palavras-chave:
Critério de ruptura de Hoek-Brown; método dos elementos finitos; elastoplasticidade; algoritmo de integração de tensão; cavidade cilíndrica, túnel; maciço rochoso

1. Introduction

The application of the finite element method (FEM), which considers a continuous media, to analyze the mechanical behavior of a rock mass has been restricted to hard rock or non-fractured rock mass. Due to the increasing number of geotechnical works carried out on fractured rock masses, it has become necessary to use a constitutive model that takes into account the geological condition of a rock mass. Such a model is capable of providing more realistic results when using a displacement formulation of FEM.

In the early 1980's, the Hoek-Brown failure criterion was developed for hard rock (Hoek and Brown 1980). Since then, several versions have been published in order to include the influence of geological conditions on the failure parameter of rock masses (Hoek and Brown 1988; Hoek et al 1992, 1995, 2002).

The use of a Hoek-Brown failure criterion, as the yield function in an elastic-plastic analysis, leads to the application of an incremental iterative procedure at the global level of a FEM analysis and the application of a stress integration scheme at the local level (Sharan 2003; Sharan 2005; Clausen and Dumkilde 2008; Wang and Yin 2011).

2. Hoek-Brown elasto-plastic model formulation

The equilibrium equations of a mechanical problem, in static condition, describe a non-linear equation system when adopting an elasto-plastic stress-strain-strength relationship (Teixeira et al2012). During a given equilibrium path, the variation on the displacement, strain, and stress fields depends on the stress and strain levels and their history through the equilibrium path.

Based on the displacement finite element formulation, the equilibrium equations can be written as:

Where, _ext is the external nodal force vector that represents the global arrangement of the element external nodal force vector ^e _ext defined as:

in which ^e _s represents the parcel of external nodal force due to surface load; ^e _b represents the parcel of external nodal force due to body force, ^e _; and represents the parcel of external nodal force due to non-null prescribed displacements, δ. _int is the internal nodal force vector that represents the global arrangement of the element internal nodal force vector ^e _int equivalent to the stress state σ in a given element that is defined as:

is the cinematic matrix which depends on the strain-displacement relationship.

Due to the non-linear nature of the equation system represented by Equation (1), an incremental-iterative procedure should be used in order to obtain the displacement, strain, and stress fields. Then, starting from a given equilibrium configuration n, where the stress and strain states are known, a predicted incremental solution in terms of the global displacement (Δ^o _n) is obtained. This predicted approximation should be corrected by successive iteration (δ Δn+1 (Crisfield 1991). In this strategy, the problem solution is obtained by updating the nodal displacement vector () until reaching a new equilibrium configuration ) in each new equilibrium configuration, by doing:

Where,

iter is the necessary iterative cycle number to reach convergence at the current step, while Δλ is the increment of load factor, which can be automatically defined starting from the initial trial provided by the user (Nogueira 1998; Oliveira 2006; Simão 2014). _ep is the global stiffness matrix that represents the global arrangement of the element stiffness matrix ^e _ep defined by:

where, _ep is the elasto-plastic constitutive matrix which depends on the stress-strain-strength relationship. According to the Modified Newton-Raphson iterative scheme the global stiffness matrix is kept constant during the iterative cycles.

The vector _ext ^k represents the external nodal force applied at each load step and kept constant throughout the iterative cycles, according to the Newton-Raphson iterative scheme. This vector is updated at the beginning of a given step load by:

where _{ext n}is the external nodal force vector at a given equilibrium configuration n. The internal nodal force vector _int ^k is evaluated at each iterative cycle depending on the stress state evaluated at this iterative cycle, σ ^k.

At the end of each iterative cycle, a convergence state of the solution is verified by using a criterion that relates the Euclidian norm of the unbalance nodal force vector with the Euclidian norm of the external nodal force vector. Thus, for a given tolerance and at each increment, the iterative scheme ensures the overall balance by satisfying the compatibility conditions, boundary conditions and constitutive relationships.

This iterative scheme involves the stress state evaluation at each iterative cycle. Then, in each element, the stress vector σ ^k is obtained by:

Where,

and, Δ ^k is the incremental displacement vector at element level and updated at the current iterative cycle

By adopting linear elastic, perfectly plastic (which is free of hardening during the plastic flow) and with the associated plasticity (in which the potential plastic and yield functions are the same) constitutive model, the elasto-plastic constitutive matrix can be written as:

where, _e is the elastic constitutive matrix which depends on the young modulus (E) and Poisson coefficient (ν). The vector a is the gradient of the yield function (F) which depends on the failure criterion used.

This paper adopts the generalized Hoek-Brown failure criterion (Hoek et al 2002) written in terms of principal stress, σ₁ and σ₃, as:

In which σ_ci is the uniaxial compressive strength of intact rock, and m_b, a and s are constants defined as:

Where, m_i is a constant of Hoek-Brown criterion intact rock; D is a disturbance coefficient which varies from 0.0 for undisturbed in situ rock mass to 1.0 for very disturbed rock mass, and GSI is the Geological Strength Index, which takes into account the geological condition of the rock mass.

As no hardening occurs in a perfectly plastic constitutive model and the yield concept merges with the failure concept, the yield function can be written as:

or, by using the Hoek-Brown failure criterion (Equation 15):

where, I is the first invariant of the stress tensor; I_2D is the second invariant of the deviator stress tensor, and θ is the Lode angle that depends on the second and the third (I_3D) invariant of deviator stress tensor (Owen and Hinton 1989).

In theory, at each increment and for a selected tolerance, the iterative scheme satisfies the global equilibrium equations, the compatibility and boundary conditions, and the constitutive equations. However, the constitutive equation integration is not trivial, even if one knows the incremental strain magnitude on each iterative cycle, what is still unknown is the way it varies across the incremental path. It is necessary then to use an accurate stress integration algorithm.

It is possible, by adopting the explicit stress integration scheme as suggested by Sloan et al (2001), to increase the precision of the stress calculation. According to these authors the increment of stress can be divided into two parcels: elastic, Δσ _e, and elastic plastic, Δσ _ep, such as:

where, α is a scalar that varies from 0 to unity. For α = 0 the strain increment generates an elastic plastic stress increment. For α = 1 the strain increment generates a purely elastic stress increment. The α scalar value is obtained by solving iteratively:

where FTOL is the tolerance suggested by Sloan et al, (2001) as 10^-5. Once defined, the α scalar, the stress state upon the yield surface, is updated according to:

The increment of elastic plastic stress is obtained by:

where ΔT^j is a scalar known as increment of pseudo-time (T) that varies from zero to unity and is evaluated while taking into account the local error committed during the stress integration. This error cannot exceed a STOL tolerance. The first trial is conducted by adopting ΔT = 1. The procedure is controlled by pseudo-time T (0≤T≤1) at the end of each subincrement, ΣΔT = T = 1. The stress state is updated at the end of each sub increment starting from stress state upon the yield function (σ ⁰ = σ _int). The procedure described in this paper was implemented into ANLOG (Nogueira 2010) by Simão (2014).

3. Cylindrical cavity in elasto-plastic rock mass

This item describes the numerical simulation of a cylindrical cavity in a semi-infinite rock mass. The main goal of this simulation is to establish the plastic and elastic zones around the cavity and its stress distribution, highlighting the influence of the pressure acting internally on the cavity wall.

The problem is depicted in Figure 1 and consists of a cylindrical cavity of internal radius, ri, deeply conducted in a rock mass considered homogeneous, isotropic, and subjected to a isotropic initial stress state, σ₀. The stress-strain behavior of the rock mass is represented by using non-linear elastic-perfectly plastic, with associate plasticity, based on the Hoek-Brown criterion constitutive model as presented in the previous item. The elasto-plastic transition radius, R, defines a zone with plastic behavior (R-r_i) around the cavity and the external radius, r_e, defines the dominium of the problem.

Figure 1
Deep cylindrical cavity in rock mass

Sharan (2003, 2005) and Park and Kim (2006)presented an analytical solution for this problem considering the rock mass, under isotropic initial stress state, as elastic-brittle-plastic material which presents a sudden loss of strength after reach its maximal value. In this paper, this solution is adapted to consider an elastic perfectly plastic material which does not present this sudden loss of strength. Then, the radial and circumferential stress distribution on the plastic zone, r_i<r<R, is given by:

On the elastic zone, R < r < r_e, the radial and circumferential stress distribution is given by:

where, σ_R is a radial stress on the transition radius, R, defined as:

in which

Figure 1 presents the finite element mesh used in the numerical simulation; it is composed of 220 quadrilateral quadratic isoparametric elements (Q8) and 661 nodal points. The following constitutive parameters were adopted: E = 5.5GPa; ν = 0.25; σ_ci = 30MPa; m_b = 1.7; s = 0.0039; a=0.5 (which correspond to a GSI = 50 and m_i = 10, approximately).

The study uses a modified Newton-Raphson incremental iterativ e procedure with automatic increment of load (I_d = 10, miter = 20; toler = 0.1%; Δλ₀ = 0.01; Δλ_min = 10^-6; Δλ_max = 10^-2) and a Forward Euler stress integration (FTOL = 10^-5 and STOL = 10^-2).

Figure 2a illustrates the analytical and numerical (y = 0) results along the radial direction in terms of the radial (σ_r) and circumferential (σ_θ) stresses, considering an initial isotropic stress state (σ₀) of 30MPa and a null internal pressure (p_i). Table 1presents the normalized elasto-plastic transition radius (R/r_i) and stresses. As can be observed, numerical and analytical solutions agree strongly.

Figure 2
Stress distribution

Table 1
Elasto-plastic transition radius and stress (pi = 0MPa and ri = 5m)

Solution	R/ri	σr/σ0	σ0/σ0	σz/σ0
Analytical	2.833	0.526	1.474	0.500
Numerical	2.893	0.568	1.400	0.492

Figure 2b shows an elastic analytical solution provided by Kirsch (Poulos and Davis 1972) of a circular opening considering a horizontal stress coefficient of one. In this case, the circumferential stress decreases along the radial direction, while the radial normal stress increases in this direction. The elastic analysis overestimates the circumferential stress on the wall cavity. The elasto-plastic solution presents an abrupt decrease in the circumferential stress on the elasto-plastic transition zone.

Figure 3 presents a stress distribution around the cavity in terms of the isocurve of stress components. The highest vertical stress level is observed near the lateral wall of the cavity while the highest horizontal stress level is observed near its roof (Fig. 3)

Figure 3
Stress distributions – P_i = 0

Figure 4 shows the influence of the internal pressures acting on the internal wall cavity on the elasto-plastic transition radius. No plastic zone is observed from the internal pressure on the order of magnitude around half of the initial isotropic stress. The highest elasto-plastic transition radius is observed in an unsupported excavation, represented in this paper by a null internal pressure. In this case the radius depends on the property's material and the cavity radius.

Figure 4
Elasto-plastic transition radius - P_i varying

Figures 5a and 5b show the influence of the GSI on the elasto-plastic transition radius and on the normal radial and circumferential stresses at this point. As was expected, the elasto-plastic transition radius decreases as the GSI increases. Related to normal stresses, the radial stress decreases as the GSI increases, while the circumferential stress increases as the GSI increases.

Figure 5
Influence of GSI

5. Conclusion

The results presented in this paper demonstrated the importance of using the FEM and elastic-plastic constitutive model to simulate the opening cavity in rock masses. It was shown, for instance, that the stress distribution on the support structure changes significantly. By using the computer program ANLOG, it is possible to perform parametric studies for a wide variety of materials and geometrical configuration to improve the design of tunnel support.

6. Acknowledgments

The authors are grateful for the financial support received from CAPES, CNPq, and FAPEMIG. The first author is grateful for the grant received by CAPES. They also acknowledge Prof. John P. White for the editorial review of this text.

7. References

• CLAUSEN, J., DAMKILDE, L. An exact implementation of the Hoek-Brown criterion for elasto-plastic, finite element calculations. International Journal of Rock Mechanics & Mining Sciences, v. 45, n. 6, p. 831-847, 2008.
• CRISFIELD, M. Nonlinear finite element analysis of solids and structures. John Wiley and Sons, 1991.
• HOEK, E. Strength of rock and rock mass. ISRM New journal, v. 2, n. 2, p. 4-16, 1994.
• HOEK, E., WOOD, D., SHAH, S. A modified Hoek-Brown criterion for jointed rock mass. Proc. In: ROCK CHARACTERIZATION SYMPOSIUM, EUROCK'92, London, p.209-214, 1992.
• HOEK, E., BROWN, E.T. Underground excavation in rock London: Institution of Mining and Metallurgy, 1980.
• HOEK, E., BROWN, E.T. The Hoek-Brown failure criterion - a 1988 update. In: Rock Engineering for underground excavation, Proc. In: CANADIAN ROCK MECHANICS SYMPOSIUM, 15. Toronto: University of Toronto, p. 31-38, 1988.
• HOEK, E., CARRANZA-TORRES, C., CORKUM, B. Hoek-Brown failure criterion - 2002 edition, In: NORTH AMERICAN ROCK MECHANICS SYMPOSIUM, 15. Proc., TUNNELING ASSOCIATION OF CANADA CONFERENCE, 17. p.267-271, 2002.
• HOEK, E., KAISER, P.K, BAWDEN, W. F. Support of underground excavation in hard rock. Rotterdam: Balkema, 1995.
• NOGUEIRA, C. L. Análise não linear de escavações e aterros. Rio de Janeiro: Departamento de Engenharia Civil, Pontifícia Universidade Católica do Rio de Janeiro (PUC-Rio), 1998. 250 f. (Tese de Doutorado)
• OLIVEIRA, R. R. V. Análise elastoplástica via MEF de problemas em solos reforçados. Ouro Preto: Escola de Minas, Universidade Federal de Ouro Preto (UFOP), 2006. 143 f. (Dissertação de Mestrado em Engenharia Civil). OWEN, D.R.J. , HINTON, E. Finite elements in plasticity: theory and practice. Swansea, U.K.: Pineridge Press, 1980.
• POULOS, H. G., DAVIS, E. H. Elastic solution for soil and rock mechanics John Willey Ed., 1972.
• PARK, K-H., KIM, Y-J. Analytical solution for a circular opening in an elastic-brittle-plastic rock. International Journal of Rock Mechanics & Mining Sciences, v. 43, n.4, p.616-622, 2006.
• SHARAN, S. K. Elastic-brittle plastic of circular opening in Hoek brown media. International Journal of Rock Mechanics & Mining Sciences, v. 40, n. 6, p. 817-824, 2003.
• SHARAN, S.K. Exact and approximate solutions for displacements around circular openings in elastic-brittle-plastic Hoek-Brown rock. International Journal of Rock Mechanics & Mining Sciences, v. 42, p. 542-549, 2005.
• SLOAN, S. W., ABBO, A. J., SHENG, D. Refined explicit integration of elasto-plastic models with automatic error control. Engineering Computations, v. 18, n. 1/2, p. 121-154, 2001.
• TEIXEIRA, M.B., NOGUEIRA, C.L., OLIVEIRA FILHO, W.L. Numerical Simulation of Hillside Mine Waste Dump Construction. REM: Revista Escola de Minas, v. 65, n. 4, p. 553-559, 2012.
• SIMAO, J. T. Um modelo numérico para análise elastoplástica de maciços rochosos com base no critério de ruptura de Hoek-Brown. Ouro Preto: Escola de Minas, Universidade Federal de Ouro Preto (UFOP), 2014. 90 f. (Dissertação de Mestrado em Engenharia Mineral).
• WANG, S. & YIN, S. A closed-form solution for a spherical cavity in the elastic-brittle-plastic medium. Tunneling and underground space technology, n.26, p.236-241, 2011.

Publication Dates

• Publication in this collection
  Apr-Jun 2015

History

• Received
  19 Mar 2014
• Accepted
  18 Aug 2014