Numerical model for the zonal disintegration of the rock mass around deep underground workings

Theoretical and Applied Fracture Mechanics 67–68 (2013) 65–73

Contents lists available at ScienceDirect

Theoretical and Applied Fracture Mechanics journal homepage: www.elsevier.com/locate/tafmec

Numerical model for the zonal disintegration of the rock mass around deep underground workings S.C. Li, X.D. Feng ⇑, S.C. Li Geotechnical & Structural Engineering Research Center, Shandong University, No.17923 Jingshi Road, 250061 Jinan, PR China

a r t i c l e

i n f o

Article history: Available online 11 December 2013 Keywords: Zonal disintegration Mechanical damage Dynamic process FLAC3D

a b s t r a c t Deep underground excavation can produce zonal disintegration in the surrounding rock under certain conditions. This phenomenon has mainly been studied in a qualitative manner using in situ investigations or laboratory experiments. Researchers have also derived analytical solutions, which resulted in the formation conditions for zonal disintegration. However, the analytic solutions developed to date are not suitable for practical engineering applications because of the stringent model requirements and complicated boundary conditions. To investigate the mechanism of zonal disintegration, a numerical method is proposed to model the phenomenon. The following concepts are incorporated into FLAC3D via its built-in FISH. First, coal mine (or tunnel) excavation is considered as a dynamic process. Second, the element failure criteria are developed based on the maximum tensile stress criterion and strain energy density theory. Third, the mechanical damage is modeled through a decrease in the multi-step elastic modulus, i.e., the nonlinear stress–strain behavior is approximated by the multi-linear elastic softening model. Two practical cases with zonal disintegration are simulated using the proposed method. Both simulations predict the same number of fractured zones at the same location as those obtained through in situ monitoring. Ó 2013 Elsevier Ltd. All rights reserved.

1. Introduction In the field of rock mechanics, investigations on the evolution of the stress–strain state of the rock mass around the underground workings have primarily used elastoplastic mechanics within the framework of continuum mechanics. As a result, the following regions, from the tunnel outward, apply: fissured zone, plastic zone, and undisturbed or elastic zone. However, for deep rock engineering, geo-mechanical phenomena differ greatly from the shallow rock phenomena, especially when the in situ stress is greater than the uniaxial compressive strength of the rock mass. For example, during the excavation of a tunnel in deep rock, alternate regions of fractured and relatively intact rock mass appear around or in front of the working area, which is called zonal disintegration in related research [1,2]. Academics, primarily from Russia [3–9], have validated the existence of zonal disintegration using both in situ monitoring and model tests. Shemyakin et al. [3–5] observed zonal disintegration in the laboratory and qualitatively analyzed its mechanism. The phenomenon of zonal disintegration was also observed in the gold mines of South Africa and in the deep mines of Ukraine. Similarly, the ‘‘onion-skin-shell-structure’’ phenomenon,

⇑ Corresponding author. Tel.: +34 682443093. E-mail address: [email protected] (X.D. Feng). 0167-8442/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.tafmec.2013.11.005

i.e., systems of concentric fractures surrounding an excavation, which was first identified by Müller et al. [10], was observed using in situ ultrasonic wave velocity tests in the Jinchuan nickel mine in China [11]. Researchers, including Zhou et al. [12,13], Wang et al. [14–17], and Qian et al. [18], have investigated the mechanism of zonal disintegration in deep rocks by modeling the excavation as a dynamic process. The prerequisite for zonal disintegration is that there is a massive release of the elastic energy stored in the surrounding rock as the excavation proceeds, which then causes changes in the stress field of the surrounding rock mass; these changes are affected by the unloading time and release of energy. When the stress in the surrounding rock is greater than the ultimate strength of the rock mass, tensile and compressive failure occurs, which may produce zonal disintegration under certain conditions. Zonal disintegration is a unique phenomenon of deep rock engineering, and a change in stress during excavation is a necessary condition for it to occur. At present, the research related to zonal disintegration of deep rocks focuses primarily on qualitative analysis using field observations and model tests. Some scholars have conducted quantitative research, and the formation conditions of zonal disintegration have been obtained analytically [13]. However, the analytical solution cannot be applied directly in engineering practice since it usually has very stringent requirements on the model and boundary

66

S.C. Li et al. / Theoretical and Applied Fracture Mechanics 67–68 (2013) 65–73

conditions. Field tests and model tests are both very expensive and difficult to repeat. Therefore, it is more convenient to simulate zonal disintegration numerically. In recent years, studies have been performed on the numerical simulation of zonal disintegration. For instance, Qian et al. [18] adopted weak elements to model the growth and coalescence of cracks within the rock mass in the Jinping II Hydropower Station and determined the magnitude and distributions of the fractured zones. Tang et al. [19] performed numerical tests on the mechanism and evolution laws of fracture spacing (zonal disintegration) with RFPA. Gao et al. [20] analyzed the generation and evolution of zonal disintegration of surrounding rock mass in roadways using FLAC. However, to date, zonal disintegration has not been simulated numerically from an energy perspective. Based on the above analysis, in this study, a numerical method is proposed to model zonal disintegration. The element failure criterion is established based on the maximum tensile stress criterion and strain energy density (SED) theory. Elastic damage mechanics is applied to simulate the softening behavior of the rock mass. The response of the rock mass to excavation is modeled as a dynamic process, and the patterns of zonal disintegration can be determined by the unloading time. The zonal disintegration code is developed by the FISH in FLAC3D (FISH is the built-in programming language of FLAC3D that enables the user to define new variables and functions. These functions may be used to extend FLAC3D’s usefulness or add user-defined features. In this study we mainly use FISH to implement the element failure criteria (See Section 2.2) and the unloading curve (See Eq. (21))), which is used to model the number and distribution of the fracture zones in the deep rock mass.

most of the dissipated energy is transformed into the kinetic energy of the rock blocks, subsequently leading to rock burst. Based on the above analysis of the stress–strain curve of a rock mass and motivated by the integration principle, a multi-linear stress–strain curve can be used to approximate the true stress– strain curve, as illustrated in Fig. 2. Stress may increase linearly with strain up to point B. The initial elastic modulus is equal to E1. Starting at point B, there is a permanent degradation of the elastic modulus because of the energy dissipation. The initial elastic modulus, E1, decreases to the effective elastic modulus, E2. To obtain more accurate results, the true stress strain curve is approximated using a multiple reduction of the elastic modulus during the yielding stage, shown as segment BU in Fig. 2. Consequently, the nonlinear computational result is obtained by multi-linear elastic model analysis. The macrocrack formation is initiated when the stress increases to the point of ultimate strength, U. Then, the strain may increase as the stress decreases to zero at point F, as the macrocrack expands to become a macroscopic fracture. During the above process, there is both recoverable elastic deformation and irrecoverable damaged (or yielding) deformation. The irrecoverable damage deformation causes the elastic modulus to decrease. With constant energy dissipation, there is also a decreased effective elastic modulus, such as E3, E4, . . . , E⁄. The computation precision depends on the number of iterations of the effective modulus. The approximation will more closely match the true stress–strain curve as the number of iterations increases. 2.2. The strain energy density criterion Strain energy per unit volume is called SED, which can be expressed as follows:

2. Energy dissipation theory for rock failure

ðdW=dVÞ ¼

Z

Stress (MPa)

40

Strain-softening: Initiation of macrocracks

B Elastic plasticity: Initiation of microcracks

30

ð1Þ

where rij and eij are the stress and strain components, respectively. According to SED theory, crack initiation occurs when the SED (dW/dV) reaches a critical value, (dW/dV)c, a quantity defined as the area under the true stress–strain curve. Because the theory can be generally applied to any material, it can also be applied to strain-softening materials [21]. The initial sections, OUC, from Fig. 2, constitute a typical step in the multi-linear computation and can be computed using the bilinear strain-softening constitutive model [21], as shown in Fig. 3.

T

40 35

Stress (MPa)

The complete stress–strain curve, particularly the segment near the peak strength, closely reflects the characteristics of rock deformation and failure. As shown in Fig. 1, curve OA indicates the closing of pre-existing microcracks. Segment AB is a straight line that indicates the linear-elastic characteristics of the rock mass. From points B to U, the rock mass becomes elastoplastic and material damage occurs with energy dissipation. The dissipated energy is used primarily for the generation of new microcracks, and the modulus of the material decreases simultaneously. The rock mass starts to experience strain softening after the peak point, U. The propagation and coalescence of microcracks result in macrocracking. At this stage, energy dissipation is primarily caused by the propagation of the macrocracks and block movements. For brittle materials, if the energy stored in the rock mass is large and little energy is dissipated during the process of macrocrack growth, then

U

rij deij

0

2.1. Rock failure and energy dissipation

50

eij

30

U

25

3

C E

2

E

E

B

20

E

4

E

1

5

15

E

F

20 10 0 0.0

*

5

Linear elasticity

F

A

0 0.2

0.4

0.6

0.8

Strain (%) Fig. 1. Typical stress–strain curve for rock.

1.0

1.2

0

2

4

6

8

10

12

-4

Strain (10 ) UT The stiffness loss induced by damage Fig. 2. Values of effective elastic modulus taking into account energy dissipation.

67

S.C. Li et al. / Theoretical and Applied Fracture Mechanics 67–68 (2013) 65–73

SED theory considers both global and local energy fields instead of focusing on the stress field of the crack tip to avoid the singularity of the stress, which is of vital significance in solving such problems as crack propagation and material failure. In the numerical simulation, it is easy to determine the damaged elements. Motivated by the idea of the integration method, the use of a multilinear calculation to approximate the nonlinear behaviors avoids the need for nonlinear computations.

4 3

Stress (MPa)

30

U C

25

E

20

3

2.3. Elastic damage constitutive equation

E

15

4

It is generally believed that the nonlinearity of the rock stress– strain curve is caused by the initiation and propagation of microcracks instead of plastic deformation of the intact rock material. The brittleness is apparent, particularly in tension. Therefore, the elastic damage mechanics constitutive equation is suitable for describing the mechanical properties of the mesoscopic elements. Starting in the 1980s, Li et al. [22], Carpinteri [21], Zhu et al. [23], Lan et al. [24], and Lai et al. [25] proposed various models, including the elastic damage model and the elastoplastic damage model. According to equivalent strain theory, the constitutive law of the damaged materials can be obtained as follows:

5

B 0

0

2

4

6

8

F 1

12

14

16

18

-4

Strain (10 ) Area OUCB Absorbed SED (dW dV ) Area OCF Decreased Ultimate SED ( dW dV )*u Area OUC Dissipated SED ( dW dV ) d

r ¼ ð1  DÞEe

Area OCB Recoverable SED ( dW dV ) r Area BCF Additional SED (dW dV ) a Fig. 3. Strain softening constitutive model taking energy dissipation into account.

This model accounts for mechanical damage by decreasing the elastic modulus. For an undamaged material element, the initial ultimate value of the SED, (dW/dV)u, is equal to the area DOUF. At point C, damage occurs and the effective elastic modulus decreases from E3 to E4. The dissipated SED, (dW/dV)d, equals the shaded area DOUC, whereas the recoverable SED, (dW/dV)r, is equal to the area DOCB. The area DBCF represents the additional SED,  (dW/dV)a. The decreased ultimate SED, ðdW=dV Þu , can be expressed as follows:

ð7Þ

where E is the elastic modulus of the undamaged material and D is the damage variable. The damage variable can grow from D = 0 (representing the undamaged state) to D = 1 (completely damaged state). At the start, the meso-elements comply with elasticity, which can be expressed by their elastic modulus and Poisson’s ratio. With an increase in stress, the damage threshold is reached and element damage occurs. 2.4. Incremental representation of the elastic damage constitutive equation The strain increment is decomposed into the elastic strain increment and damage strain increment as follows [24]:



ðdW=dVÞu ¼ ðdW=dVÞu  ðdW=dVÞd ¼ ðdW=dVÞr þ ðdW=dVÞa ð2Þ Based on Fig. 3, the relations between the respective energies can be obtained as follows:

Dei ¼ Deei þ Dedi

ð8Þ

The elastic stress increment due to damage is as follows:

Dri ¼ Sdi ðDeei Þ

ð9Þ

ðdW=dVÞ ¼ 1=2ðre þ ru e  reu Þ

ð3Þ

ðdW=dVÞr ¼ 1=2ðreÞ

ð4Þ

where Sdi is the linear function of the elastic strain increment, Deei . The explicit formula of elastic damage can be derived from Eqs. (7) and (9) as follows:

ðdW=dVÞd ¼ ðdW=dVÞ  ðdW=dVÞr ¼ 1=2ðru e  reu Þ

ð5Þ

~ ij ¼ 2GDeij þ a2 Dekk dij Dr

ð10Þ

ðdW=dVÞu ¼ ðdW=dVÞu  ðdW=dVÞd ¼ 1=2ðru ef  ru e þ reu Þ

ð6Þ

a2 ¼ K  2=3G

ð11Þ

where ru is the ultimate strength, eu is the ultimate strain, and ef is the fracture strain. Two threshold criteria for damage are adopted. One is the maxi~ 3 , of the mum tensile stress criterion, i.e., when the effective stress, r ~ 3 P rt , tensile damage element reaches the tensile strength, rt, or r occurs. The other threshold criterion is the SED criterion, i.e., when the absorbed SED is less than the critical energy density, or (dW/ dV) < (dW/dV)c, there is no damage or energy dissipation, so the effective elastic modulus and yielding SED are equal to their initial values. However, if ðdW=dVÞ P ðdW=dVÞc , damage occurs from the coalescence and propagation of cracks, and the material exhibits strain-softening behavior with energy dissipation. The effective elastic modulus and critical SED both decrease simultaneously. Finally, if ðdW=dVÞ P ðdW=dVÞu , the material is completely damaged and the effective elastic modulus returns to zero.

where dij is the Kronecker symbol and G and K are the effective shear modulus and effective bulk modulus caused by the damage, respectively, which are as follows:

G¼

ð1  DÞE ; 2ð1 þ mÞ

K¼

ð1  DÞE 3ð1  2mÞ

ð12Þ

The evolution equation of the damage variable, D, is defined as follows [25]:

  m  1 ~ei D ¼ 1  exp  m ec

ð13Þ

where ~ei is the effective principal strain at time t + Dt, which is determined by the effective principal stress, ec is the ultimate compressive strain, rc is the ultimate compressive stress, and m can be expressed as follows:

68

S.C. Li et al. / Theoretical and Applied Fracture Mechanics 67–68 (2013) 65–73

  1 Eec m ¼ ln

Start

ð14Þ

rc

The effective principal stress is obtained by the effective stress balance as follows:

r~ 3i  ~I1 r~ 2i þ ~I2 r~ i  ~I3 ¼ 0

Input parameters

ð15Þ

Stress state based on the total strain increment

~ i is the effective principal stress (i = 1, 2, 3) and where r eI 1 ; eI 2 ; and eI 3 are the first, second, and third effective stress invariants modeled on the first, second, and third stress invariants I1, I2, and I3. When the effective stress state at a point is determined, the three quantities, eI 1 ; eI 2 ; and eI 3 , become constant and the effective principal stress is obtained. ~ Ii , at Based on Eqs. (10) and (14), the effective principal stress, r time t + Dt is derived as follows:

r~ Ii ¼ r~ oi þ Dr~ i

The principal stresses

Effective principal stresses

ð16Þ

N

Maximum tensile stress criterion

~ Ii is the effective principal stress at time t + Dt and r ~ oi is the where r effective principal stress at time t . Assuming that the new stress vector satisfies the Mohr–Coulomb plastic flow rule, based on the constitutive theory of FLAC3D, the effective stress in the shear failure zone is as follows [24]:

Y

Energy failure criterion

9

r~ N1 ¼ r~ I1  ks ða1  a2 Nw Þ > = r~ N2 ¼ r~ I2  ks a2 ð1  Nw Þ > r~ N3 ¼ r~ I3  ks ða1 Nw þ a2 Þ ;

ð17Þ

N

Y Mark elements with “cfail”

Mark elements with “tfail”

Damage evolution

Stress state renewal

where

ks ¼

pffiffiffiffiffiffiffi

r~ I1  r~ I3 þ 2c Nu ða1  a2 Nw Þ  ða1 Nw þ a2 ÞNu

Nu ¼

1 þ sin u ; 1  sin u

Nw ¼

1 þ sin w 1  sin w

ð18Þ Convergence by FLAC3D

ð19Þ N

a1 ¼ K þ ð4=3ÞG; a2 ¼ K  ð2=3ÞG

ð20Þ

and c is residual cohesion, u is the residual friction angle, and w is the residual dilatancy angle. The new effective stress in the tensile failure zone is as follows:

Y End

Fig. 4. Procedure for elastic-damage constitutive modeling using FLAC3D.

9

r~ N1 ¼ r~ I1  ðr~ I3  rt Þ aa21 > = r~ N2 ¼ r~ I2  ðr~ I3  rt Þ aa21 > ; r~ N3 ¼ rt

ð21Þ

In the FLAC3D calculation, Hooke’s law is used to determine the stress components and total strain increments in the elastic state. If the calculated elastic stress state meets the element failure criterion, the element will be marked using the tags ‘‘tfail’’ or ‘‘cfail’’, which indicate that element damage has occurred. Then, the elastic damage constitutive model is used to calculate the stress of each load step. At the onset of plastic flow, the Mohr–Coulomb criterion is used to calculate the stress. This flow is implemented in FLAC3D as shown in Fig. 4. 3. The numerical formulation of zonal disintegration After excavation, redistribution of the in situ stress conditions within the rock mass surrounding a tunnel or mine roadway will result in a new equilibrium state being established. A complete numerical simulation must model all stages of the adjustments made by the rock mass. The initial state is static and becomes dynamic if the energy released is very large. For example, blasting would usually release more energy than is released in the course of machine excavation. The implementation of the dynamic

calculation proceeds as follows. The reverse load, F, which is perpendicular to the tunnel profile, is applied immediately after excavation. F is specified as follows:

( F¼

pt=TÞ F m  1:0þcosð2 t 6 2T 2

0

t > 2T

ð22Þ

where Fm is the peak load, T is the unloading cycle, and t is the unloading time. Then, the dynamic unloading calculation begins, in which the stress state is judged against the damage criteria. The final equilibrium state is achieved with a transformation from the dynamic state to the static state. The simulation steps are shown in Fig. 5. The above concept is incorporated into FLAC3D using the FISH for predicting zonal disintegration. Numerical results are provided to demonstrate the utility and robustness of the proposed method. As shown in the above analysis, zonal disintegration occurs during the dynamic process. The released energy from the excavation disturbance is sufficient to result in the dynamic behaviors of the surrounding rock mass and is of vital importance to produce zonal disintegration.

S.C. Li et al. / Theoretical and Applied Fracture Mechanics 67–68 (2013) 65–73

Equilibrium of the in situ stress (static process)

The surrounding rock is in the elastic state

The excavation unloading is considered a dynamic process to simulate the true construction process. A reverse load is applied on the contour of the tunnel when the excavation is completed. The energy loss should be included in the process.

Energy dissipation of the surrounding rock Fig. 5. Procedure for numerical modeling progress for determination of zonal disintegration.

4. Numerical simulation of the zonal disintegration of surrounding rock 4.1. Simulation of a Oktyabr’skil mine in situ observation experiment In the in situ observation experiment [3] used for comparison purposes, the numerical model has dimensions of 40  40  1.0 m, with a straight-wall-top-semi-circular-arch tunnel profile with 2 m high walls. The model consists of 12,451 elements and 21,356 nodes. The surrounding rock is sandy mudstone and has the following properties: initial elastic modulus, E = 3 GPa, Poisson’s ratio, l = 0.20, density, q = 2225 kg/m3, buried depth = 957 m, initial hydrostatic in situ stress = 21.3 MPa, excavation unloading time = 5 ms, unloading peak strength = 21.3 MPa, rc = 13.8MPa, ec = 6.0  103, and rt = 2.5 MPa. The unloading curve (Eq. (21)) is shown in Fig. 6 for an unloading time of 10 ms. The surrounding rock has three fractured zones,

69

as shown in Fig. 7, where ‘‘cfail’’ denotes compressive yielding failure and ‘‘tfail’’ denotes tensile failure. The tensile fissured zone near the side wall of the tunnel is 0.8 m wide, and the yielding zone around the arch crown and arch foot spreads to 1.8 m. The two remaining fissured zones are tensile and compressive yielding failure zones with widths of 1.3 m and 0.7 m, respectively. The average spacing distance is 3.8 m between the first and second fissured zones and 1.0 m between the second and third fissured zones. The widths of the fissured zones in the Russian Oktyabr’skil mine in situ observation [3] are respectively 0.9, 1.2, and 1.0 m from the periphery of the tunnel towards infinity, and the spacing distances are 4.0 m and 1.1 m between the first and second fissured zones and between the second and third fissured zones, respectively. Therefore, the computed results agree well with the in situ observations, indicating that the proposed method is valid.

4.2. 3-D simulation of zonal disintegration around a deep roadway of the Dingji mine 4.2.1. Introduction to the in situ monitoring In November 2007, Li et al. [26] used the mine borehole television imager (Type No. KDVJ-400) to observe the internal fissured zones of surrounding rock in a roadway in the 11-2 mining area of the Dingji Mine, located in the Huainan mining area, Anhui Province, China. The buried depth of the observed roadway is approximately 955 m. The mine borehole television imager includes four parts: the underground panoramic camera probe, the system controller, a dedicated TV cable, and a color monitor. The imager (Type No. KDVJ-400) was designed by the China University of Mining and Technology to observe cracks and damage in the rock surrounding boreholes. The inner wall structure of the boreholes can be clearly displayed on the LCD screen. Three 32 mm diameter by 10 m long observation boreholes were drilled to monitor the distribution of fractured zones in a part of the tunnel. Li et al. [26] published the in situ monitoring results for this location, as shown in Fig. 8. The average radii of the four

Fig. 6. Curve for the excavation unload time vs. load at 10 ms.

70

S.C. Li et al. / Theoretical and Applied Fracture Mechanics 67–68 (2013) 65–73

FLAC3D 3.00 Step 7976 Model Perspective 17:22:40 Thu Dec 27 2007 Center: X: 0.000e+000 Y: 0.000e+000 Z: 5.000e-001 Dist: 1.306e+002

Rotation: X: 90.000 Y: 0.000 Z: 0.000 Mag.: 1.25 Ang.: 22.500

Block Group 1 cfail tfail

Itasca Consulting Group, Inc. Minneapolis, MN USA

Fig. 7. Case 1: Fractured zones around the Oktyabr’skil mine tunnel.

Depth (m)

Intact zone Fissured zone

Fig. 8. Case 2:Distribution of the fractured zones in section A of Dingji mine roadway [26].

fissured zones are respectively 0–1.7 m, 2.3–3.4 m, 4.5–5.9 m, and 7.3–7.5 m. 4.2.3. 3-D simulation Rock strata around the roadway in the 11-2 mining district of the Dingji mine are approximately horizontally bedded sandstone. The physical and mechanical parameters of the rock mass are as follows: initial elastic modulus, E = 6.5 GPa, Poisson’s ratio, l = 0.26, density, q = 2,480 kg/m3, initial hydrostatic in situ stress = 23.5 MPa, excavation unloading time = 5 ms, unloading peak strength = 23.5 MPa, rc = 12 MPa, ec = 1.80  103, and rt = 1.5 MPa. Residual cohesion and residual friction angle are as-

sumed to be 0.1 MPa and 30°, respectively. The dimensions of the numerical model are 25  25  15 m, divided into 86,925 elements and 93,696 nodes. The diameter of the upper semi-circle is 5.0 m, and the wall height is 1.38 m. Normal displacement constraints are applied to the left, right, and bottom boundaries, and the upper boundary is a free surface. All of the displacements in the z-direction are fixed. The computed results are as follows. When the excavation unloading time equals 2.5 ms, the first tensile failure zone close to the roadway appears, with an average width of 1.4 m, and the second compressive yielding zone is 3.1 m away from roadway periphery, as shown in Fig. 9. Then, when the unloading time

S.C. Li et al. / Theoretical and Applied Fracture Mechanics 67–68 (2013) 65–73

71

Fig. 9. Case 2: Failure states near Dingji mine roadway for an unloading time of 2.5 ms. (Group ‘‘1’’ denotes the surrounding rock and Group ‘‘2’’ the roadway.)

reaches 5.0 ms, the width of the first fissured zone spreads to approximately 2.0 m, and the third tensile failure zone appears at approximately 4.0 m from the roadway periphery (Fig. 10). When the unloading time reaches 6.5 ms, the fourth tensile failure zone appears at approximately 7.3 m from the roadway periphery (see Fig. 11). As shown in Fig. 11, there are four fissured zones surrounding the roadway when unloading is finished, with the first being the most severely damaged and the fourth being relatively smaller and more dispersed. The first zone ranges from 0.0 m to 2.0 m. The second fissured

zone ranges from 3.2 m to 3.7 m, and small rupture elements are observed in the 2.5–2.7 m area. The third fissured zone ranges from 4.2 m to 5.4 m. The fourth fissured zone ranges from 7.3 m to 7.7 m, with small rupture areas observed from 6.5 m to 6.8 m. The numerical simulation results are compared with the in situ observations [26], as shown in Table 1. The comparisons indicate that the number, width, and position of the fissured zones are consistent with the in situ observations, which indicates the accuracy and feasibility of the proposed analysis.

Fig. 10. Case 2: Failure states near Dingji mine roadway for an unloading time of 5.0 ms.

72

S.C. Li et al. / Theoretical and Applied Fracture Mechanics 67–68 (2013) 65–73

Fig. 11. Case 2: Failure states near Dingji mine roadway for an unloading time of 6.5 ms.

Table 1 Ranges of the disintegration zones (m) determined using the numerical model compared with in situ observations (unit: m). Zone

Position

In situ observation [26]

Observed width

Numerical simulation

Calculated width

I

Inner Outer Inner Outer Inner Outer Inner Outer

0.00 1.73 2.27 3.40 4.50 5.93 7.33 7.53

1.73

0.00 2.04 3.20 3.71 4.25 5.43 7.34 7.67

2.04

II III IV

1.13 1.43 0.20

0.51 1.18 0.33

5. Conclusions (1) Based on the assumptions that the initial in situ stress of deep rock mass is hydrostatic and the surrounding rock mass is isotropic, the numerical results obtained are in good agreement with the in situ monitoring results, indicating that the proposed method could generally simulate the generation and evolution of zonal disintegration within deep isotropic rock mass. (2) Zonal disintegration is determined not only by the uniaxial compressive strength and the in situ stress of the rock mass but also by the excavation unloading time. The excavation unloading time has a significant influence on the morphology of zonal disintegration. (3) The zonal disintegration pattern of the rock mass, including the width and quantity of the fissured zones, can be obtained using the proposed numerical method and failure criteria, which would be very useful in the design of support systems, rock bolting, etc. (4) Because the development of the numerical method is based on the FISH of FLAC3D, the method has wide application prospects, in contrast to analytical methods. It can be used directly in practical engineering problems such as design of support systems in deep mines.

Acknowledgments The authors are grateful for the support provided by the 973 Program (No. 2010CB732002), the General Programs of National Natural Science Foundation of China (No. 51179098, 51134001), and the Ph.D. Programs Foundation of the Ministry of Education of China (No. 20120131110031), which made this research possible. References [1] Q.H. Qian, The characteristic scientific phenomena of engineering response to deep rock mass and the implication of deepness, J. East China Inst. Technol. (Nat. Sci.) 27 (1) (2004) 1–5 (in Chinese). [2] Q.H. Qian, The key problems of deep underground space development, in Proceedings of the Key Technical Problems of Base Research in Deep Underground Space Development-the 230th Xiangshan Science Conference, Beijing, 2004, pp. 1–5 (in Chinese). [3] E.I. Shemyakin, G.L. Fisenko, M.V. Kurlenya, et al., Zonal disintegration of rocks around underground workings, part I: data of in situ observations, J. Mining Sci. 22 (3) (1986) 157–168. [4] E.I. Shemyakin, G.L. Fisenko, M.V. Kurlenya, et al., Zonal disintegration of rocks around underground workings, part II: rock fracture simulated in equivalent materials, J. Mining Sci. 22 (4) (1986) 223–232. [5] E.I. Shemyakin, G.L. Fisenko, M.V. Kurlenya, et al., Zonal disintegration of rocks around underground mines, part III: theoretical concepts, J. Mining Sci. 23 (1) (1987) 1–6.

S.C. Li et al. / Theoretical and Applied Fracture Mechanics 67–68 (2013) 65–73 [6] E.I. Shemyakin, M.V. Kurlenya, V.N. Oparin, et al., Zonal disintegration of rocks around underground workings, part IV: practical applications, J. Mining Sci. 25 (4) (1989) 297–302. [7] M.A. Guzev, A.A. Paroshin, Non–euclidean model of the zonal disintegration of rocks around an underground working, J. Appl. Mech. Tech. Phys. 42 (1) (2001) 131–139. [8] F. Borzykh, Features of the zonal disintegration of roof rocks and a coal seam around mine workings, J. Mining Sci. 26 (5) (1990) 418–427. [9] V.N. Reva, Stability criteria of underground workings under zonal disintegration of rocks, J. Mining Sci. 38 (1) (2002) 31–34. [10] L. Müller, G. Sauer, M. Vardar, Dreidimensionale Spannungsumlagerungen prozesse im Bereich der Ortsbrust, Rock Mech. 7 (1978). [11] M.C. He, H.P. Xie, S.P. Peng, et al., Study on rock mechanics in deep mining engineering, Chin. J. Rock Mech. Eng. 24 (16) (2005) 2803–2813 (in Chinese). [12] H.W. Zhou, H.P. Xie, J.P. Zuo, Developments in researches on mechanical behaviors of rocks under the condition of high ground pressure in the depths, Advan. Mech. 35 (1) (2005) 91–99 (in Chinese). [13] X.P. Zhou, Q.H. Qian, Zonal fracturing mechanism in deep tunnel, Chin. J. Rock Mech. Eng. 26 (5) (2007) 877–885 (in Chinese). [14] M.Y. Wang, Z.P. Zhou, Q.H. Qian, Tectonic deformation and failure problems of deep rock mass, Chin. J. Rock Mech. Eng. 25 (3) (2006) 448–455 (in Chinese). [15] M.Y. Wang, C.Z. Qi, Q.H. Qian, Study on deformation and motion characteristics of blocks in deep rock mass, Chin. J. Rock Mech. Eng. 24 (16) (2005) 2825–2830 (in Chinese). [16] M.Y. Wang, Z.P. Zhou, Q.H. Qian, Tectonic, deformation and failure problems of deep rock mass, Chin. J. Rock Mech. Eng. 25 (3) (2006) 448–455 (in Chinese). [17] M.Y. Wang, H. Song, D.L. Zheng, et al., On mechanism of zonal disintegration within rock mass around deep tunnel and definition of ‘‘deep rock engineering’’, Chin. J. Rock Mech. Eng. 25 (9) (2006) 1771–1776 (in Chinese).

73

[18] Q.H. Qian, X.P. Zhou, H.Q. Yang, et al., Zonal disintegration of surrounding rock mass around the diversion tunnels in Jinping II hydropower station, Southwestern China, Theoret. Appl. Fract. Mech. 51 (2) (2009) 129–138. [19] C.A. Tang, Y.B. Zhang, Discussion on mechanism and evolution laws of fracture spacing in rock mass, Chin. J. Rock Mech. Eng. 27 (7) (2008) 1362–1369 (in Chinese). [20] F.Q. Gao, H.P. Kang, J. Lin, Numerical simulation of zonal distrigation of surrounding rock mass in deep mine roadways, J. China Coal Soc. 35 (1) (2010) 21–25 (in Chinese). [21] A. Carpinteri, Mechanical Damage and Crack Growth in Concrete: Plastic Collapse to Brittle Fracture, Martinus Nijhoff Publishers, Dordrecht, 1986. [22] S.C. Li, W.S. Zhu, Damage constitutive equations for energy dissipation and its applications to stability analysis of surrounding rock mass of caverns, Chin. J. Rock Mech. Eng. 24 (15) (2005) 2646–2653 (in Chinese). [23] W.C. Zhu, C.A. Tang, T.H. Yang, et al., Constitutive relationship of mesoscopic element used in RFPA2D and its validations, Chin. J. Rock Mech. Eng. 22 (1) (2003) 24–29 (in Chinese). [24] H. cLan, J.G. Yao, H.X. Zhang, et al., Development and application of constitutive model of jointed rock mass damage due to mining based on FLAC3D, Chin. J. Rock Mech. Eng. 27 (3) (2008) 572–579. in Chinese. [25] Y. Lai, Y.X. Zhang, Study on macro- and meso-damage composite model of rock and crack propagation rule, Chin. J. Rock Mech. Eng. 27 (3) (2008) 534–542 (in Chinese). [26] S.C. Li, H.P. Wang, Q.H. Qian, et al., In-situ monitoring research on zonal disintegration of surrounding rock mass in deep mine roadways, Chin. J. Rock Mech. Eng. 27 (8) (2008) 1545–1553 (in Chinese).