Applications of numerical limit analysis (NLA) to stability problems of rock and soil masses

ARTICLE IN PRESS

International Journal of Rock Mechanics & Mining Sciences 43 (2006) 408–425 www.elsevier.com/locate/ijrmms

Applications of numerical limit analysis (NLA) to stability problems of rock and soil masses A.F. Duranda, E.A. Vargas Jr.b,, L.E. Vazc a

Civil Engineering Laboratory, State University of Northern Rio de Janeiro (UENF), Av. Alberto Lamego, 2000, Campos dos Goytacazes, CEP 28013-602, Brazil b Department of Civil Engineering, Catholic University of Rio de Janeiro (PUC-Rio), Rua Marqueˆs de Sa˜o Vicente 225, Ga´vea, CEP 22453-060, Rio de Janeiro, Brazil c Department of Civil Engineering, Federal University of Rio de Janeiro (UFRJ), Cidade Universita´ria, CEP 21495-970, Rio de Janeiro, Brazil Accepted 31 July 2005 Available online 3 November 2005

Abstract This work presents formulations and results obtained with computer implementations of an alternative to the more standard techniques for the determination of the state of collapse of geotechnical structures in rock or soil masses. Examples of normally available and used techniques for those purposes are limit equilibrium based procedures and elasto-plastic finite elements. As an alternative to these techniques, the present paper describes Numerical Limit Analysis (NLA). The fundamentals for limit analysis, summarized in the so-called bound theorems, have been known for decades. Analytical solutions obtained with limit analysis are however limited in scope and are seldom used in the engineering practice. NLA on the other hand, by solving the limit analysis equations through numerical methods are general and applicable to a wide range of problems. The paper presents a discussion on available alternatives for the formulation of NLA specialized for the determination of collapse load factors of geotechnical structures in/on rock (fractured or not) and soil masses. Rock masses in particular are modelled as standard continua, Cosserat equivalent continua and true discontinua formed by discrete blocks. Finite elements are used for the solution of NLA equations of standard continua and Cosserat continua. The paper presents derivation of the pertinent equations, the numerical formulations used and details of their numerical implementation in computer programs. Attempt was made to validate all the implementations through existing analytical solutions. The obtained results permit to state that NLA is a promising and very often advantageous numerical technique to establish collapse states of geotechnical structures in rock and soil masses. r 2005 Elsevier Ltd. All rights reserved. Keywords: Limit analysis; Numerical limit analysis; Fractured rocks; Discrete blocks; Cosserat continua; Standard continua

1. Introduction In the design stage of a number of geotechnical problems such as the ones found in bearing capacity of foundations, retaining structures, slope stability and underground excavations, a primary objective consists in the determination of a collapse load, a maximum load the geotechnical system is able to support before it collapses. These loads are generally determined by limit equilibrium methods or elasto-plastic finite element analyses. In the present work, Corresponding author. Tel.: +55 21 3114 1188; fax: +55 21 3114 1195.

E-mail address: [email protected] (E.A. Vargas Jr.). 1365-1609/$ - see front matter r 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijrmms.2005.07.010

use is made of numerical limit analysis (NLA), an alternative, often very advantageous technique in relation to the ones described but in general seldom used in practice. The paper presents a theoretical background of limit analysis and possibilities of its numerical implementation. Main emphasis of the paper concerns applications of the technique for both continuum and discontinuum problems. In the latter case, more relevant to Rock Mechanics situations, the medium can be represented both by discrete blocks (true discontinua) and by Cosserat-type continua. Initially, the paper formulates the general limit analysis problem. Subsequently, it describes its specialization and numerical implementation for analysis of

ARTICLE IN PRESS A.F. Durand et al. / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 408–425

standard continuua, Cosserat continua and true discontinuua (discrete blocks) problems. Application and validation examples are presented and discussed. Finally, a general discussion on the possibilities of the use of NLA in practical problems is presented. 2. Fundamentals of limit analysis The interest in the study of collapse of engineering structures began with Coulomb in 1773 [1]. Coulomb developed plasticity concepts applied to soils as well as the concept of a plastic limit as applied to retaining structures. Later, in 1857, Rankine introduced the concept of slip lines in order to interpret the plastic equilibrium of soils. Based in such studies, the concept of limit analysis evolved, more or less heuristically, until approximately 1950 in different areas of engineering. In 1952, Drucker and Prager, in a study of plastic materials obeying Mohr–Coulomb’s yield criterion [2,3], defined bounds (upper and lower) for the collapse load. The collapse load in stability problems can be obtained in an independent way through the application of the theorems of the limit analysis (upper bound or lower bound). When the collapse load is obtained from an established statically admissible stress field, it is considered to be a lower bound to the true collapse load. In the same way, if the collapse load is obtained by means of an established kinematically compatible failure mechanism, it is considered to be an upper bound to the true collapse load. In practice, the establishment of kinematically admissible failure mechanisms is in general easier when compared with the establishment of statically admissible stress fields. In this way, various analytical solutions are available today using various forms of simplified collapse mechanisms [4]. 2.1. Limit analysis theorems The application of the limit analysis theorems is possible for solid materials that present the following properties [1]: 1. The plastic behavior of the material is perfect or ideally plastic, i.e., the yield surface is fixed in the stress space. 2. The yield surface is convex and the rates of plastic deformation are obtained from a yield function through an associated flow rule. 3. Changes in the geometry of the body are considered insignificant when the loads reach the limit load state. The principle of virtual work can therefore be applied. Central in the theory of limit analysis is the determination of a collapse load factor, roughly speaking an analogous measure to the factor of safety in limit equilibrium. In the lower bound formulation, the collapse load factor, or (or simply) the load factor can be defined as a multiplying factor (a scalar) by which the external loads have to be multiplied in order that the structural system reaches

409

collapse. Similarly, in the upper bound formulation, the load factor is defined as a multiplying factor (also a scalar) by which the external work (work performed by the external loads) has to be multiplied in order that the structural system reaches collapse. Formally, the lower bound (or static) and the upper bound (kinematic) theorems can be stated as in the following: Lower bound theorem (static theorem). A load factor ls that corresponds to a statically admissible stress field, one that satisfies (a) equilibrium equations in the domain, (b) equilibrium equations on the boundary and (c) nowhere in the domain the yield function is violated, will not exceed the true load factor of the structure. Upper bound theorem (kinematic theorem). The kinematic load factor lk as determined by equating the rate of external work with the rate of internal dissipation of energy along a kinematically admissible velocity field (_u), one that satisfies (a) the velocity boundary condition and (b) the compatibility relations between strain rates and velocity, is not less than the true collapse load factor. According to the principles of continuum mechanics, the statements contained in both upper and lower bound theorems can be stated mathematically by the following equations: Given f t

in O ðbody forcesÞ; on Gt ðboundary forces; tractionsÞ;

Find l, r, u_ , e_ , e, c_ such that the following conditions are satisfied: Static equilibrium: rT r ¼ lf in O; rg ¼ lt on Gt ;

(1)

Yield criterion: f ðrÞp0

in O,

(2)

Kinematic consistency: e_ ¼ ru_ in O; u_ ¼ 0 on Gu ;

(3)

Flow rule: e_ p ¼ g_

qf qs

g_ ¼ 0 g_ X0

if f ðrÞo0; if f ðrÞ ¼ 0;

(4)

where f is the vector of body forces, t the boundary forces, tractions, g the unit vector normal to surface Gt , r the stress field, u_ the velocity (displacement rate) field, e_ p the plastic strain rates field, l the collapse load factor and g_ the plastification factor. The complete solution of the problem of establishing collapse of a system considering both statically and kinematically admissible fields must use Eqs. (1)–(4). The problem however can be solved by considering either a statically admissible field or a kinematically admissible field. When the collapse load of the system is obtained

ARTICLE IN PRESS 410

A.F. Durand et al. / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 408–425

through conditions represented by Eqs. (1) and (2), the formulation is denominated static (lower bound). On the other hand if the collapse load of the system is obtained through the use of Eqs. (3) and (4), the formulation is denominated kinematic (or upper bound). Alternative formulations, denominated mixed formulations, are the ones that consider at the same time all Eqs. (1)–(4). However, obtaining solutions to all equations is in general very cumbersome and difficult. 3. Numerical limit analysis of geotechnical problems Lower ðls Þ and upper ðlk Þ bound estimates of the collapse load factor can be obtained by both analytical and computational methods. Finding analytical solutions for practical engineering soil and rock mechanics problems is, however, a difficult task. Chen [1], Chen and Liu [4] and Finn [5] have presented purely analytical applications of limit analysis to some practical geotechnical problems. In the case of lower bound solutions, statically admissible stress fields have to be assumed, whereas in the case of upper bound solutions, a kinematically admissible velocity field must be assumed. In the latter case, this can be done by establishing a failure mechanism. In the seventies, numerical methods, the finite element method in particular, together with mathematical programming techniques, were brought into use to solve limit analysis problems. The simultaneous solution of Eqs. (1)–(4), even numerically, is no easy task due to the complexity in finding simultaneously statically admissible stress fields and kinematically admissible velocity fields. Formally, for any of the abovementioned formulations, the set of corresponding equations constitute a mathematical programming (optimization) constrained problem. This optimization problem can be either linear or non-linear. In the particular case of limit analysis this distinction will depend on the form of the yield condition (2) used. A number of alternative formulations and numerical implementations are available in the technical literature [3,6–20]. Basically, the strategies described in abovementioned papers differ in the way the equilibrium and compatibility equations are considered and assumptions made about the stress and velocity fields. In the sequence, based on the characteristics of the existing formulations, an attempt is made on classifying them. Initially, the strategies are classified according to the method used for discretizing the equilibrium and compatibility equations. In sequence the strategies are classified according to the assumptions made regarding stress and velocity fields.

mined directions, the corresponding strategies are denominated strong formulations of limit analysis. In contrast, weak formulations can be obtained via the use of the principle of virtual work. 3.2. Classification according to assumptions made about admissibility of stress and velocity fields When the solution of the problem posed guarantees that, in all points of the domain, the stress field is statically admissible, the formulation is denominated as lower bound. Alternatively, when the solution of the problem guarantees that the velocity field is kinematically admissible, the formulation is denominated as upper bound. When, on the other hand, the formulations do not guarantee these conditions in the domain, they are designated as approximate formulations. In this case, it is expected that the solution converges to the true solution when the mesh is refined but the convergence does not respect bounds. Besides the two broad classes of formulations described above, one also distinguishes the so-called mixed formulations, so denominated because they consider as unknowns both the stress and the kinematic fields satisfying partially or totally, respectively, the static and kinematic conditions. It is, therefore, possible to classify some of the published formulations for the problem of limit analysis. The formulations presented by [3,7,11] could be classified as strong formulations with statically admissible stress fields. The formulations presented by [7,12,19] could be classified as strong formulations with kinematically admissible velocity fields. Examples of weak formulations with statically approximated stress fields are described by [6,20]. 3.3. Formulations used in the present work The present work does not intend to present a discussion on pros and cons of the different formulations. Instead, two particular formulations are described and used throughout the paper, emphasis being on their applicability to stability problems in Rock and Soil Mechanics. A strong formulation with statically admissible stress field is applied to systems of discrete rigid blocks simulating a fractured (blocky) rock mass. A weak formulation with a statically approximate stress fields is applied to the solution of rock and soil mechanics problems where the geological medium is modelled by both conventional and Cosserat continua [17,21]. 4. Numerical formulation of limit analysis

3.1. Classification according to the method used for discretizing the equilibrium and compatibility equations When the discretization of the equilibrium or compatibility equations is obtained through the explicit application of the differential equations of equilibrium or compatibility (Eqs. (1) and (3)), in points and predeter-

4.1. Rock mass modelled as a discontinuum (discrete blocks) This section describes the use of a strong formulation associated with a statically admissible stress field to perform limit analysis of a discontinuum simulating a

ARTICLE IN PRESS A.F. Durand et al. / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 408–425

The assembly of the described matrices and vectors are described next according to the diagrams shown in Fig. 1. A more complete derivation of the matrices presented below can be found in [23].

y Block i xi

Block i C

xc

fn ci

α

fs cj

xj

Block j

fs yi

yc

yj

x

(a)

ci

4.1.1. Assembly of matrix H (equilibrium of blocks) The assembly of matrix H for the whole structure is obtained by incorporation of each individual Hi matrix for each block. With reference to Fig. 1, the assembly of each Hi for a block i is given by

C

fn cj

α

Block j

(b)

Fig. 1. Geometry of block contacts: (a) global reference system; (b) normal and tangential forces in contact C. Angle a is the angle between the block surface and the horizontal.

fractured (blocky type) rock mass. The medium is composed of a system of discrete rigid blocks. The 2

cos a1

6  sin a 1 6 6 6 6 6 L¼6 6 6 6 6 4

(8)

Hi ¼ AL,

where L is the local to global reference system transformation matrix, and A the equilibrium matrix considering forces and moments at the centroid of each block i, and L and A are given by 3

sin a1 cos a1

0 cos a2

sin a2

 sin a2

cos a2 :::

0

cos ak  sin ak

2 6 A¼4

1

0

0 1 ðyc1  yi Þ ðxc1  xi Þ

l

7 7 7 7 7 7 , 7 7 7 7 sin ak 7 5 cos ak 2k
2k

1

0

:::

1

0 ðyc2  yi Þ

1 ðxc2  xi Þ

::: :::

0 ðyck  yi Þ

procedure adopted follows the work of Livesley [22] who stated that the collapse behavior of a structure formed by discrete, rigid blocks (Fig. 1) can be formulated by equilibrium statements for each block and by strength constraints at each block interface. In the context of the present work, a linear Coulomb type criterion is used for that matter. The mathematical programming (optimization) problem can in this case be stated as: Maximize

411

0

(9)

3

7 1 5 ðxck  xi Þ

,

(10)

3
2k

where k is the number of contact points for each block i. The global matrix H for the whole assembly of blocks can therefore be assembled: 2 3 H1 0 6 7 H2 6 7 H¼6 . (11) 7 ::: 4 5 0 HNB 3NB
2NTC

(5)

subject to Hrc ¼ lq equilibrium condition,

(6)

Rrc prf

(7)

yield condition,

where NB is the total number of blocks of the system, NC the total number of contacts of all blocks, H the equilibrium matrix which relates the forces at each block contact with the forces and moments and each block centroid (3NB
2NC), R the matrix of coefficients of the failure (yield ) criterion (2NC
2NC), q the vector of forces and moments at each block centroid (3NB), rc the vector of forces at each contact (2NC) and rf the vector of failure criterion constants (2NC).

4.1.2. Assembly of vectors q, rc and rui Vectors qi , rc and uui are given by 9 8 c 9 8 fss1 > fchc1 > > > > > > > > > > c> > > > > > > fn1 > fchc1 > > > > > > > > > > > > > > > > > c c > > > > fs fch 9 8 > > > 2 > 2> > > > > > > > F > > > c= > > > ix = c= < < < fn2 fch2 c f ri ¼ ; qi F iy ; ri ¼ , : > : > > > > > > > > ; :M > > > > > > > i > > > : > > : > > > > > > > > > > > > > > > c > > c> > > > > > > fs fch > > > k > k> > > > > > > > ; : fnc ; : fchc > k k 2k

2k

(12)

ARTICLE IN PRESS A.F. Durand et al. / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 408–425

412

where fsc1 is the shear force of contact 1 of block i, fnc1 the normal force of contact 1 of block i, F ix the resultant force in the x direction at the centroid of block i, F iy the resultant force in the y direction at the centroid of block i, M i the resultant moment at the centroid of block i, and fchc1 the equivalent force due to cohesion at contact 1. 4.1.3. Assembly of matrix R (strength restrictions) Using a Coulomb type criterion, matrix Ri for each block i containing k contacts can be represented: 3 2 1 m 0 7 6 1 m 7 6 7 6 7 :: Ri ¼ 6 , (13) 7 6 7 6 1 m 5 4 0 1 m

x3 σ33 +

σ31 +

∂σ31 dx3 ∂x3

∂σ σ32 + 32 dx3 σ11 ∂x3

σ23 +

σ12

σ21 σ13 +

σ22

∂σ33 dx3 ∂x3

∂σ13 dx ∂x1 1

σ23

∂σ23 dx ∂x2 2 σ22 +

σ13

σ12 +∂σ12 dx1 ∂x1

σ21 +

σ21 dx σx2 2 x2

σ31 σ11 +

∂σ22 dx σx2 2

∂σ11 σ dx1 32 ∂x1

2k
2k

σ33

where

x1

m ¼ tan f,

(14)

f being the friction angle at the block contacts. The global constraint matrix R for the whole assembly of blocks can therefore be assembled: 2

R1

6 6 R¼6 4

0

7 7 7 5

R2 :: 0

3

RNB

.

(15)

2NTC
2NTC

The mathematical programming problem generated by Eqs. (5)–(7) constitutes a linear problem as both the equilibrium and the yield conditions at the contacts and the objective function are linear. The variables of the problem are vectors rc and scalar l. Details of numerical implementation and solution of the problem are described in Section 4.4. 4.2. Numerical formulation for limit analysis of conventional continua A weak formulation with statically approximate stress fields is used in this work to obtain numerical limit analysis solutions to soil/rock stability problems where the geological materials involved can be represented as conventional continua. The derivation details of the formulation are presented in sequence. Some advantages of this formulation when compared to others [3,6–20] are: the numerical implementation is simple, the estimates for the collapse load factor are in a high degree mesh independent and they converge to the exact solution when the mesh is refined [16,20].

Fig. 2. Stress field for an infinitesimal element in a 3D conventional continuum.

In reference to Fig. 2, one obtains the static equilibrium equations: 9 8 s11 > 3> 2 > > > q q q > > > > > 0 0 0 > s 22 > 7 6 qx1 > > qx qx > > 2 3 7> > 6 > > > 7> 6 < 7 s33 = 6 q q q 7 6 0 0 0 7> s > 6 qx2 qx1 qx3 12 > 7> 6 > > > 7> 6 > > > q q q 5> 4 > > s > > 23 0 0 0 > > > > qx3 qx2 qx1 > > ; : s31 8 9 8 9 q1 > > 0 > > > = < > = > < > þ q2 ¼ 0 . ð16Þ > > > > > ; : > ; > :q > 0 3 Considering 2 q q 0 0 0 6 qx1 qx2 6 6 6 q q q 0 =¼6 6 0 qx qx qx 2 1 3 6 6 q q 4 0 0 0 qx3 qx2 9 8 s11 > > > > > > > > > > > > s 22 8 9 > > > > > > q1 > > > > > > = < s33 > = < > and q¯ ¼ q2 , r¼ > > > s12 > > > > > > ; :q > > > > > 3 > > > > > > > s23 > > > > > ; : s31

3 q qx3 7 7 7 7 0 7 7, 7 7 q 5 qx1

ð17Þ

ARTICLE IN PRESS A.F. Durand et al. / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 408–425

where r is the stress field vector, q the applied body forces vector, = the differential operator matrix. Eq. (16) may be written as

Simplifying Z  Z Z T T T ðBu Hs Þ dO ¼ l Hu f dO þ Hu t dGt O

=r þ q ¼ 0.

(18)

O

Gr^ ¼ lp,

e_ ¼ =T u_ ,

where Z G ¼ ðBTu Hs Þ dO;

where 8 9 e_11 > > > > > > > > > > _ e > > 22 > > > > > = < e_33 > ; e_ ¼ > e_12 > > > > > > > > > > > > e_23 > > > > ; : e_ >

8 9 > = < u_ 1 > u_ ¼ u_ 2 . > ; : u_ >

(20)

3

O

Gt

In Eq. (21), the term on the left-hand side represents the virtual work performed by the stresses in the domain on virtual strain rates. The first term on the right-hand side represents the virtual work of the body forces distributed in the domain and the second term represents the virtual work of the tractions applied at the surface. Particular forms of (21) can be generated when the collapse factor multiplies only the body forces and, therefore, only these forces will be maximized (volume loads are active and the external surface loads are passive loads) and when the collapse factor multiplies only the surface forces, therefore, only these forces will be maximized (external surface loads are active and the volume loads are passive loads). Considering the numerical discretization for the velocity and stress fields, we have u_ ¼ Hu u^

(22)

where r is the finite element stress field, u_ the finite element velocity field, r^ the nodal stress vector, u^ the nodal velocity vector and Hs , Hu are the stress field and velocity field interpolation matrices, respectively. The vector for the velocity deformation rates ð_eÞ is defined by e_ ¼ r u u_ ¼ r u Hu u^ ¼ Bu u^ ,

(23)

where the matrix Bu is obtained from Bu ¼ r u Hu .

(24)

The substitution of Eqs. (22) and (23) into (21) yields Z  Z Z T T T T T d^u Bu Hs r^ dO ¼ ld^u Hu f dO þ Hu t dGt . O

Z p¼ O

HTu f

Z dO þ Gt

HTu t dGt

 . (27)

In the above equations, e_ is the deformation rates vector and u_ is the velocities vector. The discrete equilibrium equations used in the weak mixed formulation for limit analysis can be derived from the principle of virtual velocities field:
Z  Z Z T T T ðd_e rÞ dO ¼ l ðd_u fÞ dO þ ðd_u tÞ dGt . (21)

^ r ¼ Hs r;

(26)

O

31

O

(25)

Gt

which in compact form may be written as

The deformation rates are related to the velocity field by [24] (19)

413

O

Gt

G is the equilibrium matrix and p is the vector of applied nodal forces. Using the discrete form of the equilibrium equations derived above and the failure criterion defined for each element of the finite element mesh, the collapse load factor problem can be formulated as a lower bound limit analysis problem. Once the stress field is not an admissible one, the problem is not strictly a lower bound limit analysis problem and there is no guarantee of convergence to a lower bound solution. On the other hand, as the mesh is refined, the stress field can approach an admissible one and the estimate converges to the true collapse load factor. The compact form of the optimization problem generated is in this case: Maximize

l

(28)

subject to Gr^ ¼ lp ^ f ðrÞp0

equilibrium condition, failure criterion for each element;

(29) (30)

^ is the function which represents the failure where f ðrÞ criterion for the material. ^ 4.2.1. Failure criteria for a conventional continuum f ðrÞ Failure criteria for geological materials are a function of the stress vector. In the present work the following failure criteria were used and implemented: Mohr–Coulomb (31), Hoek and Brown (32) and Drucker–Prager (33) failure criteria: ^ M2C : s1 ½1  sinðfÞ  s3 ½1 þ sinðfÞ  2c cosðfÞp0, f ðrÞ (31) ^ H2B : ½s1  s3 2  ½ss2c þ msc s3 p0, f ðrÞ pffiffiffiffiffiffiffiffi ^ D2P : J 2D  aJ 1  kp0, f ðrÞ

(32) (33)

where s1 is the maximum principal stress, s3 the minimum principal stress, J 2D the second invariant of the deviatoric stress tensor, J 1 the first invariant of the stress tensor, sc the uniaxial compression strength of the material, c and f are the cohesion and friction angle of the material for the Mohr–Coulomb yield function, respectively, s and m are the material parameters of the Hoek and Brown failure criterion and a and k the material parameters of the Drucker–Prager failure criterion.

ARTICLE IN PRESS A.F. Durand et al. / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 408–425

414

σ22 +

µ32 +

x2

∂µ32 ∂x2

∂σ22 ∂x2

∂σ21 σ21 + dx ∂x2 2

dx2

σ12 + F2

µ31 σ11

dx2

9 8 9 8 F1 > >0> > > > = < > = > < F þ ¼ 0 . 2 > > > > > ; : > ; > :M > 0 3

dx2

∂σ12 ∂x1 µ31 +

F1

dx1 ∂µ31 dx1 ∂x1 ∂σ11 σ11 + dx1 ∂x1

M3 σ12

dx1 σ21 µ32 σ22

x3

x1

Fig. 3. Stress field of a microstructure in a Cosserat representation.

For the 3D cases analysed, only the Drucker–Prager criterion (33) was used and its form in 3D is the same as in 2D. Eqs. (28)–(30) constitute a non-linear mathematical programming problem as the yield condition (failure) is a non-linear relationship. Details of numerical implementation and solution of the problem are described in the next section. 4.3. A numerical formulation for limit analysis of Cosserat continua The relevance of Cosserat continua [25] for the modelling of the mechanical behaviour of (fractured) rock masses has been stressed by a number of authors (see, e.g., Chappel [26]). A particular and relevant feature of Cosserat continua is its ability to incorporate extra degrees of freedom when compared to the standard continuum. These extra degrees of freedom appear due to the fact that, in Cosserat continua, a microstructure having dimensions is taken into account in each point in the domain. These microstructures are able to transmit moments. Fig. 3 shows the stress state existing at the level of a microstructural unit in a Cosserat continuum in 2D. There one is able to notice that besides the usual stresses sij , the so-called couple stresses (distributed areal moments mij ) are introduced [21,27]. With reference to Fig. 3, the static equilibrium equations are 9 8 s11 > 3> 2 > > q q > > > > > 0 0 0 0 > > > s 22 7 6 qx > > qy > > > 7> 6 > > > 7> 6 < s 12 = q q 7 6 7 6 0 0 0 0 7> s > 6 qy qx 21 > 7> 6 > > > 7> 6 > > > 4 q q 5> > > m > 31 > 0 0 þ1 1 > > > > qx qy > > ; : m32

ð34Þ

In analogy to (18) one may write

r c rc þ qc ¼ 0,

(35)

where 2

q 6 qx 6 6 6 c r ¼6 0 6 6 4 0

0

0

q qy

0

q qy

q qx

0

0

0

þ1 1

q qx

0

3

7 7 7 0 7 7 7 7 q 5 qy

(36)

and sc is the Cosserat stress field vector and qc the applied Cosserat body forces vector. The superscript c stands for Cosserat continua. Another particular and relevant feature of Cosserat continua is the fact that s126¼s21, meaning that the stress tensor is asymmetric. This happens because of the intrinsic microstructure dimensions and the couple moments included in the formulation. The asymmetry in the stress tensor generated has considerable influence on the failure mechanisms which arise in the medium. The corresponding Cosserat deformations are [28] 9 8 g_ 11 > > > > > > > > > > > _ 22 > g > > > > > > > > > = < g_ 12 > c_ c ¼ > > > > > g_ 21 > > > > > > > > > _ k > > 31 > > > > > > ; : k_ 32 3 2 q 0 0 7 6 qx 7 6 7 6 q 6 0 0 7 7 6 qy 7 6 7 6 78 9 6 q 6 0 1 7 > > u_ > 7> 6 qx 7< = 6 7 v_ , ¼6 ð37Þ 7> > 6 q > 7> 6 ; : 0 þ1 7 o 6 qy 7 _ 6 7 6 6 q 7 7 6 0 0 6 qx 7 7 6 7 6 4 q 5 0 0 qy where u_ is the microstructure velocities in x direction, v_ the _ the microstrucmicrostructure velocities in y direction, o ture rotation velocities, k_ 31 the curvature velocities in ‘‘x’’ direction and k_ 32 the curvature velocities in ‘‘y’’ direction. The appearance of curvatures arises from the couple stresses and this means that the medium is able to ‘bend’

ARTICLE IN PRESS A.F. Durand et al. / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 408–425

as a structural element such as a beam or a plate. Considering an approximate stress field as well as an approximate velocity field within a finite element, one has the following relations: c

c

c

c

r ¼ Hs r^ ,

(38)

u_ ¼ Hu u^ ,

(39)

where the vectors representing the stress and the deformation rate fields within the finite element for the Cosserat medium in a state of plane strain are, respectively, 9 8 s11 > > > > > > > > > > > s 22 > > > > > > > > > > = < s12 > c , r ¼ > > > > > s21 > > > > > > > > m13 > > > > > > > > > ; : m23 9 8 g_ 11 > > > > > > > > > > > _ g > > 22 > > > > > > ( ) > > = < _ g12 > g_ c ¼ ru u_ c ¼ ru Hu u^ c , g_ ¼ ð40Þ ¼ > > _ g k_ > > 21 > > > > > > > > > > > k_ 13 > > > > > > > ; : k_ 23 where rc is the Cosserat stress field within the finite element, r^ c the vector of Cosserat nodal stresses, Hs the matrix of the stress interpolation functions,cu_ c the Cosserat _ velocity field within the finite element, u the vector of Cosserat nodal velocities, Hu the matrix of the velocities interpolation functions and gcij the vector of the Cosserat deformation rates. The implemented numerical formulation makes use of the weighted residual method applied to the equilibrium equations: Z d_ucT ðrs rc þ fÞ dO ¼ 0 (41) O

Separating terms Z Z ðd_ucT rs rc Þ dO þ ðd_ucT fÞ dO ¼ 0. O

(42)

Z O

BTs Hs dO r^ c ¼ 

Z O

HTu f dO

415

(44)

or, simply Cr^ c ¼ q,

(45)

where BTs ¼ rTs Hu , Z C¼ BT Hs ðxÞ dO,

(46) (47)

O

Z q¼

Hu ðxÞT f dO,

(48)

O

C is the equilibrium matrix and q is the vector of applied nodal forces. Once again, good estimates are obtained when the mesh is refined relaxing the constraints on the stress field. The compact form of the optimization problem generated is in this case: Maximize

l

(49)

equilibrium condition,

(50)

subject to Cr^ c ¼ lq f ðr^ c Þp0

failure criterion for each element:

(51)

The load factor l multiplies the active external loads. In (50), all the nodal, surface and volume external loads are active. 4.3.1. Failure criteria for Cosserat continuum Here, the generalized Cosserat continuum presented in the previous Section is specialized for the case of a rock mass containing two families of discontinuities [29,30], as shown in Fig. 4. Families a and b divide the medium into parallelograms with sides ‘a and ‘b, which will, thereafter, constitute the basic structure (microstructure) of Cosserat continuum. This basic structure allows for three failure mechanisms, sliding and separation at the fractures and rotations of the basic structures formed as shown in Fig. 5. The model includes two families of fractures and three mechanisms of failure for each family (sliding, separation and toppling) [29,30]. Failure by sliding must satisfy the Coulomb criterion (52). Failure by separation is limited by

O

Considering Eqs. (38)–(40): Z Z cT c cT T d^u Hu rs Hs r^ dO ¼ d^u HTu f dO O

(43)

O

and eliminating d^ucT on both sides of (43): Z Z c T ðHu rs ÞHs dO r^ ¼  HTu f dO, O

Z O

O

ðrTs Hu ÞT Hs dO r^ c ¼ 

Z O

HTu f dO,

Fig. 4. Rigid blocks separated by two families of discontinuities having inclinations ya and yb from the positive x axis.

ARTICLE IN PRESS A.F. Durand et al. / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 408–425

416

relationship. Details of numerical implementation and solution of the problem are described in the next section. 4.4. Numerical implementation and solution of the mathematical programming problems

Fig. 5. Failure mechanisms: sliding, separation and rotation [29].

l

l

σ ~ o α

Cm τα mα σα

Fs

o mα

o

Fn

mα Fn

Fig. 6. Geometry of the toppling condition [31].

the tensile strength of the fractures according to Eq. (53). Failure by toppling (54) is illustrated by the rigid body diagram of Fig. 6. This figure shows that the moments at the faces of the blocks are the ones responsible for an eventual toppling failure of the block in relation to edge ‘‘o’’. f 1 : jss j þ sn tan f  cs p0,

(52)

f 2 : sn  T 0 p0,

(53)

‘ f 3 : jmj þ sn  cm p0, 2

(54)

where ss is the shear stress, sn the normal stress, T0 the tensile strength of fracture, f the friction angle at the face of block, cs the (stress) cohesion at the fracture plane, cm the moment cohesion of the block, ‘ the length of a face of the block and m the moment per unit area of the block face. One may observe that a new property of the system, designated as moment cohesion (cm), appears. In principle, it can be regarded roughly as an analogous measure, in terms of moments, to the standard (stress) cohesion cs. Its physical meaning, however, is not yet fully understood. It may be related to different characteristics of the rock mass such as the existence of filling in the discontinuities, the existence of intact rock bridges or even some degree of interlocking of the rock blocks. There are, to date, no testing procedures designed especifically for its determination although backanalyses of actual failure cases may serve that purpose. The analysis presented in Section 5.3 is an example of that. Eqs. (49)–(51) constitute a non-linear mathematical programming problem because although the yield conditions at the fractures are linear (52)–(54), the yield condition (failure) for the intact rock is a non-linear

This section describes details of the numerical implementation through the development of a number of computer programs for the solution of the equations described in Sections 4.1, 4.2 and 4.3. For the assembly of Eqs. (5)–(7), related to limit analysis of discrete blocks, program LIMAG-BLOK was developed [32]. In this program, the coordinates of the gravity centres of the blocks and the contact points between blocks (two contact points between two blocks) are determined automatically by program rigid block model (RBM) [33]. Three plane equilibrium equations relating the contact forces acting on the boundary of the blocks to the applied forces acting on the centroid of each block are generated by the program. The yield conditions at each contact between blocks are also defined, relating the contact forces to the strength parameters of the discontinuities. The complete set of equations defines a linear programming problem, which is finally solved by the commercial program LINDO [34] to obtain the collapse load of the system. The solution of the pertinent equations provides a value for the load factor, the forces at the contacts and velocities in the domain. The velocities are not determined directly but indirectly through the dual statement of the mathematical programming problem [35]. For the assembly of Eqs. (28)–(30), related to limit analysis of geotechnical problems of a standard continuum, program limit analysis for geotechnics (LIMAG) was implemented. LIMAG is associated with a preprocessor for automatic generation of meshes (in 2 and 3D). The interpolation functions considered for the stress and velocity fields are represented by matrices Hr, and Hu (Eq. (22)). Different interpolation functions were implemented but the most often used ones were a constant distribution for stresses and a bilinear distribution for velocity fields. The reason for using these types of element formulation was that it provided the most consistent results both in terms of load factors and stress distribution. It is not known at this moment the theoretical reasons for it to occur. It is, in the opinion of the authors, a point to be investigated in more detail. On the other hand, it was observed that when using bilinear stress and velocity fields the finally obtained distribution of the stress field was sometimes inconsistent although the corresponding load factor was close to the one obtained with constant stress and bilinear velocity fields. In LIMAG, the yield criteria of Mohr–Coulomb, Hoek and Brown and Drucker and Prager were implemented (in 2D and 3D) according to the forms shown in Eqs. (31), (32) and (33), respectively. For the solution of the non-linear optimization problem, LIMAG provides the users the option of using programs LINGO [34], MINOS [36] and LANCELOT [37]. These

ARTICLE IN PRESS A.F. Durand et al. / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 408–425

programs are well known commercial packages for the solution of general optimization, linear and non-linear problems. For the assembly of Eqs. (49)–(51) related to limit analysis of a Cosserat continuum, the program LIMAG/COS was implemented. This program considers the same finite elements existing in LIMAG with the addition of extra interpolation fields for couple stresses and rotations. In all cases described above related to continuum problems, the solution of the pertinent equations provides a value for the load factor, the stresses and velocities in the domain. In the case of Cosserat continua, couple stresses and rotations are included. The velocities are not determined directly but indirectly by means of the dual statement of the mathematical programming problem [35]. 5. Validation and illustrative examples 5.1. Load capacity problem in shallow foundations This example deals with the determination of the collapse load of shallow foundations on rock or soil. These materials are, for that matter, considered as a conventional continuum in the numerical model (see Section 4.2). The numerical limit analysis (NLA) problem is represented in this case by Eqs. (28)–(30). In the formulation, the action of the shallow foundation on rock or soil is included in vector p. The numerical solution of the optimization problem determines the nodal stress vector in equilibrium with the amplified active boundary forces and satisfying the failure criterion. Two finite element models were used for this numerical limit analysis, namely a 2D model and a 3D model. The 2D model uses 4-noded finite elements while the 3D model uses 8-noded finite elements for describing the virtual velocity field. In both cases, a constant stress distribution within the element is considered. Based on the mesh discretization, the equilibrium equations, which correspond to the equality constraints in the optimization problem, are generated automatically. For the generation of the inequality constraints, Mohr–Coulomb, Hoek–Brown and Drucker– Prager failure criteria were applied depending on the materials involved in the analysis. 5.1.1. Load capacity of strip footings using the Mohr– Coulomb failure criterion Bowles (1977) [38] presents an analytical solution for this problem considering that the material satisfies the Mohr– Coulomb failure criterion. According to Bowles, the collapse load factor for a shallow strip foundation can be determined analytically from the following equations: qu ¼ cN c

(55)

    f p tan f 2 N c ¼ cot f e tan 45 þ 1 , 2

(56)

417

where qu is the load capacity. The parameters used in the validation problem were cohesion c ¼ 0:5 MPa and angle of friction f ¼ 151. Substituting these values into Eqs. (56) and (55), the collapse load qu ¼ 5:488 MPa for the strip shallow foundation is obtained. Fig. 7 depicts the deformed finite element mesh, the vertical stress and the velocity vector distributions in the imminence of collapse, given by the code LIMAG. The collapse load factor obtained for this problem with NLA was l ¼ 5:49, which corresponds to a collapse load of 5.49 MPa. This value is quite close to the result obtained with the analytical solution. In the sequence, a 3D load capacity problem is considered. It involves the determination of the load factor of a rectangular footing (with a smooth base) as a function of its geometry. The geometry of the footing is expressed by the relationship B/L where B is the smallest and L the largest dimension of the rectangle. Due to the numerical instabilities encountered with the Mohr–Coulomb strength criterion when solving 3D problems, Drucker and Prager strength criterion was used instead. Parameters a and k (33) adopted in the analysis were defined by fitting the Drucker and Prager criterion to the Mohr–Coulomb strength criterion, leading to following expressions [39]: 2 sinðfÞ a ¼ pffiffiffi ; ½3 þ sinðfÞ 3

6c cosðfÞ k ¼ pffiffiffi . 3 ½3 þ sinðfÞ

(57)

In Fig. 8 the collapse mechanism and the vertical stress field given by the code LIMAG for a shallow rectangular foundation with a relation B=L ¼ 1=3 is shown. The soil parameters used in the analysis were cohesion ¼ 1:0 MPa and friction angle ¼ 351. Fig. 9 shows results of a parametric study for the determination of the load factor of a rectangular footing having a smooth base as a function of its geometry. Along with the results obtained with code LIMAG for NLA, results using a 3D pseudo-analytical solution proposed by de Beer in 1970 as presented by Das [40] (Eqs. (58) and (60)) are also plotted. The figure also includes results obtained using an analytical solution proposed by Das [40] for square footing (B=L ¼ 1:0) (Eq. (60)). It may be observed in this figure that the solution for the collapse load factor via NLA for the footing with B=L ¼ 5 is very similar to the solution for the strip footing. qu ¼ cN c F cs , F cs ¼ 1 þ

B Nq , L Nc

qu ¼ 1:3cN c ,

(58) (59) (60)

where qu is the load capacity of footing and Nc and Nq are the load capacity factors. The figure shows that the load factors obtained with NLA are smaller than the ones predicted by de Beer’s solution. They compare well, however, with Das’s solution

ARTICLE IN PRESS 418

A.F. Durand et al. / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 408–425

Fig. 7. Limit analysis of a shallow strip footing. Solution given by the optimizer MINOS using a mesh of 1680 elements and Mohr–Coulomb strength criterion (c ¼ 0.50 MPa, f ¼ 151). Collapse load factor l ¼ 5.49. (a) Deformed finite element mesh, (b) vertical stress distribution (MPa) in the imminence of collapse and (c) velocity vector distribution in the imminence of collapse.

for a square footing. Furthermore, NLA solution for a square footing with B=L ¼ 1=5 matches the solution for a strip footing (2D) as could be expected. Fig. 10 illustrates the failure mechanism and the vertical stress distribution for the case B=L ¼ 1=5, one which in practice corresponds to a strip footing problem.

Fig. 8. Rectangular rigid footing (smooth base). Deformed mesh and vertical stress distribution (MPa) obtained using optimizer LANCELOT; inferior circle approximation parameters (a ¼ 0.1853, k ¼ 0.7941 ), collapse load factor l ¼ 49.700.

5.1.2. Load capacity of shallow foundations on Hoek– Brown materials This example deals with the determination of load capacity of foundations in 2D on materials obeying the Hoek–Brown yield criterion (32), extensively used in Rock Mechanics. Fig. 11a shows the geometry of the problem to be solved as well as the finite element discretization with 4-noded isoparametric elements for the virtual velocity field. A constant stress distribution is assumed once again in this case. In the Fig. 11b, the deformed mesh given by the LIMAG code solution is depicted. In the same figure, the

ARTICLE IN PRESS A.F. Durand et al. / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 408–425

419

85 80 75 q u (MPa)

70 B/L =1/1

65 60 B/L = 1/2 B/L = 1/3

55 50

B/L = 1/5

45 40

0.2

0.4

0.6

0.8

1

1.2

B/L

De Beer (1970) LIMAG

Square footing (Das, 1995)

Fig. 9. Rigid rectangular footing (smooth base). Comparison of collapse load factors using different solutions (c ¼ 1.00 MPa, f ¼ 351).

Fig. 11. Strip, flexible shallow footing (smooth base) with inclined loads (b ¼ 201). Solution given by the optimizer MINOS using Hoek–Brown yield criterion (m ¼ 0.621, s ¼ 0.0044 , sc ¼ 20 MPa), collapse load factor lq ¼ 6.81. (a) Initial mesh (400 nodes and 378 element) and (b) vertical stress distribution (MPa).

LIMAG

Serrano & Olalla (1994)

qu (MPa)

20

Fig. 10. Shallow strip footing (smooth base). Deformed mesh showing failure mechanism of case B/L ¼ 1/5 from Fig. 9 and vertical stress distribution (MPa); yield function parameters: c ¼ 1.00 MPa, f ¼ 351; collapse load factor l ¼ 46.679.

vertical stress distribution is also shown. The collapse load factor in the imminence of collapse for this problem is l ¼ 6:81. Serrano and Olalla [41] presented a solution for this problem using the slip lines theory. Their solution for the load factor in this case was l ¼ 5:72 MPa. Fig. 12 shows the load capacity (qu) as a function of the angle of inclination of the applied load. In the same figure, the results obtained by Serrano and Olalla [41] are included. One notices a good match between the two solutions as the angle (b) varies between 01 and 301 although they depart somewhat from one another for b larger than 201. This discrepancy may due to a number of factors, not perfectly identifiable at present. However, the fact that both solutions, the one presented in the present paper and theme of Serrano and Olalla’s were numerically obtained (bringing with it issues such as discretization, solution procedures, etc.) may be an important reason.

15 10 5 0 0

5

10

15

20

25

30

ß (°) Fig. 12. Load factors for a shallow and strip footing (smooth base) using Hoek–Brown yield criterion (m ¼ 0:621, s ¼ 0:0044, sc ¼ 20 MPa). Analytical [41] and limit analysis solutions.

Additionally, the theoretical principles of limit analysis are not coincident with the ones upon which the slip lines theory is based [47]. 5.2. Excavation front of an underground excavation [42] Mu¨lhlhaus [42] presents a simplified analytical solution for the stability analysis of an excavation front in underground structures. Two concentric spheres represent the boundaries of the problem as shown in Fig. 12. This assumption simplifies the solution of the 3D problem as it transformed into one of spherical symmetry. In this case, only one equilibrium equation remains (in polar coordinates). Mu¨hlhaus used the Mohr–Coulomb strength

ARTICLE IN PRESS A.F. Durand et al. / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 408–425

420

criterion, proposing the following system of equations to be solved: qsr s r  sy þ2 ¼ 0, qr r

(61)

sy ¼ Gp sr  sc ,

(62)

Gp ¼

1 þ sinðfÞ , 1  sinðfÞ

(63)

where sy is the circumferential stress, sr the radial stress and sc the uniaxial compressive strength of the material. The solution of the differential equation (61) in connection with Eqs. (62) and (63) yields the radial stress ss in Eq. (64). This stress is the necessary stress to cause yield of the material between spheres S1 and S2 (see Fig. 13). The material parameters used in this example are: cohesion equal to 0.10 MPa and friction angle f equal to 51. Substituting these values in (64) the value 1.258 MPa is obtained for the stress ss . ( )  sc H 0 2ðGp 1Þ ss ¼ 1þ2 0 1 . (64) Gp  1 D For the numerical solution of this problem via NLA, spheres S1 and S2 were substituted by cubes circumscribed to the respective spheres. The Drucker–Prager failure criterion was applied as an approximation of the Mohr– Coulomb criterion. Two different solutions, each one corresponding to a different Drucker–Prager failure surface, were obtained. One surface was generated by circumscribing circles to the hexagon sections to the Mohr–Coulomb surface on the deviatoric planes and the other surface by inscribing circles to the same sections. Once the rules for defining the two Drucker–Prager surfaces are fixed, the corresponding parameters are easily obtained from the Mohr–Coulomb parameters. Two limit analysis problems were then solved, one corresponding to the circumscribed surface and the other to the inscribed surface leading to two different collapse load factors, 1.306 and 1.189, respectively. Fig. 14a and b presents the failure mechanisms and the stress distribution for the problem, calculated by means of

Fig. 14. Muhlhaus’s problem [42] of an excavation front. Solution obtained using optimizer LANCELOT. (a) vertical stress distribution (MPa), circumscribed failure surface, ss  l ¼ 1:306 and (b) vertical stress distribution (MPa), inscribed failure surface, ss  l ¼ 1:189.

the program LIMAG using the two different Drucker– Prager failure surfaces. The failure mechanisms and vertical stresses are in both cases very similar. The magnitudes of the vertical stresses are somewhat higher when the analysis is performed with the circumscribed Drucker–Prager failure surface. The theoretical radial stress obtained with Mu¨hlhaus’s solution (Eq. (64)) lies between the two values obtained numerically ðss ¼ 1:258 MPaÞ.

5.3. Ladanyi and Archambault [43,44] test results of a physical model

Fig. 13. Geometry of the stability problem of the excavation front showing the two spheres (S1 and S2) [42].

In 1972, Ladanyi and Archambault built a physical model for the study of strength characteristics of a fractured rock mass with interlocking. The model was formed by 1800 blocks having square section, each one measuring 1:27 cm
1:27 cm
6:30 cm, sawn from commercial concrete bricks. The blocks were placed forming a persistent (primary) and a non-persistent set of discontinuities (secondary) as illustrated in Fig. 15.

ARTICLE IN PRESS A.F. Durand et al. / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 408–425

421

Table 1 Parameters used for persistent and non-persistent discontinuities Discontinuities

1—Persistent set 2—Non-persistent set

Parameters f (1) [44]

cs (MPa)

cm (MPa m)

39 39

0.0 0.5

10.0 10.0

0 Fig. 15. Configuration for the physical model by Ladanyi and Archambault [44].

1 f1 Bf C @ 2 Ap0;

(67)

f3

f ðr^ c ÞIntact Rock p0:

(68) c

Fig. 16. Collapse loads obtained by limit analysis of the experiments carried out by Ladanyi and Archambault (1972) [43] and approximate analytical solution [44].

The tests carried out consisted in confined compression tests of the assembly of blocks. Ladanyi and Archambault obtained collapse curves for the experiments as function of angle b and, the confining stress [44] as seen in Fig. 16. The observed failure mechanisms in the experiments were sliding along the persistent fractures, and coupled mechanisms of sliding, toppling and failure of intact rock when the non-persistent fractures were involved. Limit analysis of the same geometry was performed by modelling the rock mass as a Cosserat continuum as described in Section 4.3. The mathematical programming problem formulated in order to obtain the collapse load, using the Cosserat representation is stated by Eqs. (65)–(68): Maximize

Z¼l

(65)

subject to Cr^ c ¼ qG jp þ qG ja l,

(66)

where l is the load factor to be maximized, r^ is the vector  of stresses in the domain, C the equilibrium matrix, qG p the permanent surface loads (referred to the confining stress s3 ) or passive loads, qG a the unit active or variable surface forces (referred to stress s1 ) or active loads, f1, f2, f3 are the linear failure criteria of fractures in shear, tension and toppling, respectively, f ðr^ c ÞIntact Rock is the Mohr– Coulomb failure criterion for the intact rock. For the intact material, modelled with the Mohr– Coulomb criterion, the parameters used in the numerical analyses were cohesion 5.5 MPa and friction angle 451. These values were taken from Ladanyi and Archambault [43,44]. Table 1 presents strength parameters for the persistent and non-persistent discontinuities (see Section 4.3 for the definition of these parameters). The values of cm and cs were not explicitly available in [43] and they were obtained through adjustment in order to match the experimental results. The obtained results are shown in Fig. 16 and correspond to one of the tests carried out by Ladanyi and Archambault [43]. In these tests, the confining stress was equal to 0.35 MPa. The figure shows that the collapse loads obtained numerically are close to the experimental results in the range of angles 10ob1o751. In this range of angles, experiments indicated that the failure mechanisms presented the formation of kink bands and shear zones as shown in Fig. 16. In both cases rotation of the blocks present must have had a large influence on the results. The numerical model appeared, however, to be non-conservative in angles b close to 01 and 901. For these angles, the numerical model predicts failure at the intact rock, and the discontinuities have no influence. The figure also shows Hoek’s predictions [44] solely based in sliding modes at the persistent and non-persistent fractures. These predictions overestimate the strength of the block system due to the non consideration of block rotations. The results obtained showed that for a rock mass formed by blocks, the inclusion of rotational modes could be very important in the prediction of collapse. Furthermore, the results show that limit analysis of the Cosserat representation of a fractured rock mass may provide reliable

ARTICLE IN PRESS A.F. Durand et al. / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 408–425

422

estimates for the collapse of the rock structures when rotational modes of failure are relevant. 5.4. Barla et al. [45] test results of a physical model Barla et al. [45], used a physical model to determine the angle of inclination y at collapse of a block configuration shown in Fig. 17a. Displacements were monitored by sensors and a video camera so that the collapse angle could be precisely determined at y ¼ 9:851. An elastoplastic analysis done by the authors provided an approximated value of y ¼ 9:01. The present work presents a limit analysis of the same block configuration used in the physical model (120 square blocks) having a friction angle of 281 and specific weight of 28 kN/m3. These values are given in [45]. Two procedures were adopted in order to perform NLA of the block geometry depicted in Fig. 17a: a discrete block approach (see Section 4.1) [23] and using Cosserat continuum (see Section 4.3) to simulate the blocks as an equivalent continuum. In the discrete approach all blocks composing the physical model were included in the numerical model.

The procedure followed in performing NLA in both cases consisted in maximizing the specific weight of block assembly and varying the inclination angle of the set of blocks. When the load factor reached unity, the corresponding inclination angle was the one being sought. Table 2 shows the obtained relationship between the load factor and the inclination angle. Fig. 17b shows the results obtained with both formulations. Regarding the values obtained with the Cosserat continuum representation, it appears that the load factor obtained numerically depends somewhat on the value of the moment cohesion (cm) defined in Section 4.3. However, for 0:0ocm ðMPa mÞo0:2, the range of inclination angles (in degrees) corresponding to failure was 8:5oyð1Þo11:0, reasonably close to both the experimentally determined value and the numerical solution obtained with elastoplastic analysis. The inclination angle obtained using a discrete approach [23] is also close to the experimental results.

5.5. Problem of symmetric triangular roof prism in an underground excavation In this example, NLA in its discrete approach is applied to the prediction of the stability of rock blocks on the roof of underground excavations. The formulation presented in Section 4.1, expressed by the Eqs. (5)–(7), is used in this case. In Fig. 18, the geometry of the discrete problem to be solved via NLA is shown. The problem consists in the evaluation of the stability of a symmetric and prismatic block (block 3), located on the roof of an underground excavation whose dead weight may lead to the collapse of

θ (a)

Table 2 Load factor as function of the inclination of the slope

4

y1 l

3.5 Friction angle φ = 28° 28

Load Fator λ

3

3.01 3.33

9.251 1.07

9.751 1.02

10.01 0.988

2.5 2 1.5 1 0.5 0 2

4

6

8 Angle θ°

Blocks Cosserat Continuum (Macias. 1997) (cm = 0.0 MPa-m)

10

12

14

Cosserat Continuum (cm = 0.20 MPa-m)

(b) Fig. 17. (a) Slope geometry for Barla et al. [45] physical model. (b) Load factor versus slope inclination y, for discontinuum and Cosserat continuum formulations.

Fig. 18. Symmetric prism on the roof of an underground excavation (a ¼ 301) submitted to a normal confinement force N.

ARTICLE IN PRESS A.F. Durand et al. / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 408–425

the system. The friction forces acting on both faces of the block and the confinement force produced by the rock mass represent the stabilizing forces. Two different forms were adopted for the confinement stress acting on the faces of the prism. As a first situation, the confinement stress is considered acting perpendicular to the lateral faces of the prism (Fig. 18). In sequence, the confinement stress is considered acting horizontally (Fig. 19). The limit equilibrium solutions for both cases are presented by Brady and Brown [46]. Eqs. (69) and (70) show the value of the maximum force due to the dead load of the prism in the imminence of collapse. P‘ ¼ 2N sec f sinðf  aÞ,

(69)

P‘ ¼ 2H tanðf  aÞ,

(70)

where N is the force acting perpendicular to the lateral faces of the prism, H the horizontal force applied in the lateral faces of the prism, f the friction angle of the fracture and 2a the angle of the symmetric prism.

Fig. 19. Symmetric prism on the roof of an underground excavation (a ¼ 301) submitted to a horizontal confinement force H.

423

Considering an angle a ¼ 301 for the prism, a friction angle f ¼ 401, cohesion c ¼ 0:0 MPa, force N ¼ H ¼ 10 kN and substituting these values in Eqs. (69) and (70), load factors equal to 4.534 kN and 3.526 kN are obtained, respectively. For NLA, the forces N and H, are applied at the centroid of the blocks 1 and 2. The results given by the program LIMAG-BLOK [32] are summarized in Tables 3 and 4. It was found that the results concerning load factors depend on the specification of the type of degrees of freedom allowed for blocks 1 and 2. The results correspond to three different specifications for the degrees of freedom in blocks 1 and 2. As an example from Tables 3 and 4 (case II), blocks 1 and 2 are free to displace in the vertical and horizontal directions (displacements u and v are free) but unable to rotate. The collapse load factors obtained with NLA for this case are l ¼ 4:54 and l ¼ 3:52, very close to the values obtained using Eqs. (69) and (70). An explanation is required for the results obtained in cases I and III from Tables 3 and 4. In case I, blocks 1 and 2 are free to rotate and displace. The load factor obtained in both cases show that block 3 is inherently unstable in this situation, any load will displace it from the rock mass. In case III the explanation is somewhat more involved. A basic underlying assumption of limit analysis is the validity of an associate flow rule (see Section 2.1). Regarding the analysed problem, this means that for block 3 to slide out of its initial position, blocks 1 and/or 2 must have degrees of freedom that permit the dilatancy which appears as a consequence of the flow rule. In case III, no movement is allowed for these blocks and, therefore, an infinite load factor is obtained in the analysis. Fig. 20,

Table 3 Load factor for the detachable roof block of Fig. 18 determined by the program LIMAG-BLOK Case

I II III

DOF of blocks 1 and 2 u

v

$

Free Free Fixed

Free Free Fixed

Free Fixed Fixed

Brady and Brown (1999)

l limit analysis

– 4.53 –

1.67 4.54 Infinity

Angle a ¼ 301, friction angle f ¼ 401, cohesion c ¼ 0.0 MPa and normal force N ¼ 10 kN. Displacements u and v correspond to respectively the horizontal and vertical directions. Block rotations are designated by w.

Table 4 Load factor determined by the program LIMAG-BLOK for the detachable roof block shown in Fig. 19 Case

I II III

DOF of blocks 1 and 2 u

v

$

Free Free Fixed

Free Free Fixed

Free Fixed Fixed

Brady and Brown (1999)

l limit analysis

– 3.53 –

Inconsistent 3.52 Infinity

Angle a ¼ 301, friction angle f ¼ 401 and horizontal confinement force H ¼ 10 kN. Displacements u and v correspond to respectively the horizontal and vertical directions. Block rotations are designated by w.

ARTICLE IN PRESS 424

A.F. Durand et al. / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 408–425

Fig. 20. Results of the numerical simulation corresponding to the case shown in Fig. 18. Velocities (normalized) at the imminence of failure according to boundary conditions of case II of Table 3.

obtained from the kinematics of the failure mechanism, shows that indeed for block 3 to displace down (in case II), block 2 had to displace along the vertical direction. 6. Discussion and conclusions The present work presented formulation, implementation and validation of numerical limit analysis (NLA) procedures for the study of stability problems in soil and rock masses. These media were considered as both standard and Cosserat continua and as a discontinuum (discrete blocks). NLA provides a value for the collapse load factor, roughly equivalent to the safety factor obtained in limit equilibrium analyses. Regarding analysis of continuum problems (standard continuum and Cosserat continua), the numerical formulation used proved to be relatively easy to implement, allowing for different interpolation functions for stresses and velocities in the elements. The results obtained showed a clear superiority of the constant stress-bilinear velocity element regarding the consistency of the stress and velocities distribution. The developed programs for this type of analysis, LIMAG and LIMAG/COS determine directly the collapse load factor and the stress field. The solution of the dual problem gives the nodal velocities, which permits the inference of the failure mechanisms involved. LIMAG showed to be very efficient when applied to 2D problems (shallow foundations) considering three different yield/failure criteria (Mohr–Coulomb, Hoek and Brown and Drucker–Prager). Its application to 3D problems was satisfactory only when Drucker–Prager failure criterion was used. For practical purposes, this should not be a problem as Mohr–Coulomb criterion can be approximated by a Drucker–Prager criterion. The application of LIMAG/COS to stability problems of fractured media modelled as a Cosserat continua was particularly successful. As demonstrated by other authors, the representation of fractured media by Cosserat continua is a clear option when modelling such materials. Furthermore the application/validation examples presented in the paper demonstrated that NLA is a feasible, precise and fast

option to determine collapse load factors in such media. NLA analysis of Cosserat continua was indeed able to incorporate and reflect the dominating failure mechanisms of sliding and toppling. The definition of moment cohesion for such media remains a problem. Although only 2D problems were analysed, the extension to 3D should pose no problems. The application of NLA to discrete blocks simulating a discontinuum rock mass is relatively simple provided the geometry of the blocks is known. The developed program to perform NLA of discrete blocks, program LIMAG BLOK deals automatically with the generation of contacts between blocks. The generated problem is a linear problem of easy solution. Two examples (in 2D) were presented to demonstrate the capabilities of NLA when applied to discrete blocks. In the first (see Section 5.4), NLA results using the discrete approach were compared to results obtained with a physical model and also with NLA results using Cosserat continuum. The second example consisted in the determination of the stability of a block at the roof of an underground excavation (see Section 5.5) to which an analytical solution is available. Results obtained with the discrete approach of NLA in both cases were quite satisfactory. Extensions to 3D problems in this case are also relatively straightforward. However, some theoretical issues remain and should be addressed in the future. One such issue that deserves investigation concerns the influence observed in the analyses carried out of the type of element on the quality/coherence of the stress fields obtained. Simpler elements (constant stress) appeared to give better results than more sophisticated ones. Another issue concerns the suitability of the existing optimizers to the solution of large problems especially in 3D situations. As the size of the problem, translated into number of elements, grow, the computational effort for solution grows almost exponentially. It is believed that especially designed optimizers will have to be developed for that purpose. However the overall results reported in the present work permit to state that NLA can be regarded as a very promising tool to use in the study of stability problems in soil and rock masses when represented as continua and discontinua. Acknowledgements Most of the work presented in the present paper was carried out during the development of the PhD dissertation of the first author at the Department of Civil Engineering of Catholic University of Rio de Janeiro, Brazil. During that time he received a grant from CNPq, the Brazilian National Research Council whose support the authors wish to acknowledge. References [1] Chen WF. Limit analysis and soil plasticity. Amsterdam: Elsevier; 1975.

ARTICLE IN PRESS A.F. Durand et al. / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 408–425 [2] Drucker DC, Prager W. Soil mechanics and plastic analysis on limit design. Q Appl Math 1952;10:157–65. [3] Lysmer J. Limit analysis of plane problems in soil mechanics. J Soil Mech Found Div ASCE 1970;96:1311–34. [4] Chen WF, Liu XL. Limit analysis in soil mechanics. Amsterdam: Elsevier; 1990. [5] Finn WDL. Application of limit plasticity in soil mechanics. J Soil Mech Found Div ASCE 1967;101–20. [6] Anderheggen E, Kno¨pfel H. Finite element limit analysis using linear programming. Int J Solids Struct 1972;8:1413–31. [7] Bottero A, Negre R, Pastor J, Turgeman S. Finite element method and limit analysis theory for soil mechanics problems. Comput Methods Appl Mech Eng 1980;22:131–49. [8] Christiansen E. Computation of limit loads. Int J Numer Anal Methods Eng 1981;17:1547–70. [9] Munro J. Plastic analysis in geomechanics by mathematical programming. In: Martins B, editor. Numerical methods in geomechanics. Dordrecht: D. Reidel Publishing Company; 1982. p. 247–72. [10] Casciaro R, Cascini L. A mixed formulation and mixed finite elements for limit analysis. Int J Numer Anal Methods Eng 1982;18:211–43. [11] Sloan SW. Lower bound limit analysis using finite elements and linear programming. Int J Numer Anal Methods Geomech 1988;12: 61–77. [12] Sloan SW. Upper bound limit analysis using finite elements and linear programming. Department of Civil Engineering and Surveying, University of Newcastle, Australia: Report No. 025.09.1987; 1987. [13] Arai K, Jinki R. A lower-bound approach to active and passive earth pressure problems. Soils Found 1990;30(4):25–41. [14] Assadi A, Sloan SW. Undrained stability of shallow square tunnel. J Geotech Eng Div ASCE 1991;117(8). [15] Chuang PH. Stability analysis in geomechanics by linear programming I: formulation. J Geotech Eng Div ASCE 1992;118(11). [16] Singh DN, Basudhar PK. Optimal lower bound bearing capacity of strip footings. Soils Found 1993;33(4). [17] Arau´jo LG. Numerical study of stability problems in rock mass using limit analysis. PhD thesis, Pontifı´ cia Universidade Cato´lica do Rio de Janeiro, Rio de Janeiro, Brazil; 1997 [in Portuguese]. [18] Yu HS, Salgado R, Sloan SW, Kim JM. Limit analysis versus limit equilibrium for slope stability. J Geotech Geoenviron Eng ASCE 1998;January:1–11. [19] Lyamin AV, Sloan SW. Upper bound limit analysis using linear finite elements and non-linear programming. Int J Numer Anal Methods Geomech 2002;26:181–216. [20] Pontes IDS. Non-linear limit analysis in geotechnics problems. PhD thesis, COPPE-UFRJ (Federal University of Rio de Janeiro), Rio de Janeiro, Brazil; 1993. 172pp. [in Portuguese]. [21] Durand AF. Applications of limit analysis as conventional continuum and Cosserat continua to geotechnics problems. PhD thesis, Pontifı´ cia Universidade Cato´lica do Rio de Janeiro, Rio de Janeiro, Brazil; 2000 [in Portuguese]. [22] Livesley RK. Limit analysis of structures formed from rigid blocks. Int J Numer Method Eng 1978;12:1853–71. [23] Macı´ as AJE. Computational implementations for the study of the stability of fractured rock masses. Msc dissertation, Pontifı´ cia Universidade Cato´lica do Rio de Janeiro, Rio de Janeiro, Brazil; 1997 [in Portuguese]. [24] Chou PC, Pagano NJ. Elasticity tensor, dyadic, and engineering approaches. New York: Dover; 1992.

425

[25] Cosserat E, Cosserat F. The´orie des corps de´formables. Paris: Hermann et fils; 1909. [26] Chappell BA. Load distribution and redistribution in discontinua. Int J Rock Mech Min Sci Geomech Abstr 1979;16:391–9. [27] Unterreiner P. Contribution a` l’e´tude et a` la Mode´lisation Numerı´ que deˆs Sols Cloue´s: Application au Calcul en De´formation des Ouvrages de Soute`nement, 2 vols. PhD Thesis, E´cole Nationale des Ponts et Chausse´es, Paris; 1995. 775pp. [28] Eringen AC. Mechanics of micromorphic continua. Mechanics of generalized continua—IUTAM symposium. Freudenstadt & Stuttgart, Kro¨ner: Springer, 1968. p. 18–35. [29] Dawson EM. Micropolar continuum models for jointed rock. PhD thesis, University of Minnesota, Minneapolis, USA; 1995. [30] Besdo D. A plasticity theory of frictionless sliding systems of rocks basing on the Cosserat model. In: Desai, et al., editors. Second conference on constitutive laws for engineering materials: theory and applications; 1987. p. 761–8. [31] Durand AF, Vargas Jr, EA, Vaz LE, Macias JA. Stability of fractured rock masses by limit analysis. In: Ayres da Silva LA, Quadros EF, Gonc- alves HHS, editors. Design and construction in mining, petroleum and civil engineering. Brazil: Sa˜o Paulo; 1998. p. 341–9. [32] Durand AF, Vargas Jr, EA, Vaz LE. Analysis of the stability of discrete blocks through limit analysis (NLA). In: 12th Panamerican conference on soil and rock mechanics, 2003, vol. 1, Cambridge, MA, USA; 2003. p. 1099–1104 [in Spanish]. [33] Maini T, Cundall P, Marti J, Beresford P, Last N, Asgian M. Computer modeling of jointed rock masses. Technical Report N-78-4, US Army Engineer Waterways Experiment Station, Vicksburg, MI; 1978. [34] Schrage L. LINDO—user’s manual. The Scientific Press; 1991. [35] Arora JS. Introduction to optimum design. New York: McGraw-Hill; 1989 624 pp. [36] Murtagh BA, Saunders MA. MINOS 5.1 user’s guide. Technical report SOL 83-20R, Stanford University, California; 1987. [37] Conn AR, Gould NIM, Toint Ph L. LANCELOT—a Fortran package for large-scale nonlinear optimization (Release A). Berlin: Springer; 1992. 330pp. [38] Bowles JE. Foundation analysis and design. 2nd ed. New York: McGraw-Hill; 1977 750pp. [39] Chen WF, Han DJ. Plasticity for structural engineers. New York, Berlin, Heidelberg: Springer; 1988. [40] Das B M. Principles of foundation engineering. 3rd ed. Boston: PWS Publishing Company; 1995. [41] Serrano A, Olalla C. Ultimate bearing capacity of rock. Int J Mech Sci Geomech 1994;31(2):93–106. [42] Mu¨hlhaus HB. Lower bound solutions for circular tunnels in two and tree dimension. Rock Mech Rock Eng 1985;18:37–52. [43] Ladanyi B, Archambault G. Direct and indirect determination of shear strength of rock mass. AIME Annual Meeting, Las Vegas, NV; 1980, Preprint 80–25. [44] Hoek E. Strength of jointed rock masses. 23rd Rankine Lecture. Ge´otechnique 1983;33(3):187–223. [45] Barla G, Bruneto MB, Gerbaudo G, Zaninetti A. Physical and mathematical modeling of a jointed mass for the study of block toppling. In: Myer, Cook, Goodman, Tsang, editors. Conference of fractured and jointed rock masses; 1995. p. 647–53. [46] Brady BHG, Brown ET. Rock mechanics for underground mining. UK: Chapman & Hall; 1999. [47] Mendelson A. Plasticity: theory and application. Florida: Krieger Publishing Company; 1983.