Simulation of failure around a circular opening in rock

Simulation of failure around a circular opening in rock

International Journal of Rock Mechanics & Mining Sciences 39 (2002) 507–515 Simulation of failure around a circular opening in rock A. Fakhimia,b, F...

851KB Sizes 56 Downloads 532 Views

International Journal of Rock Mechanics & Mining Sciences 39 (2002) 507–515

Simulation of failure around a circular opening in rock A. Fakhimia,b, F. Carvalhoc, T. Ishidad, J.F. Labuze,* a

Department of Mineral Engineering, New Mexico Institute of Mining and Technology, Socorro, USA b Department of Civil Engineering, Tarbiat Modarres University, Iran c Department of Mechanical Engineering, Pontificia Universidade Catolica do Rio de Janeiro, Brazil d Department of Civil Engineering, Yamaguchi University, Ube, Japan e Department of Civil Engineering, University of Minnesota, Minneapolis, MN 55455, USA Accepted 7 April 2002

Abstract A biaxial compression test was performed on a sandstone specimen with a circular opening to simulate a loading-type failure around an underground excavation in brittle rock. The axial force and displacements were monitored throughout the failure process, and microcracking was detected by the acoustic emission technique. To model the observed damage zone around the opening, the distinct element computer program, particle flow code (PFC2D), was used. The numerical model consisted of several circular elements that can interact through contact stiffness, exhibit strength through contact bonds and particle friction, and develop damage through fracture of bonds. For the determination of micro-mechanical parameters needed in the calibration process of the computer program, only the macroscopic parameters of Young’s modulus, Poisson’s ratio and uniaxial compressive strength were used. It is shown that PFC2D was capable of simulating the localization behavior of the rock and the numerical model was able to reproduce the damage zone observed in the laboratory test. r 2002 Elsevier Science Ltd. All rights reserved. Keywords: Biaxial test; Acoustic emission; Underground opening; Distinct element method

1. Introduction There has been a substantial effort by engineers to model localization and softening behavior of rock around underground excavations [1]. The development of damage zones from localization of deformation is a challenging problem because of mesh sensitivity involved in the numerical modeling [2]. To remedy this problem, some researchers have proposed the idea of a non-local continuum [3]. In this approach, stress at a point is a function of a weighted average of strains around that point. However, the characteristic size of the averaging zone for strain is an important issue [4]. Preliminary experimental evidence shows that the characteristic size is related to the grain size and fabric of the rock [5]. Another useful approach in modeling the softening behavior of rock is the distinct element method [6,7]. In this approach, the intact rock is modeled as an assembly *Corresponding author. Tel.: +1-612-625-9060; fax: +1-612-6267750. E-mail address: [email protected] (J.F. Labuz).

of small particles or cylinders, which interact with each other through normal and shear springs that can break. The particles can be bonded to each other to model rock-like materials or be unbonded to represent granular materials. In this way, no mesh discretization is involved and, therefore, no mesh sensitivity is present. This approach can be used to model the static or dynamic, large deformation mechanical behavior of rocks. A micro-mechanical discontinuum program from Itasca Consulting Group (Minneapolis, MN) called the particle flow code (PFC2D) can represent rock by a random assembly of circular particles (grains) bonded together in a dense packing [8]. The resulting distinct model can be regarded as a solid that exhibits deformability through contact stiffness, strength through contact bonds and particle friction, and damage through fracture of bonds. The PFC2D model produces complex macroscopic behavior such as strain softening and dilation by simulating the micro-mechanisms without resorting to complex constitutive laws. Increases in computing power make it feasible to model entire boundary-value problems with distinct particles such

1365-1609/02/$ - see front matter r 2002 Elsevier Science Ltd. All rights reserved. PII: S 1 3 6 5 - 1 6 0 9 ( 0 2 ) 0 0 0 4 1 - 2

508

A. Fakhimi et al. / International Journal of Rock Mechanics & Mining Sciences 39 (2002) 507–515

that the constitutive behavior is built into the model; that is, the constitutive response arises from the interaction of a large number of particles, each of which obeys simple contact laws. This approach holds promise for classes of problems where microcracking and localization must be dealt with to design structures. In this paper, a plane strain (biaxial) test was performed on a rectangular prism of Berea sandstone. The specimen was machined with a 14 mm diameter hole to simulate an underground opening, although failure was achieved through loading rather than excavation of the hole. Acoustic emission (AE) was detected by placing sensors on the faces of the specimen exposed to confining pressure. By recording the microseismic events during application of axial load, the failure response of the structure was monitored. In order to simulate the experiment, PFC2D was used. Upon calibration of the model, a numerical test on a specimen with a hole was conducted. The results with the simulated model showed very good agreement with the damage zone observed in the laboratory.

positioned on each 100 mm  100 mm face, and their average was used as the feedback signal within the closed-loop system. Axial displacements were also monitored. Oil pressure was kept at 7.5 MPa throughout the experiment, and a vertical actuator was moved so as to keep a constant lateral displacement rate of 0.001 mm/s.

2. Experimental model

Mij ¼ Cijkl ckl ¼ Cijkl bk nl DA;

The mechanical properties of the Berea sandstone used were a uniaxial compressive strength of 40 MPa, Young’s modulus of 16 GPa, Poisson’s ratio of 0.28, and density of 2600 kg/m3. A rectangular prism of this rock containing a hole was tested using the University of Minnesota Plane-Strain Apparatus [9]. To model an underground opening, a 14 mm diameter circular hole was cut through the width of the specimen, which measured 100 mm along the plane-strain direction with a cross section of 100 mm  40 mm (Fig. 1). Lateral displacements were measured by two LVDTs, one

2.1. Source characterization AE was used to monitor the failure process during the test, such that both locations and source mechanisms were determined. Damage within a quasi-brittle material such as the Berea sandstone is associated with microcracking. Assuming a point-source model, a microcrack can be described by the displacement discontinuity tensor c: ckl ¼ 12ðbk nl þ bl nk Þ DA;

ð1Þ

where b is the displacement vector, n is the microcrack normal, and DA is the microcrack area [10,11]. Furthermore, a moment tensor M is related to the displacement discontinuity tensor by ð2Þ

where C is the material’s stiffness tensor. Interchanging b and n has no effect on M; consequently, some additional information must be used to distinguish these two vectors. Under its eigenvector coordinate system, the moment tensor for an isotropic material can be written as Mk ¼ lcii þ 2mck ;

ð3Þ

where l and m are the Lame constants, ck and Mk are the eigenvalues of the tensors cij and Mij ; and cii ¼ c1 þ c2 þ c3 : For isotropic solids, the eigenvector coordinate systems of M and c coincide. The solution

Fig. 1. Test setup: (a) specimen geometry and AE monitoring; (b) loading configuration.

A. Fakhimi et al. / International Journal of Rock Mechanics & Mining Sciences 39 (2002) 507–515

of the eigenvalue problem for cij yields 2

c1 ¼ jbjDA cos a; c2 ¼ 0; 2

c3 ¼ jbjDA sin a;

jbjDA ¼ c1  c3 ; pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a ¼ 7tan1 ð c3 =c1 Þ;

ð4Þ

where the crack mode angle a is the angle between the normal to the displacement discontinuity n and the first eigenvector of M: The angle 2a is measured from the normal n and the displacement vector b: Obviously, if a ¼ 01; then it is a pure tensile mechanism; if a ¼ 451; it is pure shear. Note that the intermediate value of c is zero, so the moment tensor must have the form l M2 ¼ ðM1 þ M3 Þ: ð5Þ 2ðl þ mÞ The displacement vector u at any point X and time t due to a suddenly produced displacement discontinuity in an infinite domain is conveniently given in terms of the moment tensor [12]: ui ðX ; tÞ ¼ Mjk * Gij;k ;

ð6Þ

where Gij;k is the derivative of a Green’s function that describes the elastodynamic problem, which for AE corresponds to a point-step force at point X0 and time t0 in an infinite domain. Since only the first peak of the displacement is considered in this analysis, the analytic solution to u for a radial displacement at a point due to the P-wave is used [12]. Knowing the displacement at each transducer, the moment tensor of the source, which is assumed not to depend on time, can be obtained by minimizing the sum of the error E in the displacement: E¼

N  X 2 uk  un ðMpq ; X Þ k¼1

  k1 k 3 III þ þ g I  II þ I ; k ðk þ 1Þ2

ð7Þ

where uk is the normal displacement measured at the kth transducer, un ¼ ui ni is the displacement caused by the moment tensor of the source, n is the normal to the plane of the transducer, I–III are the invariants of the moment tensor, and g is the Lagrange multiplier. Note that the moment tensor components are calculated subjected to the restriction that the second eigenvalue of c must be zero to be compatible with the model of displacement discontinuities as the source mechanism. The components of the moment tensor are calculated so as to minimize the error between the measured and calculated displacements, subjected to the imposed restriction. 2.2. Sensor calibration If AE source characterization is to be performed, it is required to relate the transducer output to displace-

509

ments at the surface. For each event recorded during a test, a deconvolution process is necessary to estimate displacements at each transducer location [13], which can become time consuming. Some authors have avoided this process by directly relating transducer output to vertical displacements at the surface through a scalar sensitivity parameter [14,15]. This is usually done by assuming that the first peak of the displacement calculated from an analytic solution, which is due only to the P-wave, is directly proportional to the first peak of the recorded signal. This relation can be written SA ¼ up =Vp ;

ð8Þ

where up is the amplitude of the first displacement peak, Vp is the amplitude of the first voltage peak, and SA is the amplitude sensitivity parameter given in meters/ mvolt. Before the experiments, the piezoelectric (Physical Acoustics model S9225) transducers and preamplifiers were calibrated by breaking a pencil lead (0.5 mm diameter) at known positions on the specimen to produce a point-step force. The load for breaking the pencil lead was measured with a load cell (sensitivity of 2.2 N/mV), and the normal displacement was calculated from the analytic solution for a point-step force applied to a half space [16]. Only the P-wave component was considered, and it was ensured that reflected waves were not distorting the signal. The objective of calibration is to relate the transducer response to displacement taking place at their location. The true response of a sensor, however, is a function of the disturbance generated by the source event modified by the presence of the transducer. Since the undisturbed field is of interest, the loading effect becomes part of the representation. Several simplifying assumptions must be made in order to obtain a practical solution. In general, it is assumed that the wave modes are uncoupled and can be analyzed separately. Transducers are assumed to be vertically sensitive [17]. Also, it is desired to have transducers that are very small compared to the characteristic dimensions of the tests such that their size can be disregarded in the analysis. 2.3. Results A total of eight sensors, four on each lateral face of the top half of the specimen, were used to monitor the AE. During the test, 3150 events were recorded at a maximum rate of 12 events/s. The normalized load–displacement behavior with AE is shown in Fig. 2. Only the events that could be located with error o1 mm were selected for source characterization. The location of 497 events is shown in Fig. 3a. As expected, these locations are biased towards the coverage area of

A. Fakhimi et al. / International Journal of Rock Mechanics & Mining Sciences 39 (2002) 507–515

510

the sensors since the location errors usually increase with the distance from the sensors. From Fig. 3b, it is clear that many of the events were concentrated around two potentially competing damage zones. Fig. 4, which 3500

AE events

Normalized Load (MPa)

40

2800

30

2100

20

1400

10

700

0

Cummulative number of events

50

is a photograph of the specimen at the end of the test, shows the resulting failure mechanism. It can be observed that in addition to the visible cracks, two notches were also developed at the lateral boundaries of the hole due to the stress concentration. The AE events due to these notches were not detected, possibly due to the biased measurements noted earlier. It is believed that the damage in these notches and, therefore, the compliant behavior of the rock in this region was responsible for the initiation of the cracks that developed. These rupture zones, which have kinks, eventually reached to the vertical boundaries. Fig. 5 shows how the orientations of both the displacement vectors b and microcrack planes n vary with respect to the orientation of the failure plane. Vectors b and n were distinguished by examining the global displacements and rejecting the orientation of b that was incompatible based on the kinematics. It was found that most of the displacement vectors were aligned with the failure plane. The microcrack planes, however, were inclined at angles mainly between 501 and 701.

0 0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

Normalized Displacement (%)

100

100

90

90

80

80

70

70

60

60 y (mm)

y (mm)

Fig. 2. Normalized load (force/area) versus normalized displacement (displacement/height) behavior of the specimen together with the AE events.

50 40

40

30

30

20

20

10

10 0

0 0

(a)

50

10

20

30

40

50 x (mm)

60

70

80

90

100

0

10

20

30

40

z (mm)

(b)

Fig. 3. AE locations with o1 mm error. (a) All 497 AE events; the eight AE sensors, positioned in the upper half of the specimen, are indicated by the filled circles. (b) AE locations at different stages of loading.

A. Fakhimi et al. / International Journal of Rock Mechanics & Mining Sciences 39 (2002) 507–515

511

Fig. 4. Photograph of the specimen at the end of the test with the resulting failure plane.

100 Displacement vector Number of events

80

Crack orientation

60

40

20

0 10

20

30

40

50

60

70

80

90

Angle from failure plane Fig. 5. Statistics of the angles between the crack planes and displacement vectors with respect to the macroscopic failure plane.

3. Numerical model PFC2D is a distinct element computer program designed to simulate the mechanical behavior of bonded or unbonded granular materials. The soil or rock material is modeled as a collection of circular particles that can interact through normal and shear springs. Under the applied load, the bonds between particles can break and a small ‘‘crack’’ can form. By further generation of these microcracks,

a fracture can develop from the linking of individual microcracks. PFC2D is a good numerical tool for modeling unlimited dynamic or static movement of rocks. The crack pattern is developed automatically with no need for remeshing. Heterogeneity and planes of weakness, which are observed in rock, can be easily modeled in this program. Although PFC2D can simulate a particulate media, any circular element in this program does not necessarily model a particle in the real material as the two-dimensional nature of the program limits particles to disks or cylinders. So, PFC2D can be considered as an attempt to mimic the basic mechanical features of the actual material. There have been some successful attempts to model the damage zones around underground opening using PFC2D. Potyondy and Autio [18] and Cundall [7] used this program to predict the fracture zones around the Underground Research Laboratory (URL) in Manitoba, Canada. Recently, Potyondy [18] modeled the damage zones around a circular test hole in gnessic tonalite subjected to compressive loading. 3.1. Calibration tests To simulate the mechanical behavior of a particular rock, Berea sandstone, some numerical calibration tests were necessary. Calibration involves finding the micromechanical parameters of the assembly from the macroscopic properties such as Young’s modulus and

512

A. Fakhimi et al. / International Journal of Rock Mechanics & Mining Sciences 39 (2002) 507–515

0.8 0. 6

0. 4 ν'

E'/K n

0. 6

0. 2

0.4 0.2

0. 0

0

0 (a)

0.5

1

1.5

2

0 (b)

K s /K n

0.5

1

1.5

2

K s /K n

Fig. 6. Variation of: (a) normalized modulus ðE 0 =Kn Þ and; (b) Poisson’s ratio with Ks =Kn :

the uniaxial compressive strength. To perform this calibration, the micro mechanical parameters were estimated using the work of Huang [19]. The parameters in PFC2D are: normal and shear stiffnesses of the contact Kn and Ks ; the radii of the particles, the normal and shear strength of the contact bonds Tn and Ts and the coefficient of friction between particles m: By using dimensional analysis and performing unconfined biaxial tests numerically, the macroscopic properties of the sample can be calculated and compared to those of Berea sandstone. The effort was concentrated on calibrating the PFC2D material only for the unconfined compressive strength, Young’s modulus and Poisson’s ratio. As an example, the calibration procedure for the Young’s modulus and Poisson’s ratio is explained. Figs. 6a and b show the variation of normalized modulus ðE 0 =Kn ) and Poisson’s ratio with Ks =Kn ; these curves were obtained from dimensional analysis and several numerical tests [19]. E 0 and n0 are related to the real Young’s modulus E and Poisson’s ratio n of the material by the following relationships: E 0 ¼ E=ð1  v2 Þ;

v0 ¼ v=ð1  vÞ:

ð9Þ

Having E and n for the sandstone, E 0 and n0 can be calculated; n0 can be used in Fig. 6b to find Ks =Kn : For n ¼ 0:3 and therefore n0 ¼ 0:43; the appropriate value of Ks =Kn is 0.15. Using this value, E 0 =Kn is obtained from Fig. 6a. Finally, the value of Kn and Ks can be calculated, knowing the value of E 0 : A similar procedure for calibration of uniaxial compressive strength can be performed [19]. After several simulations, the following micromechanical parameters were obtained for the sandstone: Kn ¼ 64 GPa, Ks ¼ 9:6 GPa, Tn ¼ 30 MPa, Ts ¼ 150 MPa, m ¼ 0:5: A uniform random distribution of Tn and Ts were assumed for the model, with standard deviations of 3 and 15 MPa, respectively. The minimum radius of the particle assembly was Rmin ¼ 0:1 mm and Rmax =Rmin ¼ 3: The particle radii were chosen to have a uniform distribution in the sample. Although a particle in PFC2D does not need to correspond to a real grain in the actual

material, it was decided to have an average particle radius of 0.2 mm, which was equal to the average grain size in the Berea sandstone.

4. Results and discussion In Fig. 7a, the synthetic material with the hole is presented. After creating the sample, the hole was excavated and the lateral walls were adjusted to create a hydrostatic pressure equal to 7.5 MPa. Then, the upper and lower walls were moved slowly to simulate the biaxial test, while the movement of lateral boundaries were controlled to keep the lateral pressure constant and equal to 7.5 MPa .The normalized load–displacement curve for this specimen is shown in Fig. 7b, which can be compared with the results of actual lab test. Although the PFC2D specimen showed a slightly stiffer response, the peak loads matched very well. The damage zone at the end of the numerical test can be recognized by plotting the cracks in the PFC2D specimen (Fig. 8a). A microcrack in PFC2D is a line perpendicular to the centerline of the two particles with the broken bond that has a length equal to the average radii of two particles involved. The damage zone and microcrack patterns were similar to those observed from the laboratory test. The failure plane in both specimens exhibited a kink; the failure plane extended from the notches on the hole periphery in the vertical direction with a gradual incline toward the vertical boundaries. It is interesting to note that in the PFC2D material, there was some competition between the damage zones on the different sides of the specimen, but finally one rupture zone developed by the coalescence of the microcracks. This observation is consistent with the AE events observed (Fig. 3). AE were located on both sides of the hole, although the eventual failure plane was visible on one side only. The final rupture zone can be recognized in Fig. 8b, which shows the displacement field in the specimen. The notches and the displacement discontinuities in the domain are clearly visible. By following the crack

A. Fakhimi et al. / International Journal of Rock Mechanics & Mining Sciences 39 (2002) 507–515

513

50 Experiment

Normalized Load (MPa)

40

Numerical

30

20

10

0 0.0 (a)

(b)

0.1

0.2

0.3

0.4

0.5

Normalized Displacement (%)

Fig. 7. (a) Synthetic specimen with the hole. (b) Normalized load–displacement curve.

Fig. 8. Synthetic specimen with: (a) the microcracks; (b) the displacement field at the end of the numerical test.

generation in the specimen during the failure process, it was observed that the cracks started to develop from the periphery of the tunnel, where the compressive stresses were large. But even during the evolution of these notches, some cracks were formed near the eventual rupture zone. By continuation of loading, more damage was developed within the notches, and the material in these regions started to spall (Fig. 9). The more compliant behavior of the notches stimulated further crack generation close to the rupture zone in both sides of the specimen. Between the two possible failure planes, eventually one of them prevailed and continued to

Fig. 9. Damaged notches near the hole.

propagate. Finally, the specimen showed a macroscopic discontinuity along the rupture zone. Fig. 10 illustrates a close-up view of a portion of the failure plane; the microcracks (solid lines) and the displacement field (arrows) are shown. The figure was divided into three separate regions: A, B and C. In A and C, the displacement fields were almost like rigid body motions, separated by a transition zone. Region B, which included the majority of the

514

A. Fakhimi et al. / International Journal of Rock Mechanics & Mining Sciences 39 (2002) 507–515

crack extension without the need for mesh generation or remeshing. The results of a numerical simulation of a model tunnel were in very good agreement with the experimental test. It was shown that the failure load, the crack pattern, and the spalling of the opening matched very well with the laboratory results. The microcrack pattern in the numerical model clearly showed the competition of two sets of fractures, which was consistent with the locations of AE determined in the laboratory test. The orientation of the microcracks from the numerical analysis showed that they were not completely in the direction of the rupture zone. This was also observed experimentally through source characterization of AE. In addition, assuming that the displacement transition zone was a type of shear band, it was realized that the thickness of this band was in close agreement to that observed in the experiment.

Acknowledgements

Fig. 10. Close-up view of a portion of the failure plane together with the microcracks and displacement field.

microcracks, was the failure plane. From this figure, it is clear that the microcracks were not necessarily in the direction of the macroscopic failure plane. This is consistent with the AE observations (Fig. 5). In addition, the width of region B can be interpreted as the thickness of the rupture zone. This thickness was approximately equal to six times the grain size of the PFC2D material. Since the average particle radius was 0.2 mm, the thickness of the rupture zone was about 2.5 mm. Based on visual observations, the thickness of the rupture zone in the Berea sandstone specimen was about 3 mm.

5. Summary and conclusions Damage and fracture propagation around an underground excavation are important issues in rock engineering. Numerical modeling of these softening phenomena by finite element method needs special considerations for mesh sensitivity. It has been shown that by simulating the rock as an assemblage of cylindrical particles, surface spalling and crack propagation can be easily modeled. An important advantage of the distinct element method in simulating rock behavior is its ability to model finite displacement of particles and

Partial support was provided by the National Science Foundation (CMS-0070062) and Minnesota Supercomputing Institute. Dr. David Potyondy of Itasca Consulting Group helped with the implementation of PFC2D and Dr. Haiying Huang provided her Ph.D. dissertation containing the important calibration results.

References [1] Fakhimi A, Fairhurst C. A model for the time-dependent behavior of rock. Int J Rock Mech Min Sci & Geomech Abstr 1994;3(2):117–26. [2] Rice JR. The localization of plastic deformation. in: Koiter WT, editor. Theoretical and applied mechanics. Amsterdam: NorthHolland, 1976, pp. 202–20. [3] Bazant ZP, Pijaudier-Cabot G. Non-local continuum damage, localization instability and convergence. J Appl Mech 1988;55:287–93. [4] Bazant ZP, Pijaudier-Cabot G. Measurement of characteristic length of non-local continuum. Report 87–12/498 m, Center for Concrete and Geomaterials, Northwestern University, 1987. [5] Zietlow WK, Labuz JF. Measurement of the intrinsic process zone in rock using acoustic emission. Int J Rock Mech Min Sci & Geomech Abstr 1998;35(3):291–9. [6] Cundall PA. A computer model for simulating progressive, large scale movements in blocky rock systems. Symp Int Soc Rock Mech Nancy, France, 1971, pp. 11–8. [7] Potyondy D, Cundall P. The PFC model for rock: predicting rock-mass damage at the underground research laboratory. Itasca Consulting Group, Inc. Report no. 06819-REP-01200-10061-R00, 2001. [8] Itasca Consulting Group, 2001. Particle Flow Code in 2 Dimensions, version 2.0. Minneapolis, Minnesota. [9] Labuz JF, Dai ST, Papamichos E. Plane-strain compression of rock-like materials. Int J Rock Mech Min Sci & Geomech Abstr 1996;33(6):573–84.

A. Fakhimi et al. / International Journal of Rock Mechanics & Mining Sciences 39 (2002) 507–515 [10] Burridge R, Knopoff L. Body force equivalents for seismic dislocations. Bull Seismol Soc Am 1964;54:1875–88. [11] Rice JR. Elastic wave emission from damage processes. J Nondestr Eval 1980;1:215–24. [12] Aki K, Richards PG. Quantitative seismology: theory and methods. San Francisco: W.H. Freeman and Company, 1980. [13] Michaels JE, Michaels TE, Sachse W. Application of deconvolution to acoustic emission signal analysis. Mater Eval 1981;3:1032–6. [14] Scruby CB, Baldwin GR, Stacey KA. Characterization of fatigue crack extension by quantitative acoustic emission. Int J Frac 1985;28:201–22. [15] Ohtsu M. Source kinematics of acoustic emission based on a moment tensor. NDT Int 1989;22(1):14–20.

515

[16] Shah KR, Labuz JF. Damage mechanisms in stressed rock from acoustic emission. J Geophys Res 1995;100(B8):15527–39. [17] Simmons JA, Turner CD, Wadley HNG. Vector calibration of ultrasonic and acoustic emission transducers. J Acoust Soc Am 1987;82:1122–30. [18] Potyondy D, Autio J. Bonded-particle simulations of the in situ failure test at Olkiluoto. in: Elsworth D, Tinucci J, Heasley K, editors. Rock mechanics in the national interest, Proceedings of the 38th US Rock Mechanics Symposium, 2001, pp. 1553–60. [19] Huang H. Discrete element modeling of tool–rock interaction. Ph.D. dissertation, Department of Civil Engineering, University of Minnesota, USA, 1999.