UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

An analytical geomechanical upscaling approach for modeling jointed rock mass behaviour using ubiquitous.. Lavoie, Thierry 2010-01-04

You don't seem to have a PDF reader installed, try download the pdf

Item Metadata

Download

Media
24-ubc_2011_spring_lavoie_thierry.pdf [ 11.47MB ]
Metadata
JSON: 24-1.0052524.json
JSON-LD: 24-1.0052524-ld.json
RDF/XML (Pretty): 24-1.0052524-rdf.xml
RDF/JSON: 24-1.0052524-rdf.json
Turtle: 24-1.0052524-turtle.txt
N-Triples: 24-1.0052524-rdf-ntriples.txt
Original Record: 24-1.0052524-source.json
Full Text
24-1.0052524-fulltext.txt
Citation
24-1.0052524.ris

Full Text

  An Analytical Geomechanical Upscaling Approach for Modeling Jointed Rock Mass Behaviour Using Ubiquitous Joints     by   Thierry Lavoie       A THESIS SUBMITED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF  MASTER OF APPLIED SCIENCE   in   THE FACULTY OF GRADUATE STUDIES  (Geological Engineering)  THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver)         January 2011  © Thierry Lavoie, 2011   ii  Abstract Over the past few years, the increase in scale of open pit mines and a need to more accurately predict the subsidence induced by block cave mining have highlighted the need to develop new analytical techniques to replace empirical rock mass rating systems in order to better evaluate fractured rock mass properties and simulate its behaviour.  This thesis focuses on the development of an analytical geomechanical upscaling approach for modeling jointed rock mass behaviour in continuum simulations based on the information that can be derived from Discrete Fracture Network (DFN) modelling and laboratory test results.  For this research, many approaches have been evaluated using different constitutive models and techniques to derive the rock mass properties.  The Ubiquitous Joint Rock Mass (UJRM) constitutive model has proven to be an ideal tool to capture both the softening effect and directionality imposed by the discontinuities on the rock mass.  A good agreement has been observed between the outcomes from simulations using the upscaling approach in FLAC and similar 2D models run in ELFEN.   The potential for the upscaling approach to accurately reproduce fractured rock mass behaviour was further confirmed by testing its ability to reproduce the scale effect and applying it to four different slope models.  This research indicates that the developed geomechanical approach developed can reproduce the behaviour of fractured rock masses in continuum simulations while necessitating minimum preparation time, being less computationally intensive than its discontinuum counterparts and staying as close as possible to the data acquired in the field and from laboratory testing.   iii  Table of Contents  Abstract ............................................................................ ii Table of Contents ............................................................ iii List of Tables ................................................................... v List of Figures ................................................................. vi 1 Introduction ................................................................. 1 1.1 Problem Statement ........................................................................................................... 1 1.2 Research Objectives ......................................................................................................... 3 1.3 Thesis Organization ......................................................................................................... 4 2 Literature Review ........................................................ 5 2.1 Discrete Fracture Networks (DFN) .................................................................................. 5 2.1.1 Description ................................................................................................................ 5 2.1.2 Brief History ............................................................................................................. 5 2.1.3 Spatial Models .......................................................................................................... 8 2.1.4 Use of DFNs in this Thesis ....................................................................................... 9 2.2 Rock Mass Classification and Characterization............................................................. 10 2.2.1 Rock Mass Classification Systems ......................................................................... 10 2.2.2 Rock Mass Characterization and its Use in this Thesis .......................................... 13 2.3 Brittle Rock Failure Simulations ................................................................................... 15 2.3.1 Implicit Simulations................................................................................................ 15 2.3.2 Explicit Simulations................................................................................................ 17 2.3.3 Synthetic Rock Mass .............................................................................................. 20 2.3.4 Ubiquitous Joint Rock Mass (UJRM) .................................................................... 26 2.4 Summary ........................................................................................................................ 28 3 Methodology ............................................................. 29 3.1 FLAC ............................................................................................................................. 29 3.1.1 Finite Difference Grid............................................................................................. 29 3.1.2 Constitutive Models ................................................................................................ 31 3.1.3 Boundary Conditions .............................................................................................. 32 3.2 ELFEN ........................................................................................................................... 33 3.3 Phase² ............................................................................................................................. 33 iv  4 Upscaling Approach Development ............................ 36 4.1 Methodology .................................................................................................................. 36 4.1.1 Programming .......................................................................................................... 36 4.1.2 4x8 Models (Initial Testing) ................................................................................... 37 4.2 Representing the Fractured Rock Mass with Different Constitutive Models ................ 38 4.2.1 Mohr-Coulomb ....................................................................................................... 38 4.2.2 Ubiquitous Joint Rock Mass (UJRM) .................................................................... 43 5 Validation Testing of the Upscaling Procedure ......... 63 5.1 8x4m Pillars ................................................................................................................... 64 5.2 20x20m Blocks .............................................................................................................. 68 5.3 Reproduction of the Scale Effect and Failure Path ........................................................ 80 6 Application of Upscaling Approach to Pit Slope Models ............................................................................ 90 6.1 DFN Models ................................................................................................................... 90 6.2 FLAC Ubiquitous Joint Rock Mass Models .................................................................. 93 6.2.1 Geometry and Properties ........................................................................................ 93 6.2.2 Most Critically Oriented Joint with Length Component Upscaling Approach ...... 95 6.3 Comparison of FLAC model with ELFEN and PHASE² ............................................ 106 6.3.1 Horizontal Displacement ...................................................................................... 106 6.3.2 Horizontal Stresses and Tension ........................................................................... 111 7 Conclusion and Recommendations ......................... 116 7.1 Conclusion ................................................................................................................... 116 7.2 Recommendations for Further Work ........................................................................... 118 References .................................................................... 119 Appendix A: DFN Models Review .............................. 128 Appendix B: Rock Mass Rating Systems Review ....... 133 Appendix C: C++ Programming .................................. 145 Appendix D: Slope Models .......................................... 174    v  List of Tables Table 2-1 Fracture intensity, density and porosity measurements (modified after Golder Associates, 2008) ............................................................................................................................ 9 Table 2-2   History of rock mass rating systems (modified after Edelbro, 2003) (continued) ..... 12 Table 2-3 Parameters included in different system (modified after Edelbro, 2003) .................... 12 Table 4-1 Intact rock and fracture properties for the two 4x8m models ...................................... 38 Table 4-2 Intact rock and fracture properties for the Raster Models ............................................ 39 Table 4-3 Intact rock and fracture properties for the 4x8m UJRM models.................................. 44 Table 5-1 Intact rock and fracture properties for the eight 20x20m models ................................ 64 Table 5-2 Intact rock and fracture properties for the 60x60m model ........................................... 81 Table 6-1 Intact rock and fracture properties for the slope models .............................................. 94    vi  List of Figures Figure 2-1Quantification of GSI chart and block delimited by three joint sets (Cai et al., 2004) 14 Figure 2-2 Simulated fracture process of a rock specimen under uniaxial compression ............. 16 Figure 2-3 Sensitivity of failure modes to local variation for five specimens with same parameters of mechanical properties final failure modes (simulated with RFPA2D) (Tang, 2000a) ...................................................................................................................................................... 17 Figure 2-4 Confined compression test on a bonded disc sample. Cracks show as segments normal to ruptured bond. Right-hand image is a manually generated schematic of major rupture coalescence based on individual cracks in middle figure (Diederichs, 2003) .............................. 18 Figure 2-5 Simulated anisotropic effects induced in a fractured rock mass by varying the angle between the applied principal stress direction and the inclination of the predefined fractures (Elmo and Stead, 2010) ................................................................................................................ 19 Figure 2-6 Modes of failure for selected simulated pillars with 0.43m mean fracture length (Elmo and Stead, 2010) ................................................................................................................ 20 Figure 2-7 Scale effect on rock-mass compressive strength (modified after Sjöberg, 1999)....... 21 Figure 2-8 Illustration of how the same rock mass could appear blocky for an underground excavation, but disintegrated for a slope (Lorig, 2007) ................................................................ 21 Figure 2-9 Rock mass sample: before (left) and after (right) numerical triaxial test (circles indicate tensile failure locations) (Lorig, 2007) ........................................................................... 22 Figure 2-10 Stress damage in failed rock mass sample (modified after Clark, 2006).................. 23 Figure 2-11 Comparison of ubiquitous joint and Mohr-Coulomb rock mass behavior (modified after Clark, 2006) .......................................................................................................................... 23 Figure 2-12 Design chart to estimate: a) the rock mass strength and b) the secant modulus based on GSI rock mass classification system (modified after Clark, 2006) ......................................... 24 Figure 2-13 Two-dimensional SRM model with PFC2D with bumpy joints (modified after Park et al., 2004) ................................................................................................................................... 24 Figure 2-14 Use of three-dimensional SRM in evaluating the rock response around a large block cave (Lorig, 2007) ........................................................................................................................ 25 Figure 2-15 Use of three-dimensional SRM in evaluating the rock response around a large block cave (after Sainsbury et al., 2008) ................................................................................................ 26 Figure 2-16 The Palabora pit slope failure mechanism reproduced by the SRM–UJRM approach for representing jointed rock masses (after Sainsbury et al., 2008).............................................. 27 Figure 3-1 FLAC® General solution procedure (after Itasca Consulting Group, 2008a) ............ 30 Figure 3-2  Close-up of Phase2 model representation of discrete fracture network ..................... 35 Figure 4-1 Flow chart for the C++ code ....................................................................................... 36 Figure 4-2 2D 4x8m DFN models a) 10-80 fracture pattern b) 30-60 fracture pattern ................ 37 Figure 4-3 Raster Models (7200 cells) and DFN they are based on: a) 10-80 b) 30-60 .............. 40 Figure 4-4 Raster Models (7200 cells) compared to ELFEN models, where Sigma1 and Sigma3 refer to the major and minor principal stresses, respectively ....................................................... 41 Figure 4-5 Raster Models’ deformation after peak strength a) 10-80 b) 30-60............................ 42 Figure 4-6 Schematic representation of the Reduced Strength “GSI” Model, showing: a) the DFN model, and b) corresponding GSI number assigned to each cell (models are not related).. 43 Figure 4-7 10-80 Random Joint Model ........................................................................................ 46 Figure 4-8 Principal stress tensors for the 10-80 Random Joints Model ...................................... 46 Figure 4-9 Principal stress tensor for two completely random models ........................................ 47 Figure 4-10 Stress-strain curve for an unconfined 10-80 Random Joint Model .......................... 47 vii  Figure 4-11 Stress-Strain curves for 3 10-80 Random Joint Models with varying number of cells ...................................................................................................................................................... 48 Figure 4-12 10-80 Most Represented Model ................................................................................ 49 Figure 4-13 Principal stress tensors for the 10-80 Most Represented Model............................... 50 Figure 4-14 Stress-Strain curve for the unconfined 10-80 Most Represented Model .................. 50 Figure 4-15 Sample with through going fracture in triaxial compression .................................... 51 Figure 4-16 Variation of uniaxial compressive strength for a sample with the angle of inclination of the normal to the plane of weakness to the compression axis (β) ............................................ 52 Figure 4-17 Variation of uniaxial compressive strength with the angle of inclination of the normal to the plane of weakness to the compression axis (β) for four samples with different friction angle for the through going fracture ................................................................................ 53 Figure 4-18 Normalized vertical stress graph with second order trend line for 2 joint sets with different properties........................................................................................................................ 54 Figure 4-19 Variation of peak principal stress difference with the angle of inclination of the major principal plane of weakness, for the confining pressure indicated  (modified after Brady and Brown, 2006) ......................................................................................................................... 55 Figure 4-20 Normalized vertical stress for varying plane of weakness angles for two joint sets 55 Figure 4-21 10-80 Most Critically Oriented Joint Model ............................................................. 57 Figure 4-22 Principal stress tensors for the 10-80 Most Critically Oriented Joint Model ........... 57 Figure 4-23 Stress-Strain curve for the unconfined 10-80 Most Critically Oriented Joint Model 58 Figure 4-24 a) Simple DFN model with two fractures, b) representation of the model in FLAC with for the Most Critically Oriented Joint Method, and c) representation of the ideal model in FLAC based on the relative lengths of the two joints .................................................................. 58 Figure 4-25 Most Critically Oriented Joint Model with competing fractures .............................. 59 Figure 4-26 10-80 Most Critically Oriented Joint with Length Component Model..................... 61 Figure 4-27 Principal stress tensors for the 10-80 Most Critically Oriented Joint with Length Component Model ........................................................................................................................ 61 Figure 4-28 Stress-Strain curve for the unconfined 10-80 Most Critically Oriented Joint with Length Component Model ............................................................................................................ 62 Figure 5-1 Examples of the 20x20m models with fracture intensities corresponding to the:  a) 10th  percentile, b) 50th  percentile, and c) 90th  percentile............................................................ 63 Figure 5-2 10-80 FLAC models with: a) 18 cells, b) 32 cells, c) 50 cells, and d) 72 cells .......... 66 Figure 5-3 30-60 FLAC models with: a) 18 cells, b) 32 cells, c) 50 cells, and d) 72 cells .......... 66 Figure 5-4 Peak strength at different confining pressures for FLAC and ELFEN simulations for the 10-80 4x8m models ................................................................................................................ 67 Figure 5-5 Peak strength at different confining pressures for FLAC and ELFEN simulations for the 30-60 4x8m models ................................................................................................................ 67 Figure 5-6 Comparison between the 10-80 (10) and 30-60 (30) FLAC simulations ................... 68 Figure 5-7 10th percentile 20x20m models a) p10a b) p10b ......................................................... 69 Figure 5-8 50th percentile 20x20m models a) p50a b) p50b c) p50c ............................................ 69 Figure 5-9 90th percentile 20x20m models a) p90a b) p90b ......................................................... 70 Figure 5-10 20x20m p10a FLAC models: a) 100 cells, b) 225 cells, and c) 400 cells ................ 71 Figure 5-11 20x20m p10b FLAC models: a) 100 cells, b) 225 cells, and c) 400 cells ................ 71 Figure 5-12 20x20m p50a FLAC models: a) 100 cells b) 225 cells c) 400 cells ......................... 72 Figure 5-13 20x20m p50b FLAC models: a) 100 cells, b) 225 cells, and c) 400 cells ................ 72 Figure 5-14 20x20m p50c FLAC models: a) 100 cells, b) 225 cells, and c) 400 cells ................ 73 viii  Figure 5-15 20x20m p90a FLAC models: a) 100 cells, b) 225 cells, and c) 400 cells ................ 73 Figure 5-16 20x20m p90b FLAC models: a) 100 cells, b) 225 cells, and c) 400 cells ................ 74 Figure 5-17 Peak strength at different confining pressures for FLAC and ELFEN simulations for the p10a 20x20m model................................................................................................................ 76 Figure 5-18 Peak strength at different confining pressures for FLAC and ELFEN simulations for the p10b 20x20m model ............................................................................................................... 76 Figure 5-19 Peak strength at different confining pressures for FLAC and ELFEN simulations for the p50a 20x20m model................................................................................................................ 77 Figure 5-20 Peak strength at different confining pressures for FLAC and ELFEN simulations for the p50b 20x20m model ............................................................................................................... 77 Figure 5-21 Peak strength at different confining pressures for FLAC and ELFEN simulations for the p50c 20x20m model................................................................................................................ 78 Figure 5-22 Peak strength at different confining pressures for FLAC and ELFEN simulations for the p90a 20x20m model................................................................................................................ 78 Figure 5-23  Peak strength at different confining pressures for FLAC and ELFEN simulations for the p90b 20x20m model ............................................................................................................... 79 Figure 5-24 Displacement for two 20x20m FLAC ubiquitous joint models with a 60 degree fracture .......................................................................................................................................... 79 Figure 5-25 Models used for testing the scale effect (dimensions in meters) .............................. 81 Figure 5-26 Area vs. normalized UCS ......................................................................................... 82 Figure 5-27 Area vs. normalized Young’s modulus ..................................................................... 82 Figure 5-28 a) 60x60m model b) 20x20m model ......................................................................... 84 Figure 5-29 First stage of loading for the 60x60m model ............................................................ 85 Figure 5-30 Yielding in tension of the first rock bridges ............................................................. 86 Figure 5-31 Close up of the yielding in tension of the first rock bridges ..................................... 87 Figure 5-32 Fully developed shear zones going all the way through the model after peak strength was reached ................................................................................................................................... 88 Figure 5-33 20x20m model after peak strength was reached ....................................................... 89 Figure 6-1 Slope geometry and DFN boundaries (dimensions in meters) ................................... 91 Figure 6-2 Slope1 DFN ................................................................................................................ 91 Figure 6-3 Slope2 DFN ................................................................................................................ 92 Figure 6-4 Slope3 DFN ................................................................................................................ 92 Figure 6-5 Slope4 DFN ................................................................................................................ 93 Figure 6-6 Models geometry (all dimensions are in meters) ........................................................ 94 Figure 6-7 Windowed view of the Slope 1, 0.5 m mesh .............................................................. 96 Figure 6-8 Windowed view of the Slope 1, 1 m mesh ................................................................. 97 Figure 6-9 Windowed view of the Slope 1, 2 m mesh ................................................................. 97 Figure 6-10 Slope1  horizontal displacement contours, in meters, for: a) 0.5m mesh, b) 1m mesh, and c) 2m mesh, showing only minor differences ........................................................................ 98 Figure 6-11 Slope 1 horizontal stress contours, in Pa, for: a) 0.5m mesh, b) 1m mesh, and c) 2m mesh,  showing only minor differences ........................................................................................ 99 Figure 6-12 Slope1 vertical-stress contours, in Pa, for: a) 0.5m mesh, b) 1m mesh, and c) 2m mesh, showing only minor differences ....................................................................................... 100 Figure 6-13 Factor of safety analysis plasticity indicators Slope1 ............................................. 102 Figure 6-14 Factor of safety analysis plasticity indicators Slope2 ............................................. 103 Figure 6-15 Factor of safety analysis plasticity indicators Slope3 ............................................. 103 ix  Figure 6-16 Factor of safety analysis plasticity indicators Slope4 ............................................. 104 Figure 6-17 Influence of fracture orientation on horizontal displacements: a) Slope1, b) Slope2, c) Slope3, and d) Slope4. Displacements are reported in meters ............................................... 105 Figure 6-18 Plasticity indicators for a typical slope model ........................................................ 106 Figure 6-19 Slope1 horizontal displacements modeled in ELFEN. Displacements are reported in meters .......................................................................................................................................... 108 Figure 6-20 Slope1 horizontal displacements modeled in Phase². Displacements are reported in meters .......................................................................................................................................... 108 Figure 6-21 Slope1 horizontal displacements modeled in FLAC. Displacements are reported in meters .......................................................................................................................................... 109 Figure 6-22 Slope2 horizontal displacements modeled in  ELFEN. Displacements are reported in meters .......................................................................................................................................... 110 Figure 6-23 Slope2 horizontal displacements modeled in Phase². Displacements are reported in meters .......................................................................................................................................... 110 Figure 6-24 Slope2 horizontal displacements modeled in FLAC. Displacements are reported in meters .......................................................................................................................................... 111 Figure 6-25 Slope1 horizontal stresses modeled in ELFEN. Stresses are reported in Pa .......... 112 Figure 6-26 Slope1 horizontal stresses modeled in Phase². Stresses are reported in MPa ......... 113 Figure 6-27 Slope1 horizontal stresses modeled in FLAC. Stresses are reported in Pa ............. 113 Figure 6-28 Slope2 horizontal stresses modeled in ELFEN. Stresses are reported in Pa .......... 114 Figure 6-29 Slope2 horizontal stresses modeled in Phase². Stresses are reported in MPa ......... 114 Figure 6-30 Slope2 horizontal stresses modeled in FLAC. Stresses are reported in Pa ............. 115      1  1 Introduction 1.1 Problem Statement Over the past few years, research efforts have been focussed on the development of new analytical techniques to replace the use of empirical rock mass rating systems to better evaluate fractured rock mass properties and simulate its behaviour.  These efforts have been driven by the need for new modeling tools capable of simulating the types of behaviour encountered in the context of large open pits and block cave mining.  As the scale of the excavation increases, structurally controlled failures such as plane shear and wedge failure are thought to be less dominant and more complex structurally controlled failures such as step path failure and large scale toppling failure develop.  Therefore, when designing large open pit mines like Diavik, Palabora, Chuquicamata or Bingham Canyon, deeper-seated multi-bench failure must be carefully assessed.  Consequently the modeling technique(s) chosen must be able to properly capture the role of rock mass structure with respect to progressive strength reduction and the kinematic influence on the potential failure of the slope.  Furthermore, joint orientation has a significant effect on the evolution and rate of cave propagation (Sainsbury et al., 2009), and thus has a major impact on the design of a block cave mine and its profitability.  A good example of the importance of the rock mass structure is provided by the Palabora open pit/caving mine.  Preferential orientation of joints intersecting the cave volume and stress-strain interactions between the cave and overlying open pit caused the cave to deviate towards and behind the north wall of the pit leading to its failure.  This incident jeopardized the safety of critical infrastructure located behind the crest of the pit slopes and resulted in sterilisation and dilution of the ore reserve (Sainsbury et al., 2009).  Conventional continuum models based on rock mass properties derived from empirical rock mass rating systems are limited in reproducing such behaviour, at least with respect to the kinematic controls introduced by jointing and faults.  This calls their suitability into question where joints are expected to play an important role in the potential failure mode.  On the other 2  hand, continuum methods are easier to set up and run faster (i.e., are more computationally efficient), therefore allowing multiple scenarios to be modeled fairly quickly.  Discontinuum methods that explicitly model discontinuity behaviour, on the contrary, have the ability to accurately capture the effect of the rock mass structure but take more time to setup and are more computationally intensive.  If stochastic DFN modeling is being used to provide the input for the representation of the discontinuity network in the model, multiple DFN realizations will be produced requiring multiple (lengthy) discontinuum model runs to yield meaningful results.  This can be very problematic especially when modeling large scale excavations in 3-D.  A Synthetic Rock Mass (SRM) approach has been proposed by some to address the rock structure issue (Pierce et al., 2007, Mas Ivars et al., 2008, Sainsbury et al., 2009).  However, this method adds an extra modeling step, more uncertainties, and results in a loss of connection between the field and laboratory data and the strength properties used in the model.  The objective of this thesis is to investigate and develop the best means to reproduce, in a computationally efficient way, the behavior of a fractured rock mass in a continuum simulation by developing an analytical geomechanical upscaling approach using Discrete Fracture Network models (DFN) and strength properties acquired through lab testing as a starting point.  In order to reach this objective, a detailed literature review on DFN modeling, rock mass rating systems and brittle rock failure simulations has been conducted. Next, the characteristics of three different commercial codes used to model rock mass behavior (FLAC, ELFEN and Phase²) are investigated and compared. Finally, an analytical upscaling approach is developed, tested on small scale simulated rock samples and subsequently applied to a large-scale slope model.     3  1.2 Research Objectives The primary research objective of this thesis is to develop a computationally efficient geomechanical upscaling approach to reproduce the behavior of fractured rock masses in 2-D continuum simulations using DFN models and laboratory data as a starting point.  In order to achieve this goal, several secondary objectives can be identified:  1. Determine the best constitutive model to use in Itasca’s finite-difference code FLAC (Mohr-Coulomb or Ubiquitous Joint) to capture the influence of rock structure in a computationally efficient manner.  2. Develop a technique to select the angle of the weakness plans for the Ubiquitous Joint model in FLAC based on the properties of the fractures present in the rock mass studied (e.g., most represented, most critically oriented, most persistent, etc., or a combination).  3. Develop a method to derive the continuum model properties based on DFN simulations and data from laboratory testing of intact rock.  4.  Study the mesh/grid size effect for different models.  5. Verify the ability of the upscaling approach to reproduce scale effects associated with rock mass behavior.  6. Compare the results obtained with the proposed geomechanical upscaling approach in FLAC to the results of Elmo and Stead (2010) obtained using Rockfield’s FEM/DEM brittle fracture code ELFEN.  7. Investigate the limitations of modeling fractured rock masses in continuum simulations.  8. Apply the proposed geomechanical upscaling approach to large scale models and compare the results with other modeling techniques (i.e., ELFEN and Phase²). 4  1.3 Thesis Organization This thesis is divided into seven chapters:  - Chapter 1 introduces the problem and presents the research objectives.  - Chapter 2 presents the literature review on discrete fracture networks, rock mass rating systems and previous work on the simulation of brittle rock failure.  - Chapter 3 explains the distinctive characteristics of the different software used to simulate the behaviour of fractured rock masses in this research  - Chapter 4 describes how the upscaling approach was developed and the different hypotheses tested  - Chapter 5 presents the testing of the upscaling approach on two 4x8 m models, seven 20x20 m models and one 100x100 m model  - Chapter 6 presents the results of the simulation of the excavation of four slope models using the upscaling approach  - Chapter 7 summarises the conclusions of this research and provides suggestions for further research work.     5  2 Literature Review 2.1 Discrete Fracture Networks (DFN) 2.1.1 Description A DFN is a way of representing the rock mass fabric as accurately as possible by stochastically generating fractures in a 3-D volume.  This technique accounts for the spatial variability of the different parameters.  Orientation, persistence, termination and aperture are each assigned a statistical distribution and a model is chosen for the generation of the fractures in space. Since its conception in the 1970’s and early 1980’s, DFNs have been used in the investigation of slope stability, block caving and underground mining, tunneling, ground water hydrology, and geothermal and petroleum engineering. 2.1.2 Brief History Although DFNs as we know them were first developed in the 1970’s, the study of joint properties started much earlier.  At the end of the 19th century, Woodworth (1897) published an article describing joint features.  In his work, he describes fractures as ellipsoid or discoid structures (Woodworth, 1897), assumptions still used today in many DFN models (Beacher et al., 1977; Barton, 1978).  Woodworth also points out that joints generally occur in sets either parallel or slightly divergent, thus expressing the spatial variability of discontinuity orientations.  In 1965, in his PhD dissertation, Snow was probably the first to describe a model for the distribution of joints in space. The orthogonal model, as described by Snow (1965), is characterised by two or three orthogonal sets of parallel unbounded joints, with constant spacing.  In 1970, Piteau and Robertson, were one of the first to attempt to take the spatial variability of joint properties into account for slope stability problems.  Piteau (1970) suggested that before a slope stability analysis is attempted, a model of the rock mass is required in which a true and statistically-based sampling of the joint properties, including their geometrical characteristics, are  represented.   6  Snow (1968, 1970) proposed statistical distributions for joint properties.  He made the assumptions that the occurrence of joints in a borehole follows approximately a Poisson distribution and that the aperture of discontinuities follows a lognormal distribution.  He used those hypotheses to develop a method to determine discontinuity frequency from the proportion of zero discharge packer permeability tests that occur.  Bridges (1975 and 1976) demonstrated that fractures in rock can be described by a series of mathematical statistical models and discussed the presentation of fracture data for rock mechanics.  Bridges (1976) generated a joint trace plan using a computer.  The fracture centers were randomly distributed in space and the joint length followed a lognormal distribution.  Bridges (1976) had to arrange his model by hand because, according to him, in real outcrops discontinuities do not cross each other in a way that they would be arranged in a purely random model. He however suggested that such diagrams may form the basis for constructing more realistic models of fractures patterns.  Baecher et al. (1977) developed a conceptual model of joint geometry in which the center points of joints (non-ubiquitous discs in this case) are randomly and independently distributed in space forming a Poisson field, thus generating Poisson lines in a 2-D trace plan. Commonly referred to as a Baecher model, the authors suggested that the use of their model is applicable for the study of rock mass behaviour and the design of tunnel support (Baecher et al. 1977).  The following year, Barton used a similar model to confirm his findings (Barton, 1977) and reproduce the data obtain from the CSA Mine (Barton, 1978).  The joints in Barton’s 1978 model were assumed circular with a logarithmic normal diameter distribution and were randomly generated in a space inside a box.  Barton used the trace lengths generated by the model to compare them with the trace lengths from the field.  According to him, the chord length distribution derived from the modelling closely matched the sample trace length from the mine.  In 1979, Veneziano proposed a model using a combination of Poisson planes and Poisson lines to produce bounded joints on joint planes (see also Dershowitz, 1984), features previous models were incapable of reproducing. This work led to a report on risk analysis for rock slopes in open 7  pit mines under a contract from the United States Bureau of Mines (Einstein et al., 1980).  The first volume of this report reviews all joint properties and suggests, based on literature reviews and new data from 10 000 discontinuities, statistic distributions for each one of them.  The third part of the report suggests the use of stochastic models for assessing kinematic admissibility and failure mechanisms in jointed rock masses.  The authors suggested the use of the Veneziano model in their report and implemented it in two programs for slope stability assessment, JOINTSIM and SLOPESIM.  According to Einstein et al. (1980), their approach was the first to include the correct treatment of joint persistence in the calculations.  In 1984, Dershowitz created a new spatial model (Dershowitz model) and compared it with four other DFN models in his Ph.D. dissertation entitled “Rock Joint Systems” (Dershowitz, 1984).  These included the orthogonal model, the Baecher model, the Veneziano model and the mosaic tessellation model.  According to Dershowitz, the major flaw of these models was their incapacity to represent non-planar joints.  In 2002, Jing and Hudson published an article on the techniques, advances, problems and likely future development directions in numerical modelling for rock mechanics and rock engineering (Jing and Hudson, 2002).  Here they noted that a critical issue in DFN is the treatment of bias in the estimation of joints properties.  Jing and Hudson also pointed out that despite its advantages the lack of knowledge of the geometry of discontinuities limits the application of DFNs.  Thus, their adequacy highly depends on the interpretation of the natural fracture network geometry, which cannot be even moderately validated in practice (Jing and Hudson, 2002).  Recently, DFNs have been used for multiple purposes and research continues for new applications.  For instance, Pine et al. (2007) described a new discrete fracture modelling approach for rock masses that gives a greater insight into failure mechanisms in compression, tension and shear loading regimes.  Elmo et al. (2008), examined the problem of the interaction between block cave mining and a large overlying open pit by using the finite element/discrete element modeling approach (or FEM/DEM).  Additionally, Bakun-Mazor et al. (2009) examined the significance of mechanical layering in the deformation of blocky rock masses using an 8  integration of geologically-based discrete fracture models and numerical discrete element models (Bakun-Mazor et al., 2009). 2.1.3 Spatial Models According to Dershowitz et al. (1998), most fracture patterns belong to one of three mathematical models: Poisson, Geostatistical or Fractal.  Poisson models assume that fractures are randomly distributed in space implying that joint spacing follows a negative exponential distribution along a sampling line.  Geostatistical models imply that the location of a new fracture depends on the location of previous fractures; i.e., there is a higher probability for a new fracture to grow near pre-existing fractures.  This correlation diminishes with distance and finally disappears at large distances.   The Fractal models include the Box dimension, Mass dimension and Spectral models.  The Box model (Barton, 1995) describes the amount of fractured rock as a function of scale and quantifies how space-filling the fracture pattern is.  The Mass dimension describes how fracture intensity scales with the amount of rock observed. Lastly, the Spectral model describes how fracture intensity variability scales with core length, conforming to a power law (Staub et al., 2002).  Amongst the most commonly used models for fracture generation are: the Orthogonal model, the Baecher or random disc model, the Enhanced Baecher model, the nearest neighbour model and the Levy-Lee model.  Each model has very different characteristics and is used to represent different type of fracture pattern.  For example, the Beacher model is usually used for homogenous rock while the orthogonal model is used for rock masses with completely defined rectangular rock blocks.  A detailed review of models used for DFN simulations is provided in Appendix A.   9  2.1.4 Use of DFNs in this Thesis Research carried out in this thesis will make use of DFN modeling to retrieve fracture properties necessary for the implementation of the analytical geomechanical upscaling approach.  Specifically, properties such as fracture length, fracture density and dip angle will be used, coupled with data from laboratory testing on intact rock, to determine the different input properties required for the selected continuum modeling.  Table 2-1 presents the different types of measurement of fracture intensity, density and porosity used in DFN modeling.  All DFN models used for this research were provided by Elmo (2010, personal communication).   Dimension of Measurement  0 1 2 3 Dimension of Sample 1 P10 Fractures per unit length P11 Length of fractures per unit length   Linear Measures 2 P20 Fractures per unit area P21 Length of fractures perunit area P22 Area of fractures per unit area  Areal Measures 3 P20 Fractures per unit volume  P32 Area of fractures per unit volume P33 Volume of fractures per unit volume Volumetric Measures   Density  Intensity Porosity  Table 2-1 Fracture intensity, density and porosity measurements (modified after Golder Associates, 2008)    10  2.2 Rock Mass Classification and Characterization 2.2.1 Rock Mass Classification Systems  Rock mass classification involves the use of rating systems that are widely used as an engineering tool to: - provide better communication between geologists, designers, contractors and engineers; - correlate engineering observations, experience and judgment more effectively through quantitative; - provide numbers in place of descriptions for engineering calculations; and - classification helps in organizing knowledge (Singh and Goel, 1999).  In addition, several classification systems offer rapid assessment of rock mass strength that can be related directly to stability or support design (Lorig, 2007).  Rock mass rating systems rely on parameters typically weighted towards discontinuities and their properties, for example the number of joint sets, their persistence and roughness, the presence of alteration, infilling and groundwater conditions, and sometimes also the strength of the intact rock and the stress magnitude (Edelbro, 2003).  Major weaknesses of these systems include the simplification of the rock mass conditions and characteristics into a single number, they are not based on mechanics, they ignore scale effects (Lorig, 2007), they cannot describe anisotropy and time dependent behavior, and they do not consider failure mechanisms, deformation or rock support interactions (Riedmüller et al., 1999). Table 2-2 lists rock mass rating systems and Table 2-3 presents the parameters included in some of those systems.   11  Name of Classification Author and First version Country of origin Applications Remarks Rock Load Theory Terzaghi, 1946 USA Tunnels with steel support Unsuitable for modern tunneling Stand-Up Time Lauffer, 1958 Austria Tunneling Conservative New Austrian Tunneling Method (NATM) Rabecewicz, 1964/65 Austria Tunneling in incompetent (overstressed) ground Utilized in squeezing ground conditions Rock Quality Designation (RQD) Deere et al., 1966 USA Core logging, tunneling Sensitive to orientation effects Rock Classification for Rock Mechanical Purposes Patching and Coates, 1968 Canada  For input in rock mechanics  Unified Classification of Soils and Rocks Deere et al., 1969 USA Based on particles and blocks for communication  Rock Structure Rating (RSR) Wickham et al., 1972 USA Tunnels with steel support Not useful with steel fibre shotcrete Rock Mass Rating (RMR) Bieniawski, 1974 South Africa Tunnels, mines, foundations etc. Unpublished base case records Q-System Barton et al., 1974 Norway Tunnels, large chambers  Mining RMR (MRMR) Laubscher, 1977 South Africa Mining  Typological classification Matula and Holzer, 1978  For use in communication  Unified Rock Classification System (URCS) Williamson, 1980 USA For use in communication  Basic Geotechnical Description (BGD) ISRM, 1981  For general use  Rock Mass Strength (RMS) Stille et al., 1982 Sweden  Modified RMR Modified Basic RMR (MBR) Cummings et al., 1982  Mining  Simplified Rock Mass Rating (SRMR) Brook and Dharmaratne, 1985  Mines and tunnels Modified RMR and MRMR Slope Mass Rating (SMR) Romana, 1985 Spain Slopes  Slope Rock Mass Rating Robertson, 1988  Slopes Modified RMR Ramamurthy/Arora Classification (RAC)  Ramamurthy and Arora, 1993 India For intact and jointed rock Modified Deere and Miller approach Geological Strength Index (GSI) Hoek et al., 1995 Canada Mines and tunnels  Table 2-2 History of  rock mass rating systems (modified after Edelbro, 2003). Continued on next page.  12  Name of Classification Author and First version Country of origin Applications Remarks Rock Mass Number (N) and Rock Condition Rating (RCR) Goel et al., 1996 India  Stress-free Q Rock Mass index (RMi) Palmström, 1995 Norway Rock engineering, communication, characterization  Chinese SRMR (CSMR) Chen, 1995 China Slopes  Modified Rock Mass Classification (M-RMR) Unal, 1996    Index of Rock Mass Basic Quality (BQ) Lin, 1998    Notes: i) RSR was a forerunner to the RMR system, though they both give numerical ratings to the input parameters and summarize them to a total value connected to the suggested support  ii) The Unified Rock Classification System (URCS) is associated to Casagrandes classification system for soils in 1948 Table 2-2   History of rock mass rating systems (modified after Edelbro, 2003) (continued)    Classification system Parameters  RQD RSR Q MRMR RMS MBR SMR RAC GSI N Rmi Block size - - - - - x - - - - x Block building joint orientations - - - - - x - - - - x Number of joint sets - - x - x - - - - x x Joint length - - - - - - - - - - x Joint spacing x x x x x x x x x x x Joint strength - x x x x x x x x x x Rock type - x - - - - - - - - - State of stress - - - x x - x - - - - Groundwater condition - x x x x x x x - x - Strength of the intact rock - - x x x x x x x x x Blast damage - - - - - x - - x - - Table 2-3 Parameters included in different system (modified after Edelbro, 2003)  13  2.2.2 Rock Mass Characterization and its Use in this Thesis  Bieniawski’s (1989) RMR and Barton’s (1974) Q-system represent the two most commonly used classification systems used for geotechnical engineering design. These were based on experience gained in the design of shallow tunnels in sedimentary rocks in South Africa and crystalline rocks in Norway, respectively.  Both systems are based on assigning weighted ratings to rock mass parameters important to its behavior.  These are then applied to zones of similar geotechnical characteristics often defined by major structures such as faults or changes in rock type, although significant changes in discontinuity properties within the same unit can justify the division of the rock mass into a number of smaller structural sections. A detailed description of these and other rock mass classification systems are provided in Appendix B.  Whereas rock mass classification systems are empirical design tools, used directly to provide guidance on rock support, caveability, excavation stability, etc., rock mass characterization involves quantifying the parameters governing rock mass behavior. The objective of characterization is often to derive rock mass properties for use with numerical modeling (e.g. Hoek et al., 2002). Based on these definitions, the work carried out in this thesis utilizes and applies to rock mass characterization.  The Geological Strength Index (GSI) was introduced by Hoek et al. (1995) to complement the Hoek-Brown failure criterion, as a way to scale lab-derived rock strength properties (s, a and mi) to those applicable to a jointed rock mass (s, a and mb).  GSI values range from about 10 for very weak rock masses to 100 for intact rock, and is approximately based on RMR (although the proponents of the system stress that GSI should be assessed independently of RMR).  Whenever using the GSI tables, it is very important to carry out the classification on an undisturbed rock mass (Hoek and Brown, 1997).   Cai et al. (2004) proposed a quantitative approach to evaluating GSI.  Using block volume and joint condition as input parameters, they developed a relation allowing engineers to calculate GSI from physical properties.  Figure 2-1 presents the quantification of GSI chart.  14  For this thesis, GSI is used in an attempt to derive rock mass properties based on fracture properties retrieved from DFN modeling.    Figure 2-1Quantification of GSI chart and block delimited by three joint sets (Cai et al., 2004) 15  2.3 Brittle Rock Failure Simulations Numerous studies have been carried out investigating the modeling of brittle rock failure, fractured or intact, using a variety of numerical techniques.  These can be separated into two groups depending on whether stress-induced fracturing is represented implicitly or explicitly. Often these investigations start with the simulation of a laboratory uniaxial compression test.  2.3.1 Implicit Simulations Implicit simulations are those where fractures and the intact rock are not separately defined, for example by using different element types.  Instead, implicit models are composed of multiple cells forming a grid shaped to fit the object to be modeled.  Strength properties are assigned to every cell of the model based on the intact rock properties and a reduction factor that varies with the different techniques used.  Brittle fracturing is then simulated through yield relationships assigned to the cells, where groups of yielding cells are interpreted as implicitly representing the development of a fracture through localization and shearing. To date, implicit simulations have been used much more extensively to simulate the fracturing of intact rock under stress than the behavior of fractured rock masses under stress, which is the aim of this thesis. 2.3.1.1 Brittle Rock Failure Simulations Using FLAC FLAC is a 2-D explicit finite-difference program that can simulate the stress-strain behavior of structures built on or within soils or rock.  The materials are represented by elements that form a grid shaped to fit the object to be modeled.  Each one of those elements is assigned a constitutive model and associated properties that will define its behavior when stressed.  Constitutive models can involve elastic, plastic and/or time-dependent rheologies.  A more thorough description of the software is presented later in the methodology section.  Fang and Harrison (2001, 2002) used FLAC to develop a mechanical local degradation index and approach to model the brittle fracture of heterogeneous rocks as observed in intact rock samples (see Figure 2-2)  and intact mine pillars.  Their technique is based on the use of the FLAC Mohr-Coulomb constitutive model applied to a finely discretized problem domain using numerous small elements.  Each element is assigned slightly varying randomly distributed strength properties to represent the intact rock heterogeneity.  Once loaded, the weaker elements  failed first.  These then led to the yielding of eventual development of “fracturesfailed due to through-going fractures.lab results while showing realistic fracture propagation  Figure 2-2 Simulated fracture process of a rock specimen under uniaxial compression  2.3.1.2 Brittle Rock Failure Rock Failure Process Analysisnon-linear and discontinuous behavior related to rockseepage-stress-damage.  Tang evolution of damage and associated seismic events due to in brittle rock.  Three features distinguish his method from conventional numerical approach(1) by introducing parameter behavior even though using a elastic modulus reduction forthough computationally simpler failed elements, the seismicities associated with the progressive failure in rock can Tang (1997) used his model to simulate the progressive failure leadpillar and more extensively to model the brittle failure of intact rock samples (Tang1998, Tang et al. 2000a and 2000b) (see elements with varying strength properties.  When one of the elements fails, its properties are changed to weaker and straincritically stressed neighboring elements and the ”, which would grow until the sample was considered to have   This technique produced good results closely matching patterns.  (modified after Fang and Harrison, 2002a) Simulations Using RFPA2D   2D (RFPA2D) is a finite-element program cap heterogeneities under the effect of (1995, 1997) created the RFPA2D code to model the observed progressive failure leading to collapse heterogeneity into the model, RFPA2D can simulate noncomputationally simpler linear formulation; (2) by introducing  failed elements, RFPA2D can simulate discontinuum continuum mechanics; and (3) by recording the eventing to collapse of a crown Figure 2-3).  His models are composed of a grid of -softened mechanical properties.  Failed elements then connect and 16  able of simulating es: -linear behavior -rate of  be simulated.  and Kaiser, 17  form fractures that have an unpredictable geometry.  These fractures eventually grow and ultimately lead to the failure of the simulated test specimen.   Figure 2-3 Sensitivity of failure modes to local variation for five specimens with same parameters of mechanical properties final failure modes (simulated with RFPA2D) (Tang, 2000a) 2.3.2 Explicit Simulations Explicit simulations are those where fractures and intact rock are clearly defined as two distinct element types with different properties and behavior.  They allow for a more realistic representation of the rock mass but are more computationally expensive.  2.3.2.1 Brittle Rock Failure Simulations Using PFC Particle Flow Code (PFC) is a 2-D (and 3D) discontinuum code capable of analyzing the interaction of many discrete objects exhibiting large-strain and/or fracturing.  PFC models are composed of multiple small disc/sphere-shaped particles that can be grouped together into any shape.  They can be defined as either bonded or unbonded to replicate the behavior of cemented or granular materials in order to accommodate a wide variety of simulations — from rapid flow to brittle fracture of a stiff solid.  Diederichs (2003) used PFC to predict the brittle failure of intact rock under stresses.  His models treat the intact rock mass as a heterogeneous arrangement of discs bonded together by elastic springs.  The test sample is formed by creating a random assembly of particles with varying radii and inflating the particles until maximum contact density is achieved. Bonds between particles are then formed using contact stiffnesses and strengths assigned following a normal distribution function (see Figure 2-4).  The model is then loaded and bonds rupture.  When the density of ruptured bonds is high enough in a localized portion of the sample, the 18  cluster is considered to be a fracture.  The fractures then grow and lead to the rupture of the sample.  Diederichs’ (2003) experiments were able to yield results closely matching the failure of intact rock samples, and showed the importance of sample heterogeneity in the nucleation and development of fractures.     Figure 2-4 Confined compression test on a bonded disc sample. Cracks show as segments normal to ruptured bond. Right-hand image is a manually generated schematic of major rupture coalescence based on individual cracks in middle figure (Diederichs, 2003)  2.3.2.2 Brittle Rock Failure Simulations Using ELFEN ELFEN is a combined finite-element/discrete-element program (referred to as FEM/DEM) for the 2-D and 3-D modeling of jointed rock subjected to quasi-static or dynamic loading conditions.  Its ability to reproduce quasi-brittle fracture has been acquired through application to a large number of materials including metals, ceramics, concrete and rock (Pine et al., 2006).  In order to create an ELFEN model, pre-existing fractures (joints, faults, etc.) are first inserted as a fracture network embedded within the continuum finite-element mesh using discrete elements.  The simulation of stress-induced tensile fracturing, damage and associated softening in ELFEN is subsequently achieved by employing a fracture energy approach controlled by a designated constitutive fracture criterion. When the fracture criterion is met, a new fracture is inserted into the problem and the continuum mesh rediscretized through an adaptive remeshing algorithm. 19   Studies on brittle rock failure employing ELFEN include those by Cai and Kaiser (2004) to reproduce Brazilian testing, Stead et al. (2004) and Eberhardt et al. (2004) to simulate rock slope failure, Elmo et al. (2005), Pine et al. (2006), Elmo (2006) and Elmo and Stead (2010) to study mine pillars, and Vyazmensky et al. (2009) to simulate block caving.  Elmo and Stead (2010) used the approach developed by Elmo (2006) and Pine et al. (2006) integrating DFN modeling and ELFEN, to study the anisotropic effect of natural jointing on hard-rock pillar strength and failure.  Figure 2-5 shows the effect of varying angle between the applied principal stress direction and the inclination of the predefined fractures on the pillar peak strength while Figure 2-6 shows the effect on the mode of pillar failure.    Figure 2-5 Simulated anisotropic effects induced in a fractured rock mass by varying the angle between the applied principal stress direction and the inclination of the predefined fractures (Elmo and Stead, 2010)  20   Figure 2-6 Modes of failure for selected simulated pillars with 0.43m mean fracture length (Elmo and Stead, 2010)  2.3.3 Synthetic Rock Mass The Synthetic Rock Mass (SRM) approach combines the use of DFN simulations and PFC or FLAC (Ubiquitous Joint Rock Mass) to determine the rock mass properties at different scales.  These properties are then used to assign heterogeneous rock mass properties in a continuum code like FLAC3D to forecast displacements and potential failure of the rock mass considered. The SRM approach has been developed to overcome the limitations of empirical rock mass strength criteria like the Hoek-Brown criterion.  One of the major discrepancies between those criteria and real life rock masses is that they ignore scale effects (Lorig, 2007).  In fact, the larger the rock mass, the greater the number of potential failure paths, thus the lower the overall strength.  Sjöberg (1999) highlighted the importance of scale effects in his doctoral thesis on the analysis of large scale rock slopes (see Figure 2-7).  Lorig (2007) also mentions that a rock mass might appear blocky for an underground excavation but could appear disintegrated for a large slope  (see Figure 2-8).  SRM captures threalistic representation of the rock mass in numerical simulations.Figure 2-7 Scale effect on rockFigure 2-8 Illustration of how the same rock mass could appear blocky for an underground excavation, but  The development of the SRM method is fairly recmodels to estimate rock mass strengths for slope stability studies (observed that the collapse mechanism in the triaxial samples was mainly controlled by the failure of intact rock bridges in the rock mass (Lorig, 2007).ubiquitous joint properties to every cell of his FLAC model to describe the rock mass fabric and try to capture the structurally controlled softening of the rock mass (compared his ubiquitous joint model of a largeese important features and consequently allows for a more  -mass compressive strength (modified after Sjö disintegrated for a slope (Lorig, 2007) ent.  Carvalho et al. (2002)see Figure   Later, Clark (2006), arbitrarily assigned see Figure -scale laboratory compression test21  berg, 1999)   used 2-D UDEC 2-9).  They 2-10).  He  (50 m wide and 22  100 m high) to a continuum Mohr-Coulomb model using GSI to determine the rock mass properties.  The stress-strain curves are presented in Figure 2-11.  Clark concluded that the prediction of the peak strength was pretty similar for both models.  However, from a geotechnical point of view, the progressive softening response exhibited by the ubiquitous joint model provided a more realistic representation of the rock mass behavior.  Furthermore, the ubiquitous joint model can be set up with minimal data, represents the actual geological fabric and allows the specification of relevant support strategies based on the identification of the elastic limit of the rock mass and the creation of site specific design charts to provide initial estimates of rock mass properties (see Figure 2-12).   Figure 2-9 Rock mass sample: before (left) and after (right) numerical triaxial test (circles indicate tensile failure locations) (Lorig, 2007)   Figure 2-10 Stress damage in failed rock mass sample (Figure 2-11 Comparison of ubiquitous joint and σ modified after Clark, 2006)  Mohr-Coulomb rock mass behavior (modified  σ ε 23  after Clark, 2006)  Figure 2-12 Design chart to estimate: a) the rock mass mass classification system (In a similar approach, Park et al. (2004) sections extracted from a 3-D DFN to analyze a 30 xlevel of SKB’s Ãspö Hard Rock Laboratory (one, two or three joint sets to the intact rock on the behavior of the rock mass.  Park et al. (2004) observed that the post-peak response of the model changed from elastic for intact rock to perfectly plastic for the rock mass with three joint sets, without having to chanspecified constitutive model.  This would not be possible in other numerical models such as UDEC that rely on user-specified constitutive modelsassumptions on the behavior of the rock mass.  amount of damage developed in the rock mass increased with the number of joint setmost of the failures occurred along fractures failure.  However, Park et al.contacts due to irregularities of the scale of the individual PFC particles.Figure 2-13 Two-dimensional SRM model with PFC2D with bumpy joints (  strength and b) the secant modulus based on GSI rock modified after Clark, 2006)  used a series of PFC models generated from 2 30 m jointed rock mass from the 400see Figure 2-13).  They tested the effect of adding , requiring the user to make a priori Results from the PFC modeling showed that formed (mainly) through tension’s (2004) approach was limited by the bumpy nature of the joint   modified after Park et al., 2004)24  -D trace -m ge any user the s, and that -induced bond  25  Recent advances in PFC have overcome the limitation encountered by Park, through the development of a new smooth joint contact model and more rapid testing methodology (Lorig, 2007). As part of the industry-sponsored Mass Mining Technology project, Pierce et al. (2007) and Mas Ivars et al. (2008) developed what is now known as the SRM modeling technique using PFC3D to explicitly represent a DFN embedded within an intact rock matrix.  The SRM model is then loaded in three opposing directions and at a number of different scales to extract properties such as strength anisotropy and brittleness that cannot be derived from empirical techniques.  This process is illustrated in Figure 2-14.  It has been validated by comparing SRM results in terms of induced fractures, fragmentation and shape and advance rate of block caving to data available from Northparkes E26 Lift 2 mine.      Figure 2-14 Use of three-dimensional SRM in evaluating the rock response around a large block cave (Lorig, 2007)   26  2.3.4 Ubiquitous Joint Rock Mass (UJRM) Sainsbury et al. (2008) developed a SRM-Ubiquitous Joint Rock Mass (UJRM) approach to account for rock mass strength and anisotropy within large-scale FLAC3D continuum models. Their method differs from that of Pierce et al. (2007) and Mas Ivars  et al. (2008) in that it doesn’t use PFC3D but FLAC3D and ubiquitous joints to assess the rock mass properties from simulations of smaller model test samples (see Figure 2-15).  UJRM joint properties are assigned to each finite-difference grid cell according to Clark’s (2006) methodology, i.e. joint dip and direction are arbitrarily assigned to each cell according to predetermined proportions from stereonet data.  Sainsbury et al. (2008) applied the SRM-UJRM method to a back analysis of caving and the associated pit slope failure mechanism at the Palabora mine in South Africa to validate their method (see Figure 2-16). After calibration with the SRM sample, the matrix Young’s modulus was decreased by 50 to 70% from the intact rock properties, matrix friction, cohesion and tension were reduce to 80% of the laboratory UCS values, joint cohesion was assumed to be between 0.1 and 1% of the matrix cohesion, joint friction angles were set for each lithology in accordance with the SRM testing, and joint tensile strength was set to zero.    Figure 2-15 Use of three-dimensional SRM in evaluating the rock response around a large block cave (after Sainsbury et al., 2008)  27   Figure 2-16 The Palabora pit slope failure mechanism reproduced by the SRM–UJRM approach for representing jointed rock masses (after Sainsbury et al., 2008)  The SRM-UJRM method produces simulations that more closely match the behavior of a real rock mass than those produced using empirical rock mass strength criteria applied through elasto-plastic constitutive laws.  Its ability to reproduce scale effects and capture anisotropy and heterogeneity present in the rock mass fabric highlight its potential.  Present limitations include: limited user experience; the approach requires an extra step in the modeling process, and hence more time; the use of more than one software code is required; and the loss of connection between the final model and the field and laboratory data.  The latter implies that every time new data are acquired the SRM model has to be rerun in addition to the FLAC3D model, and more importantly, the calibration of the FLAC model is made more difficult because the exact factors that influence the reduction of the different strength properties is unknown.  The data that goes in the SRM model are measured physical properties, however, the data that comes out is the result of an unequal blending of the input properties.   28  2.4 Summary In rock engineering, numerous methods and procedures have been developed to “best” represent different rock mass properties and the heterogeneity and anisotropy often encountered in its fabric. DFN models have been developed to more realistically represent the spatial distribution of the fractures in the rock mass.  Geometrical properties such as joint persistence and dip angle are assigned statistical distributions.  Then, using the appropriate spatial model, fractures are generated in a 3-D space through a stochastic process.  This technique allows for a realistic representation of the rock mass fabric and for properties such as joint intensity to be derived from the simulations.  Rock mass characterization has evolved as a means to approximate rock mass shear strength properties through the adjustment of laboratory-derived intact rock properties using different empirical relations.  Although they have been widely used and generally yield reasonable results when used in numerical simulations, they incorporate many simplifications and assumptions that translate into deficiencies in properly capturing the influence of rock mass fabric and the behavior of a fractured rock mass under stress.   This has led to efforts to adapt different numerical techniques to simulate brittle rock failure.  These include discontinuum codes like PFC and ELFEN that allow for explicit modeling of brittle fracture initiation and failure.  It should be emphasized that even though these codes are capable of more accurately capturing the influence of rock mass fabric and brittle fracture on rock mass failure, it comes at the expense of significantly increased model run times. This has resulted in their application being largely restricted to 2-D, and as being considered impractical in the simulation of 3-D problems.     29  3 Methodology This section comprises a brief description of the three numerical modeling codes used for this research: FLAC, ELFEN and Phase². These codes were used to examine different aspects of the upscaling problem. Each code was selected based on the different ways they are able to represent and model the behavior of a jointed rock mass, with comparisons being carried out between the use of FLAC to capture the influence of rock structure in a computationally efficient manner, and brittle fracture simulations using the FEM/DEM approach in ELFEN and joint elements in Phase². 3.1 FLAC As previously described, FLAC (Fast Lagrangian Analysis of Continua) is a 2-D explicit finite difference program that can simulate the stress-strain behavior of soils and rock in response to geotechnical and engineering activities. The formulation is based on the treatment of the problem domain as a continuum that responds in accordance to one or more constitutive relationships selected by the user. Release version 6.0 was used. Figure 3-1 presents the general solution procedure as suggested by Itasca Consulting Group Inc. (2008a). 3.1.1 Finite Difference Grid In FLAC, the materials are represented by interconnected elements that form a grid shaped to fit the object to be modeled.  Elements are assigned a constitutive model and related properties that will define their behavior when stressed. This is usually done for groups of elements representing different geological units, but can also be done on an element by element basis for highly heterogeneous materials.  Grid cells are first generated as squares, but they can be later distorted so that the boundaries fit an irregular given shape.  The zones can also vary in size across the grid.  For any numerical simulations, the accuracy of the models varies with the element size.  As a general rule, finer meshes yield more accurate results, especially when localization and yield are important components of the modeled stress-strain response.  Moreover, the shape of the elements also has an influence on accuracy.  The best accuracy is achieved for elements with uniform dimensions (i.e. square zones).  Any aspect ratio (one dimension relative to the other)  above 5:1 is potentially inaccurate.  use only square elements.  In order to the mesh size effect has been carried out for all simulations run for this research. Figure 3-1 FLAC® General solution procedure (after Itasca Consulting Group, 2008a) For this reason, all simulations carried out for this research investigate the sensitivity of the results 30 obtained a study of    31  3.1.2 Constitutive Models Constitutive models act to describe, in terms of phenomenological laws, the stress-strain behavior of a material in terms of a collective behavior within a continuum. Constitutive models can be arranged into three groups: elastic, plastic and time-dependent.  The key constitutive models used for this thesis are plasticity models represented through a Mohr-Coulomb elasto-plastic shear model and a ubiquitous joint elasto-plastic model.  Plastic models, as opposed to elastic models, involve some degree of permanent deformation (failure) as a consequence of the nonlinearity of the stress-strain relations (Itasca Consulting Group, 2008c).  They also have the ability to develop localization of shear bands from their starting point as a continuum.  The latter is a crucial characteristic when modeling the brittle failure of rock. 3.1.2.1 Mohr-Coulomb The Mohr-Coulomb constitutive model describes the plastic yielding of brittle materials in shear according to the well-known Mohr-Coulomb failure criterion. Mohr-Coulomb strength is described by a cohesion component and a friction component, the latter being a linear function of the normal stress.  It is usually favored for materials with a compressive strength far exceeding the tensile strength, as in rock, soils or concrete.  It is one of the most used failure criteria in geotechnical engineering.  The key to its use here is that the development of plastic yielding involves some degree of permanent, path-dependent deformation. Once an element has reached it’s yield state, further increases in stress must be supported by neighboring elements, which in turn may yield, setting off a chain reaction leading to localization and large strain deformations. Because the representation is that of a continuum, catastrophic failure is not explicit (i.e. the grid does not shear apart) but must be interpreted using a combination of plasticity indicators, displacements and grid point velocities. More details about the Mohr-Coulomb criterion and how it is applied in FLAC can be found in the programs user guide (Itasca Consulting Group Inc., 2008c). 3.1.2.2 Ubiquitous Joint Rock Mass (UJRM) The Ubiquitous Joint Rock Mass (UJRM) constitutive model is similar to the Mohr-Coulomb model with the exception that it accounts for the presence of a plane of weakness that is ubiquitous throughout the continuum domain (Itasca Consulting Group, 2008c).  Stresses are 32  calculated for each grid zone, as well as resolved along the plane of weakness (inputted with a user-specified dip angle), to determine if failure of the element occurs either along the ubiquitous joint or through the solid, or both.  Each cell can be assigned different properties and different planes of weakness dip angles. Consequently, the lengths and spacing of the rock mass discontinuities represented by the ubiquitous joints are equivalent to the sizes of the grid zone elements.     The UJRM model has been successfully used by Clark (2006) and Sainsbury et al. (2008) to capture, to some extent, weakening and kinematic aspects of the rock mass fabric.  More details about the Ubiquitous Joint Rock Mass model and how it is applied in FLAC can be found in the programs user manual (Itasca Consulting Group Inc., 2008c). 3.1.3 Boundary Conditions As it is impossible to represent the full spatial extent of the rock mass in any numerical simulation, artificial boundaries must be added in addition to the engineering boundaries (e.g. excavation boundaries) and ground surface (if applicable).  Mechanical artificial boundaries in FLAC can be assigned either as a prescribed displacement (velocity) or stress.    Displacement boundaries cannot be controlled directly in FLAC and, in fact, they do not play a part in the calculation process (Itasca Consulting Group Inc., 2008a).  To generate a displacement boundary, grid points are simply assigned a velocity value over a certain number of time steps.  This number can be positive, negative or zero.  Once the model starts time stepping, the initial velocity is retained no matter what forces act on the grid points. To generate a prescribed-stress boundary in FLAC, forces may be applied to part of or the entire boundary.  Applied stresses may be altered at any point during the simulation.  For the modeling carried out in this thesis, in order to simulate triaxial tests, the external boundaries of the models were assigned prescribed stresses (confinement) on either side of the model and prescribed displacements for the top and bottom boundaries (axial loading).  Subsequent models analyzing slope failure mechanisms employ zero displacement boundaries 33  (normal to the boundaries) assigned to both sides and the bottom of the model.  The top of the slope model corresponds to the ground surface and therefore no artificial boundary is needed. 3.2 ELFEN ELFEN is a combined finite-/discrete-element program, or FEMDEM, for the 2-D and 3-D modeling of jointed rock subjected to quasi-static or dynamic loading conditions.  Its ability to reproduce quasi-brittle fracture has been acquired through application to a large number of materials including metals, ceramics, concrete and rock (Pine et al., 2006).  ELFEN uses fracture mechanics integrated within the FEMDEM formulation to realistically reproduce the behavior of fractured rock.  In order to create a model, a discrete fracture network is constructed and then embedded within the continuum problem domain.  The model is then discretized and meshed.  During time stepping, fracture insertion and adaptive remeshing algorithms allow existing fractures to grow and intact rock to crack if critically stressed.  Each joint in the model, existing or new, is assigned Mohr-Coulomb shear strength parameters.  A variety of constitutive models are available to describe the rock mass material including the Rotating Crack and Rankine tensile smeared crack criteria (Vyazmensky et al., 2009).  All ELFEN models used for this research were carried out by Elmo (2010, personal communication).  Outcomes were used to provide a frame of reference to analyze the simulations carried out using the geomechanical upscaling approach in FLAC. 3.3 Phase² Phase² is a 2-D elasto-plastic finite element program used to solve a wide range of mining, geotechnical and civil engineering problems (Rocscience, 2010).  It differs from FLAC in its functionality as it incorporates an automatic mesh generator.  The elements can be specified by the user as either triangular or quadrilateral and the mesh type can be set to graded uniform or radial to fit the needs of the model.    34  Like FLAC, Phase² includes a range of elastic and plastic constitutive models to define the rock mass material behavior including: Mohr-Coulomb, Hoek-Brown, Cam-Clay and others.  To be consistent with the FLAC simulations performed, all Phase² simulations carried out in this thesis use a Mohr-Coulomb constitutive model.  In addition to the artificial and real boundaries, joint boundaries can be defined in Phase².  These joint boundaries can either extend across the entire model or truncate in the middle of the continuum representation of the rock mass domain.  Joint ends can be specified as either open or closed.  A closed joint means that the end of the joint boundary is represented by only one node in the finite element mesh, and hence sliding or opening cannot take place at the joint end.  An open joint means that the joint boundary is represented by two ends thus allowing relative movement at the joint end.  If joints terminate within the rock mass it is suggested to specify the end condition as closed.  An open joint in the rock mass could lead to unexpected gaps or overlapping of material which is not kinematically feasible and thus lead to inaccuracies in the model.  On the contrary, joints ends located on external boundaries or on other joints should be defined as open joints to allow movement to occur and better capture their effect on the rock structure.  For the sake of this thesis, the fractures were imported form a DXF file and end conditions were specified as closed if terminating within the rock mass and open if terminating on another joint or at the excavation boundary.  Figure 3-2 presents a close-up example of the Phase² joint representation for one of the models constructed (and described in detail later in the thesis).  Empty circles represent open joint ends while circles with a black dot represent closed joint ends.  Like the ELFEN models, the Phase² models were carried out to provide a frame of reference to analyze the implementation of the upscaling approach.  35   Figure 3-2  Close-up of Phase2 model representation of discrete fracture network   4 Upscaling Approach Development4.1 Methodology 4.1.1 Programming To facilitate and accelerate the process of sampling the fracture properties different regions of the models, a C++ program designed to find any fracture going through a cell of an imaginary grid DFN model.  It then selects, based on different each cell of a FLAC model.  Finally, the Cinformation necessary to run the model.  developped.INPUT:-Fracture file-Model dimensions-Grid size-Intact Rock Properties-Fracture Properties-Name of file to createOUTPUT:-FLAC .dat file for the model to assign towas created (see Appendix C).  The program is superimposed criteria, the properties that will be input++ program generates a FLAC .dat file with all the Figure 4-1 presents a flow chart for the C++ code Figure 4-1 Flow chart for the C++ code Choice of the method to use (raster, most represented etc.)Determine which fracture goes through each cell (cell by cell) and the length of the fracture in that cellDetermines which fracture will represent the cell for the FLAC model based on the criteria associated with the different methods36  the over a ted into  Calculate P21, P00, P11, and P01 for each cell37  4.1.2 4x8 Models (Initial Testing) In order to develop and test different approaches in FLAC for the upscaling approach, several rock mass models have been tested in FLAC under uniaxial and triaxial compression.  To begin with, two 2-D 4x8 m DFN models of hard rock pillars, previously used in a similar study carried out with ELFEN by Elmo and Stead (2008),  were used to test the different ubiquitous joint representation criteria (see Figure 4-2).  Both models share the same intact rock and fracture properties (see Table 4-1).  They also have the same joint pattern with the exception that the joints in the second model (Figure 4-2b) are rotated by 20º from 10 and 80º to 30 and 60º.  These models were chosen because, from a rock mass classification point of view, they have equivalent rating (Q, RMR or GSI) but their behavior when loaded is expected to be very different.  Outcomes from the different FLAC 4x8 models were compared to the results obtained by Elmo and Stead (2008) with ELFEN. Note, this is a relative comparison involving ultimate pillar strength in ELFEN and mode of failure.   Figure 4-2 2D 4x8m DFN models a) 10-80 fracture pattern b) 30-60 fracture pattern   a) b) 38  Intact Rock Bulk Modulus (MPa) 40000 Shear Modulus (MPa) 24000 Cohesion (MPa) 20 Tensile Strength (MPa) 10 UCS (MPa) 100 φ (°) 45 Dilation angle (°) 5 Density (kg/m³) 3100 Fractures Cohesion (MPa) 0 Tensile Strength (MPa) 0 φ (°) 35 Dilation angle (°) 0 Table 4-1 Intact rock and fracture properties for the two 4x8m models  4.2 Representing the Fractured Rock Mass with Different Constitutive Models 4.2.1 Mohr-Coulomb The Mohr-Coulomb criterion has previously been successfully used by Fang and Harrison (2002) and Tang (1998, 2000a, 2000b) to simulate the brittle failure of rock samples under uniaxial and triaxial compression.  This constitutive model is thus a good candidate as a starting point for the development of the upscaling approach.  Two types of models have been developed, using the Mohr Coulomb model for this thesis.  They are the Raster Models and the Reduced Strength “GSI” Models.  Both types of models generated different degrees of success after being tested. 39  4.2.1.1 Raster Models The Raster Models are similar to the models used by Fang and Harrison (2002) and Tang (1998, 2000a, 2000b).  Those models are composed of many small cells that can be assigned different properties. The major advantage of these models is that, due to the small size of the mesh, they have a very good ability to accurately reproduce the fracturing process of intact rock.    Raster Models produced for this thesis, in contrast to the models used by Fang and Harrison and Tang, are not intact rock samples but include a predefined fracture pattern prior to loading.  Consequently, cells were either assigned intact rock properties or fracture properties from the beginning (see Table 4-2) depending on whether they were located over a mapped fracture or contained within the intact rock matrix.  Figure 4-3  shows the 10-80 and 30-60 Raster Models in FLAC comprised of 7200 cells.  The pink color represents cells with fracture properties whereas the red color represents cells with intact rock properties.  It can be seen that by using this technique, every single natural fracture is explicitly represented in the model. Intact Rock Bulk Modulus (MPa) 40000 Shear Modulus (MPa) 24000 Cohesion (MPa) 20 Tensile Strength (MPa) 10 φ (°) 45 Dilation angle (°) 5 Density (kg/m³) 3100 Fractures Bulk Modulus (MPa) 40000 Shear Modulus (MPa) 24000 Cohesion (MPa) 0 Tensile Strength (MPa) 0 φ (°) 35 Dilation angle (°) 0 Density (kg/m³) 3100 Table 4-2 Intact rock and fracture properties for the Raster Models 40    a)         b)      Figure 4-3 Raster Models (7200 cells) and DFN they are based on: a) 10-80 b) 30-60   41  Different mesh size, fracture properties and loading patterns have been tested for the Raster Models.  Peak strengths similar to those predicted by ELFEN were obtained for models with 7200 cells (medium to good correlation with ELFEN results) (see Figure 4-4).  Both models also exhibit realistic deformation patterns with maximum displacements located along or near the longer and more critically oriented fractures (see Figure 4-5).  It is anticipated that a better comparison between FLAC and ELFEN results for the 30-60 model could be achieved by further studying the effects of loading rate and grid size in FLAC.  However, it is argued that objective of the current study was purely to test the feasibility of the raster approach.   The goal of this thesis is to find a computationally efficient means to represent the discontinuous rock mass fabric as a continuum.  Notwithstanding the good correlation between the FLAC and ELFEN result, it is argued that the Raster Model technique is only effective for a very small mesh size, it is consequently not practical for full scale slope models and other larger simulations. Raster Models are therefore not a good candidate for the proposed upscaling approach.  Figure 4-4 Raster Models (7200 cells) compared to ELFEN models, where Sigma1 and Sigma3 refer to the major and minor principal stresses, respectively 42                Figure 4-5 Raster Models’ deformation after peak strength a) 10-80 b) 30-60  4.2.1.2 Reduced Strength “GSI” Models Reduced Strength “GSI” Models are based on the principle that a cell with a higher joint intensity should have lower strength properties.  Figure 4-6 illustrates how the Reduced Strength model works.  Like for the Raster Model, a search grid is superimposed over the DFN model.  However, this time the FLAC mesh is coarser (requiring less computational effort) and the properties that will be input are not a function of a given fracture (or fractures) being physically included within a grid cell, rather each grid cell will be assigned given properties based on the local (grid-based) fracture intensity and principal orientation.  Based on these considerations, a GSI value is assigned to each cell and the properties for the continuum are determined via the equations developed by Hoek et al. (2002).  This technique has the advantage over the Raster Model method that bigger elements can be used, thus making it suitable for large-scale models.  A literature review has been conducted to find a suitable relation between fracture intensity and GSI that could later be modified to incorporate dip angle.  Although an interesting article by Cai et al. (2003) links GSI with joint spacing and joint roughness, the search has proven unsuccessful for an automatable relation that could be used to correlate joint intensity and GSI. a) b) 43   Independently of the assumed correlation between GSI and fracture intensity within a grid cell, it is argued that the Reduce Strength model is likely to be grid size dependent, the coarser the grid size, the less the model will be able to capture anisotropic effects, since for relatively large grid sizes, it is believed this method would not differ from a more traditional equivalent continuum approach.  For these reasons, the Reduced Strength “GSI” Model technique was deemed not suitable for the development of the upscaling approach.                      Figure 4-6 Schematic representation of the Reduced Strength “GSI” Model, showing: a) the DFN model, and b) corresponding GSI number assigned to each cell (models are not related)  4.2.2 Ubiquitous Joint Rock Mass (UJRM) Building on the experience acquired from the Mohr-Coulomb models, various UJRM models were next tested.  The advantages of the UJRM constitutive model over the Mohr-Coulomb model are that: i) the weak planes embedded in the cells allows for directionality to be captured without having to use a really fine mesh, ii) properties for the models can be directly derived from field and laboratory data without the need of empirical rock mass rating systems, and iii) it can be easily automated to facilitate its application to large-scale models.  a) b) 44  Four types of UJRM models have been tested to develop the upscaling approach based on four different hypotheses on how to choose the angle of the ubiquitous joint in each cell.  These are the “Random Joint”, the “Most Represented”, the “Most Critically Oriented Joint” and the “Most Critically Oriented Joint with Length Component”.  These four techniques are explained in the following sections and applied to the simulation of the mechanical behavior of the same 4x8m pillars shown in Figure 4-2.  Table 4-3 presents the properties used for the generation of the 4x8m UJRM models.  For all four methods developed here, any cell not containing a fracture is assigned intact rock properties and is considered to represent an intact rock bridge.  Intact Rock Bulk Modulus (MPa) 40000 Shear Modulus (MPa) 24000 Cohesion (MPa) 20 Tensile Strength (MPa) 10 φ (°) 45 Dilation angle (°) 5 Density (kg/m³) 3100 Fractures Cohesion (MPa) 0 Tensile Strength (MPa) 0 φ (°) 35 Dilation angle (°) 0 Table 4-3 Intact rock and fracture properties for the 4x8m UJRM models 4.2.2.1 Random Joint The Random Joint technique is similar to the approach adopted by Clark (2006) and Sainsbury et al. (2008), with the difference that the fracture angle for each cell is randomly picked from the orientations of fractures from the DFN that fall within a particular cell and not from the distribution of angles for the DFN as a whole.  This procedure has the advantage over Clark’s method of respecting the concentration of certain joint orientations in particular areas of the DFN model.  The Random Joint technique also has the benefit of accurately representing the proportionality of every joint orientation present in the DFN.   45   Figure 4-7 presents joint angles for a 5 x 10 cell 10-80 Random Joint Model.  As the joints are randomly picked for each cell, every realization yields different results thus increasing the uncertainty for this type of model.  Figure 4-8 shows the principal stress tensors plotted as crosses aligned with the major and minor principal stress directions and scaled with respect to magnitudes.  Due to the random nature of the joint angles assigned, the stress is mostly concentrated in a single row and does not vary vertically.  This behavior is characteristic of models with no particular joint pattern.  As a comparison, Figure 4-9 shows the principal stress tensor for two completely random models with joint angles picked from a uniform distribution with a 0-90° range.  The two stress tensors from Figure 4-9 exhibit characteristics similar to the stress tensor of the Random Joint Model.  Those similarities thus infer that the latter fails to capture the directionality characteristics of the rock mass fabric.  The Random Joint Model suffers from its quality of being able to accurately represent the proportionality of every joint orientation present in the DFN.  The multiple joint orientations compete with each other and a coherent failure path following critical fracture orientations cannot develop.  A consequence of this phenomenon is the artificial strengthening of the sample.  Figure 4-10 shows the stress-strain curve for an unconstrained 10-80 Random Joint Model together with the expected peak strength value derived using ELFEN.  The peak strength of the 10-80 simulation with the Random Joint technique is more than twice that simulated by ELFEN (poor correlation with ELFEN results).  Increasing the number of cells in the model only further accentuates the strengthening effect of the Random Joint Model due in part to the increase of the competing effect between the different cells with different joint orientations (see Figure 4-11).  In addition, the post-peak behavior gradually changes from brittle to almost plastic when increasing the number of elements for the same model with the same loading conditions.  For these reasons, the Random Joint Model can be dismissed for the upscaling approach.      Figure 4-8 Principal stress tensors for the 10  Figure 4-7 10-80 Random Joint Model  -80 Random Joints Model46   Figure 4-9   Figure 4-10 Stress a)           Principal stress tensor for two completely random models-strain curve for an unconfined 10-80 Random Joint Model b) 47     48   Figure 4-11 Stress-Strain curves for 3 10-80 Random Joint Models with varying number of cells  4.2.2.2 Most Represented Due to the shortcomings of the Random Model, a procedure to choose the appropriate orientation of the weak plane based on fracture properties had to be developed.  This second type of model tested is termed the “Most Represented”.  For these models, the orientation of the ubiquitous joint is chosen by adding the length of fractures from the sets crossing through and/or terminates into a single cell and then comparing the total lengths obtained for each set.  The weak plane orientation is selected as being the average angle of the most represented family.  This method has the advantage over the Random Joint Model of being based on fracture properties and hence having only a unique solution for a given DFN and a given grid size.  Figure 4-12 presents the ubiquitous joint orientation for the 5 x 10 cell 10-80 Most Represented model.  The principal stress tensors show an almost uniquely vertical orientation with an intensity only varying from the model boundary to its center, except for the mid-bottom portion of the model (see Figure 4-13).  Although it seems like the weak plane orientations have slightly 49  more influence than the Random Joint Model on the distribution of stresses in the model, it is argued that the Most Represented Model still does not capture the directionality of the rock mass fabric for the given DFN pattern.  One of the weaknesses of the Most Represented Model is that it is not based on mechanics, but on quantity.  This problem can be illustrated by the following example. Using the Most Represented procedure, a cell containing two sub-horizontal fractures 2 m in length and one 3 m fracture angled at 30° (with respect to the loading direction) would be represented by a single sub-horizontal ubiquitous joint assigned to the cell in question. This flat weak plan would have little to no influence on the strength of the model.  On the other hand, the single 30° fracture would have an important effect on the strength of the sample, but is ignored by the selection procedure.  The limitations of the Most Represented Model for the proposed upscaling approach is further illustrated in Figure 4-14 showing a peak strength almost twice as high as the one simulated in ELFEN for the 10-80 Most Represented Model (poor correlation with ELFEN results).  This discrepancy is mostly attributable to the quantity based selection of the ubiquitous joint angle.    Figure 4-12 10-80 Most Represented Model   Figure 4-13 Principal stress tensors for the 10Figure 4-14 Stress -80 Most Represented Model -Strain curve for the unconfined 10-80 Most Represented Model 50    51  4.2.2.3 Most Critically Oriented Joint Building on the lessons learned from the Most Represented Model, a mechanistic-based approach to the selection of the ubiquitous joint angle was developed.  The concept at the base of the Most Critically Oriented Joint method is to identify, in every cell, the joint that will have the most influence on the strength of the sample.  Several procedures have been tested, but the one that had the most success involved the use of the theoretical uniaxial compressive strength graph for a jointed rock sample (see Figure 4-15 and Figure 4-16).  Figure 4-16 presents the theoretical uniaxial compressive strength for a sample with a single through going fracture with a friction angle (Φ) of 30°.  The minimum UCS value is achieved for an angle of 60° while the maximum value is obtained for angles between 0° and (Φ+5)° and between 85° and 90°.  Similar graphs can be plotted for any values of cohesion, tension and friction angle to fit the field data.  Figure 4-17 illustrates the effect of changing the angle of friction of the through going fracture on the uniaxial compressive strength.  The minimum UCS values are invariably obtained for values of Φ+30°.        Figure 4-15 Sample with through going fracture in triaxial compression   σ3 σ1 σ3 σ1 α β 52   Figure 4-16 Variation of uniaxial compressive strength for a sample with the angle of inclination of the normal to the plane of weakness to the compression axis (β) 53   Figure 4-17 Variation of uniaxial compressive strength with the angle of inclination of the normal to the plane of weakness to the compression axis (β) for four samples with different friction angle for the through going fracture   In order to use these graphs for the upscaling approach, the vertical stress values were normalized so the lowest value is equal to one and a second-order polynomial trend line was fitted to the parabolic part of the curve.  Figure 4-18 shows a normalized graph with fitted second-order polynomial trend lines for a model with two joint sets (A and B) with different properties.  Since parts of the theoretical curve are flat, in order to evaluate and rank the strength of fractures angled between 0° and (Φ+5)° and between 85° and 90°, the parabola is extended for values of β  higher than 85° and a slope with an initial value equal to the normalized vertical stress for β=90° is added for values lower than Φ+5° (see Figure 4-20).  Because fractures with a β angle of 0° or 54  90° should theoretically have no influence on the strength of the sample, they are assigned the same value, i.e. the highest value for the normalized vertical stress on the chart.  These modifications of the theoretical graphs are in accordance with observations made by Donath (1976) and McLamore and Gray (1967) on phyllite, slate and shale samples (see Figure 4-19).  The trwo graphs exhibit similar peak strength values for α=0° and α =90° (α=90- β), and the second graph present a gentler slope for α values on the right-hand side of the minimum (left hand-side of the graph if using β).  Figure 4-20 shows the curves that would be used to determine the most influential fracture, using the Most Critically Oriented Joint method, in a cell for a model with two joint sets with different properties.   Figure 4-18 Normalized vertical stress graph with second order trend line for 2 joint sets with different properties   Figure 4-19 Variation of peak principal stress diffeplane of weakness, for the confining preFigure 4-20 Normalized vertical stress  σ1-σ3  rence with the angle of inclination of the major principal ssure indicated  (modified after Brady and Brown, 2006)  for varying plane of weakness angles for two joint sets α 55     56  Once the graphs are constructed, the procedure for the Most Critically Oriented Joint method is is as follows:  i) Every fracture crossing through and/or terminating within a particular cell (if the length of the fracture in the cell is longer than 1/10 the length of the side of the cell) is assigned a normalized vertical stress value based on its properties and its β angle.  ii) The numbers are then compared and the fracture with the lowest value is chosen to represent the cell.  Figure 4-21 shows the plane of weakness orientation for the 5 x 10 cell 10-80 Most Critically Oriented Joint Model.  The principal stress tensor for this simulation shows a stress concentration typical of fractured rock (see Figure 4-22) showing that the mechanistic approach better captures the rock mass fabric than previous approaches.  The Most Critically Oriented Joint approach however also has limitations.  The technique is sensitive to mesh size and predicts lower peak strength than expected (medium correlation with ELFEN results) (see Figure 4-23).  This behavior is attributable to the fact that the choice of the joint angle is solely based on the friction of the joint.  Figure 4-24 shows a simple DFN model with two fractures and two representations of this model in FLAC.  The first representation is obtained with the Most Critically Oriented Joint method.  In this example, the short fracture is more critically oriented and is thus picked as the most influential joint in every cell it touches despite the fact that the longer joint clearly has, in this case, a greater influence on the sample’s strength.  The third graph in Figure 4-24 shows a better representation of what the model should look like in FLAC to truly capture the essence of the rock mass fabric.  In this representation, the bigger joint is preferred over the smaller joint because of its persistence.  Another problem with the Most Critically Oriented Joint approach is that it has the potential to create models where the fractures would compete with each other just like in the random model.  Because the joints are picked based on a cell by cell process, a model like the one presented in Figure 4-25 could be generated for a DFN with a less organized joint pattern.     Figure  Figure 4-22 Principal stress tensors for the 10  4-21 10-80 Most Critically Oriented Joint Model  -80 Most Critically Oriented Joint Model57  58   Figure 4-23 Stress-Strain curve for the unconfined 10-80 Most Critically Oriented Joint Model   Figure 4-24 a) Simple DFN model with two fractures, b) representation of the model in FLAC with for the Most Critically Oriented Joint Method, and c) representation of the ideal model in FLAC based on the relative lengths of the two joints  a) b) c) 59   Figure 4-25 Most Critically Oriented Joint Model with competing fractures  4.2.2.4 Most Critically Oriented Joint with Length Component In order to overcome the limitations of the Most Critically Oriented Joint approach, the length of the fracture has been added as a component of the equation used to pick the fracture that will represent the cell.  Consequently, the new equation becomes:  g1836g1867g1861g1866g1872 g1838g1857g1866g1859g1872ℎg1840g1867g1870g1865g1853g1864g1861g1878g1857g1856 g1848g1857g1870g1872g1861g1855g1853g1864 g1845g1872g1870g1857g1871g1871  In this case, the ubiquitous joint angle chosen to represent the cell is the one corresponding to the fracture with the highest value calculated from the above relationship.    60  The addition of the joint length component to the equation is supported by four major arguments:  i) As seen in Figure 4-3 the greatest displacements occur on fractures that are both critically oriented and long.    ii) Adding a length component to the equation takes care of problems similar to the one illustrated in Figure 4-24 where a small critically oriented joint is preferred over a longer, less critically oriented but more influential joint.  iii) It ensures continuity and helps avoid situations like that shown in Figure 4-25 where weak planes compete with one another.  In fact, adding the length component in the equation has the effect of putting the fractures in the model in the same frame of reference.  The Most Critically Oriented Joint with Length Component therefore is a model-based method instead of a cell-based method similar to the previous ones.  iv) It reduces the sensitivity of the models to mesh size effects because it allows less variation due to the fact that the fractures are ranked at the model scale and not at the cell scale.  Figure 4-26 shows the 10-80 model for the Most Critically Oriented Joint with Length Component approach.    The principal stress tensor for this simulation shows stress concentration similar to the Most Critically Oriented Joint Model hence showing that the method captures at least some aspects of the rock mass fabric (see Figure 4-27).  Figure 4-28 presents the stress-strain curve for the Most Critically Oriented Joint with Length Component.  There is a good agreement between the peak strength simulated in FLAC and ELFEN (good correlation with ELFEN results).   Figure 4-26 10-80 Most Critically Oriented Joint with Length Component Model Figure 4-27 Principal stress tensors for the 10  -80 Most Critically Oriented Joint with Length Component Model61   62   Figure 4-28 Stress-Strain curve for the unconfined 10-80 Most Critically Oriented Joint with Length Component Model  Based on the preceding discussion, the Most Critically Oriented Joint with Length Component has been chosen for the upscaling approach.  The following section presents the results of tests that have been performed to further investigate the advantages and limitations of the chosen procedure.   63  5 Validation Testing of the Upscaling Procedure In order to test the validity and applicability of the Most Critically Oriented Joint with Length Component upscaling approach, the procedure was applied to an 8 m high and 4 m wide pillar model for a range of mesh sizes and loading conditions.  Results are compared to outcomes for the same models run in ELFEN.  Because ELFEN has its own limitations and uncertainties, the objective of the comparison is not to obtain the same exact values for the peak strength, but rather to see if they produce similar trends in the same range.    These tests are next followed by further testing on several 20x20 m models (Figure 5-1).  Again, all models share the same fracture and intact rock properties (see Table 5-1).  However, this time, each model is based on a different realization of the same DFN with varying fracture intensities. Two models have a fracture intensity corresponding to the 10th percentile, three to the 50th percentile and two to the 90th percentile of the given P32 fracture intensity distribution (Elmo, 2010 personal communication).  Outcomes from the different 20x20m FLAC models were compared to the results obtained by Elmo (2010, personal communication).           Figure 5-1 Examples of the 20x20m models with fracture intensities corresponding to the:  a) 10th  percentile, b) 50th  percentile, and c) 90th  percentile   a) b) c) 64   Intact Rock Bulk Modulus (MPa) 30400 Shear Modulus (MPa) 20900 Cohesion (MPa) 13.1 Tensile Strength (MPa) 5.4 UCS (MPa) 105 φ (°) 63.9 Dilation angle (°) 5 Density (kg/m³) 2800 Fractures Cohesion (MPa) 0.25 Tensile Strength (MPa) 0 φ (°) 45 Dilation angle (°) 0 Table 5-1 Intact rock and fracture properties for the eight 20x20m models  5.1 8x4m Pillars When looking at the results from the 8x4 m pillar simulations with FLAC, one has to keep in mind that the elements are coarse compared to the size of the model.  As a consequence, the models are not as accurate and may exhibit unexpected behavior (i.e. a higher peak strength for the same model without confinement and with a confinement pressure of 1 MPa, see Figure 5-4).  The size of the elements is inherent to the method chosen for the upscaling approach whose ultimate goal is to be applied to much larger scale models with tens to hundreds of thousands of elements.  Here the models are purposely kept small and simple to allow rapid testing and better understanding of the different model responses.  Figure 5-2 and Figure 5-3 present the 10-80 and 30-60 FLAC models, respectively.  Due to the fact that the Most Critically Oriented Joint with Length Component upscaling approach is model based instead of cell based, it can be seen that the joint pattern remains consistent when 65  increasing the number of elements.  This characteristic helps reduce the undesirable mesh size dependency previously described.  Figure 5-4 compares the outcome of different FLAC simulations for the 10-80 sample with the ELFEN results.  Results from the FLAC and ELFEN simulations are generally in the same range and follow the same trend (good correlation with ELFEN results).  In addition, results from the FLAC simulations with 32, 50 and 72 cells yield similar values confirming that the upscaling approach produces consistent stress-strain responses for different mesh sizes.  The higher peak strength values at higher confining pressure for the 72-cell simulations are caused by the presence of an intact rock cell (i.e. cells without a fracture from the DFN mapped as crossing into it).  The difference increases with confinement because the mode of fracturing of the rock slowly shifts from being controlled by slip on ubiquitous joints to being controlled by yielding of intact rock elements.  As previously demonstrated by Fang and Harrison (2002) and Tang (1998, 2000a, 2000b), and as earlier demonstrated with the Mohr-Coulomb models in Chapter 4, failure through intact rock is better reproduced in FLAC when adopting very small mesh sizes.  In this case, the intact rock element is coarse compared to the model size due to the nature of the upscaling approach being tested and thus generates an undesirable strengthening effect.  This effect is only seen for the 72-cell simulation relative to the fracture intensity in the DFN used; the smaller the cell size, the more likely it is that a fracture won’t intersect a given cell and that intact rock properties would be assigned to it.   It can also be seen from Figure 5-4 that the 18 cells curve has a slightly different shape and constantly under predicts the strength of the sample.  This behavior is attributable to the coarseness of its elements (only 3 cells horizontally and 6 vertically).  This result therefore infers that oversized elements will yield less reliable results.  Figure 5-5 shows the outcomes of the different 30-60 FLAC models.  Again, results from the FLAC and ELFEN simulations are generally in the same range and follow the same trend (good correlation with ELFEN results).  The difference between the peak strength values from FLAC and ELFEN are significantly greater when there is no confinement because the loss of strength in the ELFEN model is mostly attributable to fragmentation under these conditions.  In other 66  words, failed blocks in ELFEN fragment and separate from neighboring blocks leading to increased sample dilation and reduced confinement (at least locally).  This behavior cannot be reproduced in FLAC because joints are not allowed to open and separate and the continuum mesh cannot break apart.  Figure 5-6 presents outcomes from both the 10-80 and 30-60 models together.  Simulations from the two models clearly exhibit different behavior and yield at different peak strength values even though the only difference is the joint orientations.  This thus shows one of the advantages of the upscaling approach over the use of empirical rock mass rating systems to assign model properties.  It also demonstrates that the upscaling approach has the ability to capture at least some of the characteristics of the rock mass fabric.                        Figure 5-2 10-80 FLAC models with: a) 18 cells, b) 32 cells, c) 50 cells, and d) 72 cells                        Figure 5-3 30-60 FLAC models with: a) 18 cells, b) 32 cells, c) 50 cells, and d) 72 cells  a) b) c) d) a) b) c) d) 67   Figure 5-4 Peak strength at different confining pressures for FLAC and ELFEN simulations for the 10-80 4x8m models  Figure 5-5 Peak strength at different confining pressures for FLAC and ELFEN simulations for the 30-60 4x8m models 68   Figure 5-6 Comparison between the 10-80 (10) and 30-60 (30) FLAC simulations  5.2 20x20m Blocks To further test the upscaling approach, it was applied to seven 20x20m models under uniaxial and triaxial simulations. Figure 5-7, Figure 5-8 and Figure 5-9 present the 20x20m DFN models tested.  The three sets of DFNs tested are based on different realizations from the same DFN with varying fracture intensities. Two models have a fracture intensity corresponding to the 10th percentile, three to the 50th percentile and two to the 90th percentile of the given P32 fracture intensity distribution     Again, these models are small (with respect to element numbers) compared to the scale of the models that are targeted for future use by the upscaling approach.  Thus, a 0.5m change in mesh size for the results presented here represents a higher relative change in scale than would be expected for a typical model and therefore a greater effect on the outcome than would otherwise be expected.   69               Figure 5-7 10th percentile 20x20m models a) p10a b) p10b               Figure 5-8 50th percentile 20x20m models a) p50a b) p50b c) p50c   a) b) a) b) c) 70             Figure 5-9 90th percentile 20x20m models a) p90a b) p90b  Results from the 20x20m FLAC models have been compared to outcomes from ELFEN.  Figure 5-10 through 5-16 present the FLAC models and Figure 5-17 through 5-18 present the peak strength for different confining pressures for the FLAC and ELFEN models.  It can be seen from these figures that the joint pattern is consistent for the different mesh sizes.  The upscaling approach, due to its model based approach, ensures that the most important features are represented regardless of mesh size.  This characteristic reduces the mesh size-based undesirable effects to their minimum.    Due to their low fracture intensity and the size of the models, both p10 models have a significant number of cells assigned intact rock properties, especially the 225 and 400 cell models (see Figure 5-10 and Figure 5-11).  This characteristic influences the predicted peak strength values because, as previously shown, shear localization and failure is only well captured in FLAC for very small element sizes.   a) b) 71               Figure 5-10 20x20m p10a FLAC models: a) 100 cells, b) 225 cells, and c) 400 cells               Figure 5-11 20x20m p10b FLAC models: a) 100 cells, b) 225 cells, and c) 400 cells  a) b) c) a) b) c) 72               Figure 5-12 20x20m p50a FLAC models: a) 100 cells b) 225 cells c) 400 cells               Figure 5-13 20x20m p50b FLAC models: a) 100 cells, b) 225 cells, and c) 400 cells  a) b) c) a) b) c) 73                 Figure 5-14 20x20m p50c FLAC models: a) 100 cells, b) 225 cells, and c) 400 cells               Figure 5-15 20x20m p90a FLAC models: a) 100 cells, b) 225 cells, and c) 400 cells  a) b) c) a) b) c) 74               Figure 5-16 20x20m p90b FLAC models: a) 100 cells, b) 225 cells, and c) 400 cells All peak strength plots show a good correlation between the FLAC and ELFEN models.  Except for the P10a and P90a models, the 400 cell models (1m mesh) generally provide a very good correlation.  The higher values for the 400 cell p10a simulation can be explained by the amount and position of the intact rock elements in the model.  As can be seen from Figure 5-10 the intact rock cells leave a sizeable gap between the upper left section and the rest of the model thus forcing the failure path through more intact rock cells.  As a result, it artificially increases the strength of the sample because it requires a lot of energy to break the 1m² elements. The plots for the P10b, P50b, P50c, and P90a DFN distributions show a good correlation between ELFEN and FLAC for higher confinement values most likely due to the failure kinematics captured in the models at lower confinement.  Because of the limitations of the ubiquitous joint model and continuum models in general, it is impossible to accurately reproduce kinematically driven events with the upscaling approach.  To illustrate this point, Figure 5-24 shows the displacements for two 20x20m FLAC models under gravity loading.  The first one is represented by a single large element with a 60° plane of weakness.  The second is comprised of 100 elements, for which 15 cells are assigned a 60° ubiquitous joint to simulate the 60° fracture a) b) c) 75  extending through the model.  Both models have a density of 2700 kg/m³ and have the same joint properties: tension=0 MPa, cohesion=0 MPa, friction angle=10°.  If the models capture the kinematics under zero confinement correctly, the upper left part of the model should slip under gravity because the friction angle is much lower than the joint angle.  However, the only movement is downward caused by gravity-driven consolidation of the model.  In order to properly reproduce the kinematics of this problem in FLAC, a 60° interface would have to be explicitly inserted in the model.  This technique is not a solution for representing fractured rock masses because the numerous interfaces required to appropriately represent the DFN would negate the computational efficiencies afforded in using a continuum-based code (i.e. FLAC). Each peak strength plot provided below shows at least two FLAC models with similar values.  This implies the existence of a mesh size effect and that the upscaling approach is only valid for a range of mesh sizes.  However, considering the size of the models relative to the element sizes, having two similar models while changing the mesh size by 0.5m suggests that the range for which the upscaling approach is valid is probably reasonably large for bigger models.  In most cases, the occurrence of different peak strength values is seen for the 100 cell models and is due to the large size of the elements.  In fact, for models with higher joint intensities, too many fractures are omitted at this scale and only the most critical fractures are accounted for.  In the 100 cell models, every element has significantly more influence on the overall strength of the model.  The critical fractures are thus given too much importance by controlling the behavior of large elements.  On the other hand, as previously discussed, Figure 5-10 and Figure 5-17 show that if the mesh size is too small for a model with low fracture intensity, more intact rock elements will be created and the model will be artificially strengthened.  Thus, a suitable mesh size for the upscaling approach must avoid overly large elements that assign too much importance to the fractures present while still capturing their localized influence, and not too small elements to avoid overpopulating the model with intact rock cells that artificially strengthen the model.  Consequently, the ideal mesh size is not a fixed number but depends on the fracture intensity and fracture size of the DFN.  Large models with persistent fractures and higher fracture intensity will perform better under a range of larger mesh sizes than small models with short fractures and low fracture intensity. 76   Figure 5-17 Peak strength at different confining pressures for FLAC and ELFEN simulations for the p10a 20x20m model   Figure 5-18 Peak strength at different confining pressures for FLAC and ELFEN simulations for the p10b 20x20m model 77   Figure 5-19 Peak strength at different confining pressures for FLAC and ELFEN simulations for the p50a 20x20m model   Figure 5-20 Peak strength at different confining pressures for FLAC and ELFEN simulations for the p50b 20x20m model 78   Figure 5-21 Peak strength at different confining pressures for FLAC and ELFEN simulations for the p50c 20x20m model   Figure 5-22 Peak strength at different confining pressures for FLAC and ELFEN simulations for the p90a 20x20m model 79   Figure 5-23  Peak strength at different confining pressures for FLAC and ELFEN simulations for the p90b 20x20m model                    Figure 5-24 Displacement for two 20x20m FLAC ubiquitous joint models with a 60 degree fracture   80  5.3 Reproduction of the Scale Effect and Failure Path Lorig (2007) points to the inability of rock mass rating systems to reproduce scale effects as being one of their biggest weaknesses.  So, for any upscaling approach to be successful, it must show weakening of the rock mass for an increase in volume size.  In order to further verify the ability of the Most Critically Oriented Joint with Length Component upscaling approach to reproduce the strength scale effect, five models of increasing size were sampled from the same DFN and tested under uniaxial compressive loading until failure.  Figure 5-25 shows the boundaries of the five models, ranging from 20 to 100 m in increments of 20 m, relative to the DFN model.  Each model was represented in FLAC using the Most Critically Oriented Joint with Length Component upscaling approach developed here and a 1x1 m mesh.  No models smaller than 20x20 m have been tested because of the limitations inherent to the upscaling approach for small models, relative to the fracture intensity, whose behavior is dominated by the presence of a small number of joints.  Table 5-2 presents the intact rock and fracture properties for the two simulations.  Figure 5-26 and Figure 5-27 respectively present the normalized UCS and normalized Young’s modulus for the five FLAC models.  Both graphs clearly exhibit a reduction of the rock mass strength properties with the increase in the DFN sampled area thus showing that the proposed upscaling approach can effectively account for scale effects.   81   Figure 5-25 Models used for testing the scale effect (dimensions in meters)  Intact Rock Bulk Modulus (MPa) 10000 Shear Modulus (MPa) 6000 Cohesion (MPa) 4.4 Tensile Strength (MPa) 1.5 φ (°) 57 Dilation angle (°) 5 Density (kg/m³) 2700 Fractures Cohesion (MPa) 0.08 Tensile Strength (MPa) 0 φ (°) 35 Dilation angle (°) 0 Table 5-2 Intact rock and fracture properties for the 60x60m model 82   Figure 5-26 Area vs. normalized UCS  Figure 5-27 Area vs. normalized Young’s modulus 83  In order to better understand how the upscaling approach captures the scale effect, a review of the plasticity indicators for the 60x60m and 20x20m model is presented (Figure 5-28). This figure shows the evolution of the shear zones in these models.  The 60x60m model was loaded under uniaxial compression until failure and three snapshots of the plasticity indicators were taken at different stages of the failure process.  Figure 5-29 shows the first stage of loading for the FLAC model.  The blue inverted Vs represent slip on the weak planes.  After a few hundred timesteps there is mostly only movement on the most critically oriented joints.  Figure 5-30 and Figure 5-31 present the second stage of shear localization and failure.  Tension has developed between fractures in the rock bridges and intact rock elements are beginning to fail in tension (pink circles).  Finally, Figure 5-32 shows the last stage of the failure path development, close to peak strength.  More intact rock cells have yielded and critically oriented fractures are now linked together and form continuous shear zones extending all the way through the model.  This behavior is consistent with observations made from loading intact rock samples and previous simulations by Fang and Harrison (2002), Clark (2006) and Tang (2000a, 2000b) where continuous shear zones developed through yielding of elements in tension prior to peak strength.  However, in this case the fracturing mechanism is controlled by the presence of the ubiquitous joints and the shear bands develop parallel to the orientation of the most critically oriented fractures instead of at the theoretical angle of 30° to the loading direction, as for intact rock samples. This shows the influence of the planes of weakness on the deformation and strength of the models. Furthermore, even if the less critically oriented joints are well represented in this model, they play little to no role in the behavior and strength of the model.  This therefore reinforces the hypothesis that, in most cases, picking only the most critically oriented joint in each cell is sufficient to capture the rock mass strength and deformation properties.  Figure 5-33 presents the plasticity indicators after the peak strength was reached for the 20x20m models.  Some slip on ubiquitous joints has occurred and a significant proportion of intact rock elements have failed in tension to form a localized shear zone.  When comparing this with the 60x60m model, the proportion and concentration of intact rock elements that have failed in tension is much smaller than the 20x20m model.  In fact, due to the smaller number of joints in 84  the 20x20m model the failure path must pass through (proportionally) more intact rock elements as well as bigger clusters of intact rock elements.  In fact, the ratio of intact rock fracture to ubiquitous joint slip for the 20x20m model is 1 compare to 2.3 for the 60x60m model.  This explains the increase in the sample strength properties with the increase in model area.  Sliding on ubiquitous joints requires less energy than breaking through intact rock elements.  As the 60x60m model contains more joints, a path going through a minimum number of intact rock elements can be more easily achieved.  However, in the case of the 20x20m model, the smaller number of fractures offers only a limited number of possible paths for the shear zone to localize along (involving weaker ubiquitous joints), resulting in paths that invariably must pass through bigger clusters of intact rock elements than in the 60x60m model.                     Figure 5-28 a) 60x60m model b) 20x20m model    a) b)  Figure  5-29 First stage of loading for the 60x60m model  85  Figure  5-30 Yielding in tension of the first rock bridges  86   Figure 5-31   Close up of the yielding in tension of the first rock bridges 87   Figure 5-32 Fully developed shear zones going all the way through the model after peak strength was reached  88  89    Figure 5-33 20x20m model after peak strength was reached   90  6 Application of Upscaling Approach to Pit Slope Models The previous chapters have described applications of the proposed approach to relatively small scale problems, discussing the critical effect of grid size on the modelling results.  In this chapter the attention is focused in applying the upscaling approach to large scale problems, such as engineered slopes.  The mechanical response of a 100m high slope with four different joint configurations has been investigated using the proposed upscaling approach and comparing the results with those obtained using both a continuum analysis with joint elements and a hybrid finite/discrete method. 6.1 DFN Models Figure 6-1 presents the external boundaries of the DFN models and principal dimensions of the slope models while Figure 6-2 to Figure 6-5 show the four joint patterns used in the simulations.  Each model is comprised of two sets of fractures, one sub vertical and one sub horizontal with different strength properties (see Table 6-1).  The sub vertical fractures have an average dip angle of -62°, -80°, 80° and 62° while the sub horizontal fractures have an average dip angle of 24°, -41°, 41° and -24° for Slope 1, 2, 3 and 4, respectively.  All joint sets have an average length of about 8m.   All four joint configurations were built based on the initial joint of (-62°/24°) by rotating one joint set at a time.   91   Figure 6-1 Slope geometry and DFN boundaries (dimensions in meters)    Figure 6-2 Slope1 DFN   92   Figure 6-3 Slope2 DFN    Figure 6-4 Slope3 DFN   93   Figure 6-5 Slope4 DFN  6.2 FLAC Ubiquitous Joint Rock Mass Models 6.2.1 Geometry and Properties To be consistent, all models simulated with all three programs used the same boundaries and excavation sequence.  Figure 6-6 shows the model limits and the 9 excavation stages (in blue).  Table 6-1 presents the intact rock and fracture properties for the slope models, with different strength properties for the sub vertical and sub horizontal fractures.  For the simulations, the horizontal stress was set at twice the vertical gravitational stress and the boundary conditions were as follows: the top is a free surface, the sides are restrained in the horizontal direction, and the bottom in the vertical direction.   94   Figure 6-6 Models geometry (all dimensions are in meters)  Intact Rock Bulk Modulus (MPa) 10000 Shear Modulus (MPa) 6000 Cohesion (MPa) 4.4 Tensile Strength (MPa) 1.5 φ (°) 57 Dilation angle (°) 5 Density (kg/m³) 2700 Sub Vertical Fractures Cohesion (MPa) 0.081 Tensile Strength (MPa) 0 φ (°) 51 Dilation angle (°) 0 Normal Stiffness (MPa/m) 40000 Shear Stiffness (MPa/m) 4000 Sub Horizontal Fractures Cohesion (MPa) 0.32 Tensile Strength (MPa) 0 φ (°) 35 Dilation angle (°) 0 Normal Stiffness (MPa/m) 40 Shear Stiffness (MPa/m) 4 Table 6-1 Intact rock and fracture properties for the slope models  95  6.2.2 Most Critically Oriented Joint with Length Component Upscaling Approach 6.2.2.1 Principal Stress Orientation As discussed in Chapter 5, the Most Critically Oriented Joint with Length Component upscaling approach depends on the orientation of the principal stresses.  It is safe to assume that the principal stress orientation will change as the slope is being excavated; therefore a methodology to select the shifting most critically oriented fracture had to be developed.  In order to capture the change in stress direction in the simulations as the slope is being excavated an elastic model was first run and the principal stress orientations recorded after each stage.  Once the principal stress orientation was resolved for each cell of the model (at each stage), the information is then used to determine the most critical fracture at every stage for every cell.  Subsequently, during the simulation, the fracture orientation for the UJRM is updated at the end of each excavation stage based on results obtained from the elastic model.  This step is essential to capture the changing influence of joints on the behaviour of the model as the stress field changes. 6.2.2.2 Mesh Size Effect In order to investigate the mesh size effect for large scale models, the excavation of Slope 1 and 2 were simulated with a 1 and 2 m mesh size.  In addition, Slope1 was run with a 0.5 m mesh.  Figure 6-7, Figure 6-8 and Figure 6-9 present a close-up of the mesh for the three Slope1 models.  As can be seen from the figures, the 0.5 m mesh is more or less a raster representation of the DFN while the 1 and 2 m mesh exhibit a smaller intact rock cell to ubiquitous joint cell ratio.  Figure 6-10, Figure 6-11 and Figure 6-12 present the horizontal displacements, horizontal stresses and vertical stresses, respectively, predicted by the FLAC simulations for Slope1 with a 0.5, 1 and 2 m mesh.  Results for Slope 2 are presented in Appendix D.  The figures clearly show that only minor differences exist between the different mesh sizes in terms of predicted stresses and displacements.  The main noticeable difference between the different simulations is the amount of horizontal tension that develops (shown in black in Figure 6-11).  As can be seen in Figure 6-10, the amount of tensile failure in the model is directly proportional to the ratio of intact rock cells to ubiquitous joint cells.  This relationship has little influence here because most of the deformation is accommodated by sliding on ubiquitous joints.  However, if a large amount 96  of tensile failure is expected in a simulation, then it would probably be better accommodated with a coarser mesh rather than a raster like model as in the case for the 4x8m and 20x20m models.  Another noticeable characteristic of the Slope1 and Slope2 models are slightly higher displacements for coarser mesh sizes.  This tendency has previously been observed with the 4x8m and 20x20m models (lower peak strength for a coarser mesh). However, the effect is marginal for the large scale models and is not worth the extra computing time necessary to run models with a finer mesh (going from a 1 m mesh to a 2 m mesh reduces the run time by a factor greater than 10).   Figure 6-7 Windowed view of the Slope 1, 0.5 m mesh  15m 97   Figure 6-8 Windowed view of the Slope 1, 1 m mesh   Figure 6-9 Windowed view of the Slope 1, 2 m mesh     15m 15m 98        Figure 6-10 Slope1  horizontal displacement contours, in meters, for: a) 0.5m mesh, b) 1m mesh, and c) 2m mesh, showing only minor differences a) b) c) 40m 99         Figure 6-11 Slope 1 horizontal stress contours, in Pa, for: a) 0.5m mesh, b) 1m mesh, and c) 2m mesh,  showing only minor differences  a) b) c) 40m 100        Figure 6-12 Slope1 vertical-stress contours, in Pa, for: a) 0.5m mesh, b) 1m mesh, and c) 2m mesh, showing only minor differences   a) b) c) 40m 101  6.2.2.3 Influence of Fracture Orientation Four FLAC models with the same intact rock and fracture properties were simulated using the upscaling approach to investigate the influence of fracture orientation on the modelled slope behaviour.  It is expected that, if the upscaling approach is properly capturing the rock mass fabric it should yield different results as a function of the orientation of the predefined joint sets.   In order to qualitatively compare the behaviour of the four slope models, factor of safety (FOS) simulations have been carried out.  Although, the FOS for the four slopes is very similar, 11.39, 11.30, 11.44 and 11.56 respectively, Figure 6-13 through Figure 6-16 clearly exhibit differences in the size and shape of the failure plane for the different slope models.  Slope4, due to the favourable orientation of its joint sets presents a thin well defined failure plane whereas the failure plane of Slope1 is more massive and failed elements reach deeper into the rock mass.  The only difference between the two slopes, and consequently the only explanation for this discrepancy, is the orientation of one of the two joint sets.  Similar observations can also be made when comparing Slope2 and Slope3 together and again, the only dissimilarity is the orientation of one joint set.  This demonstration clearly shows that both joint sets represented in the upscaling approach play a role in the strength and behaviour of the rock mass simulated.  Similar observations can also be made when analyzing displacements.  Figure 6-17 shows the horizontal displacements after the final stage of excavation of Slope1, Slope2, Slope3 and Slope4.  It can be seen that each model produces a slightly different response with Slope2 exhibiting the least amount of displacement and Slope4 the most.  These results once again confirm that the upscaling approach captures the influence of the different joint set orientations on the behaviour of the models even if one set is clearly more dominant.  The fact that Slope4 predicts the greatest displacements while Slope2 predicts the least, agrees with what would be expected given the more critical joint orientations in Slope4 relative to Slope2.  Figure 6-17 also shows more movement controlled by the ubiquitous joints closer to surface as indicated by the stepped and irregular nature of the displacement contours.  This phenomenon is a direct consequence of the increase of confinement with depth.  Close to the surface, the behaviour of the material is mostly controlled by displacement on weak planes while at depth, where confinement increases the normal stress on the joints, the rock mass behaves more like a 102  homogenous material.  Figure 6-18 shows the plasticity indicators for the slope models.  It can be seen that, as previously mentioned, displacements on the ubiquitous joints occurs close to the surface, whereas no yielding or slip on the weak planes occurs at depth.   Figure 6-13 Factor of safety analysis plasticity indicators Slope1      40m 103   Figure 6-14 Factor of safety analysis plasticity indicators Slope2  Figure 6-15 Factor of safety analysis plasticity indicators Slope3    40m         40m 104   Figure 6-16 Factor of safety analysis plasticity indicators Slope4         40m 105          Figure 6-17 Influence of fracture orientation on horizontal displacements: a) Slope1, b) Slope2, c) Slope3, and d) Slope4. Displacements are reported in meters a) b) c) d) 40m   Figure  6.3 Comparison of FLAC model with ELFEN and PHASE²In this section simulations for Slope1 and Slope2 are compared, in terms of predicted horizontal displacement and horizontal stress, to similar models run in ELFEN and Phase² in an effocompare the outcomes of the FLAC large scale models 6.3.1 Horizontal DisplacementFigure 6-19, Figure 6-20 and the excavation of Slope1 carried out respectively with ELFEN, Phase² and FLAC.  It can be seen from the figures that the FLAC and Phase² models yield vdisplacement magnitudes and model behaviour. In contrast, the ELFEN model produces slope displacements that are two to three times higher (5contours follow more closely the orientationthose in Phase² and FLAC are more representative of a continuum. These discrepancies can be explained by the way the joints are represented and are allowed to move and develop in every 6-18 Plasticity indicators for a typical slope model  with different numerical techniques. Figure 6-21 present the horizontal displacements for simulations of ery similar results in terms of -10 cm instead of 2-4 cm). Displacement  of the sub horizontal joint set in ELFEN, whereas         106  rt to .    40m 107  model.  As previously noted in the methodology section, the joints in the ELFEN model are explicitly represented and are allowed to grow and slip along their entire length. Those in Phase² are also explicitly represented but are modeled as terminating with both ends closed, thus limiting the amount of possible displacement, and cannot propagate.  In the FLAC models, the joints are implicitly represented and displacement occurs on a cell by cell basis through the ubiquitous joint formulation.  Moreover, joints are not allowed to grow or open.  In this sense, the Phase² and FLAC treatments of the DFN are more similar in the way they handle fractures, thus yielding comparable results.    As the slope is conceptual, no monitoring data is available to verify one response to the others in terms of correctness.  Regardless, the characteristics aforementioned as to how the joints are represented and modeled by the three programs signify that both continuum approaches (FLAC and Phase²) cannot fully simulate the kinematic effects properly.  In fact, as movement on joints is limited and fractures are bounded to their original size in the continuum-based treatments, toppling, bi-planer failure and other failure modes involving a significant kinematic component in their development cannot be properly reproduced.   108   Figure 6-19 Slope1 horizontal displacements modeled in ELFEN. Displacements are reported in meters   Figure 6-20 Slope1 horizontal displacements modeled in Phase². Displacements are reported in meters       40m       40m 109    Figure 6-21 Slope1 horizontal displacements modeled in FLAC. Displacements are reported in meters  Figure 6-22, Figure 6-23 and Figure 6-24 present the horizontal displacements for simulations of the excavation of Slope2 carried out respectively with ELFEN, Phase² and FLAC.  All three program yield very similar displacement values and contours.  Consequently, the upscaling approach seems to be an effective tool to simulate the behaviour of a fractured rock mass when no kinematics is involved.         40m 110   Figure 6-22 Slope2 horizontal displacements modeled in  ELFEN. Displacements are reported in meters   Figure 6-23 Slope2 horizontal displacements modeled in Phase². Displacements are reported in meters       40m      40m 111   Figure 6-24 Slope2 horizontal displacements modeled in FLAC. Displacements are reported in meters  6.3.2 Horizontal Stresses and Tension Figure 6-25 through Figure 6-30 present the horizontal stresses for of the Slope1 and Slope2 simulations carried out with ELFEN, Phase² and FLAC.  All simulations yield similar horizontal stress values and contours.  When comparing the models, it can be observed that the only major difference between them is the amount of tensile stress that develops.  Both the ELFEN and FLAC Slope1 models predict tensile stress close to the excavation boundary from the top to almost the bottom of the slope, whereas Phase2 only shows zero or tensile stress along the top of the model.  The tension in the ELFEN and FLAC models is attributable to slip on joints close to the surface, leading to extensional deformations and tensile stresses at the crest.  The absence of this characteristic in the Phase2 model is probably due to the fact that both ends of the joints in the rock mass are closed and therefore cannot accommodate as much displacement and generate tensile stress.         40m 112  The fact that all three programs predict similar stress contours and magnitudes helps to validate the upscaling approach and shows that it represents a viable alternative to more computationally demanding methods when modeling a fractured rock mass if the expected mode of failure does not involve significant kinematic controls.   Figure 6-25 Slope1 horizontal stresses modeled in ELFEN. Stresses are reported in Pa        40m 113   Figure 6-26 Slope1 horizontal stresses modeled in Phase². Stresses are reported in MPa   Figure 6-27 Slope1 horizontal stresses modeled in FLAC. Stresses are reported in Pa     40m      40m 114   Figure 6-28 Slope2 horizontal stresses modeled in ELFEN. Stresses are reported in Pa   Figure 6-29 Slope2 horizontal stresses modeled in Phase². Stresses are reported in MPa   40m      40m 115   Figure 6-30 Slope2 horizontal stresses modeled in FLAC. Stresses are reported in Pa          40m 116  7 Conclusion and Recommendations 7.1 Conclusion A geomechanical upscaling approach has been developed in an attempt to better capture the rock mass structure in continuum simulations.  Many different ways of reproducing the behavior of a fractured rock mass have been tested.  Initially, Mohr-Coulomb Raster Models have been experimented based on previous successful attempts to model the brittle failure of intact rock by Fang and Harrison (2001, 2002) and Tang (1995, 1997).  This method has proven efficient to simulate the behavior of a fractured rock mass and had the advantage of taking into account all fractures comprised in the DFN.  However, the number of elements necessary to obtain accurate results was proven prohibitively high to be used in large scale simulations.  Still using the Mohr-Coulomb constitutive model, Reduced Strength “GSI” Models have been tested to simulate the fractured rock mass.  Although they could be implemented with a larger mesh size, no satisfactory correlation between GSI and joint properties has been found making it not feasible to find an automated way to generate the continuum models.  Furthermore, this method does not provide any directionality which is an important issue when trying to capture the influence of the rock mass fabric.  The choice has then been made to use the UJRM constitutive model in order to capture the joint directionality while being able to generate models with a coarser mesh size.  Multiple approaches have been investigated to choose the orientation of the plane of weakness.  Both the random and the most represented approaches have proven inaccurate as they both greatly overestimated the strength of the samples and did not seem to accurately capture the influence of the rock mass fabric on the stress distribution.    Consequently, a mechanistic approach to the choice of the ubiquitous joint orientation has been developed to remediate to the problems encountered with the previous method.  Good results were obtained when combining the relative strength of the joints, based on their orientation 117  relative to the principal stress orientation, and a joint length component.  This method yielded peak strength values close to expected and a redistribution of the stresses characteristic of a fractured rock mass.  The Most Critically Oriented Joint with a Length Component has therefore been chosen for the upscaling approach.  This newly developed method of representing a fractured rock mass in a continuum simulation has subsequently been tested on two 4x8m pillar models and seven 20x20m models.  Results have been found to closely match similar simulations carried out with ELFEN® in terms of trend and range of values of peak strength.  The FLAC® simulations seemed however to be slightly mesh dependent and yielded higher peak strength for low confinement values, due to the incapacity of continuum codes to explicitly capture kinematics.  The upscaling approach has then been tested on five samples of different size from the same DFN to verify its ability to reproduce the scale effect.  Through an evaluation of the plasticity indicators and strength properties values of the five models, it has been found that the upscaling approach has the ability to capture, at least to some extent, the scale effect.  Finally, the proposed upscaling approach has been evaluated by simulating the excavation of the same slope model with four different joint configurations.  The outcomes of the simulations clearly showed that the influence of both joint sets was captured in every simulation even with a 2m mesh size.  Furthermore, after examining the horizontal displacements and stresses, it has been found that the mesh size effect previously found in the smaller models had only minor consequences for large scale models.  Lastly, a comparison of Slope1 and Slope2 models with similar simulations run in ELFEN® and Phase²® revealed that the upscaling approach is a viable substitute to other methods of representing a fractured rock mass as long as kinematics is not involved in the primary expected mode of failure.  Through the different tests, the upscaling approach developed for this research has proven to be a fast and objective way of assessing the rock mass properties while staying as close as possible to the data acquired in the field and from laboratory testing.  It also captures, to some extent, the scale effect and is less computationally intensive and needs less preparation than the explicit 118  models.  However, the upscaling approach, due to limitations inherent to continuum simulations, is incapable of properly capturing kinematics.  Consequently, it offers a viable alternative to other methods to model the behavior of a fractured rock mass under stress as long as kinematics is not the primary mode of failure expected 7.2 Recommendations for Further Work Some of the recommendations for further work include verifying and improving the equation used to choose the orientation of the ubiquitous joint.  The present formulation yielded good results, however it is based on a theoretical equation for the joint strength that does not comply exactly with the reality and ratio between joint length and strength that should also be further investigated.  It is also suggested that the approach be extended to 3D simulations as it has been intended for since the beginning of this research.  The ultimate goal of the upscaling approach is being able to model 3D problems while realistically capturing the influence of rock mass structure in a computationally efficient way.  Because of the constraint of time imposed by the Master’s program all simulations have been conducted in 2D for this research.  It could be also interesting to see if a better constitutive model incorporating some characteristics of the upscaling approach and the UJRM model could be developed to increase the accuracy of the method and eliminate the need to manually change the orientation of the weak plane when the orientation of the principal stress changes.  Finally, it is also recommended to test the upscaling approach by back analyzing one or more slope failures and see how it is able to capture the influence of the rock mass structure.    119  References BAECHER, G.B., LANNEY, N.A. and EINSTEIN, H.H. (1977), "Statistical description of rock properties and sampling," Proceedings of the 18th U.S. Symposium on Rock Mechanics, American Institute of Mining Engineers, 5C1-8  BAECHER, G.B. (1983), "Statistical analysis of rock mass fracturing", Mathematical Geology, Vol. 15, No. 2, pp. 329-348  BAKUN-MAZOR, D., HATZOR, Y.H., DESHOWITZ, W.S. (2009), “Modelling mechanical layering effects on stability of underground openings in jointed sedimentary rocks”, International Journal of Rock Mechanics, Mining Sciences and Geomechanics Abstracts, Vol. 46, pp. 262-271  BARTON, C. M. (1977), Geotechnical analysis of rock structure and fabric in CSA Mine NSW, Applied Geomechanics Technical Paper 26 Commonwealth Scientific and Industrial Research Organisation, Australia  BARTON, C.M., (1978), “Analysis of joint traces”, in 19th U.S. Symposium on Rock Mechanics  BARTON, N., (1988), “Rock mass classification and tunnel reinforcement selection using the Q-system”, Rock Engineering Systems for Engineering Purposes, ASTM STP 984, Louis Kirkaldie, Ed., American Society for Testing Materials, Philadelphia, pp. 59-88  BARTON, N.R., LIEN, R. and LUNDE, J., (1974), “Engineering classification of rock masses for the design of tunnel support”, Rock Mech, 6(4), pp. 189-239  BIENIAWSKI, Z.T., (1976), “Engineering classification in rock engineering”, In:  Proceedings of the Symposium on Exploration for Rock Engineering, Johannesburg, pp. 97-106  BIENIAWSKI, Z.T., (1988), “The rock mass rating (RMR) system (Geomechanics classification) in engineering practice”, Rock Engineering Systems for Engineering Purposes, ASTM STP 984, Louis Kirkaldie, Ed., American Society for Testing Materials, Philadelphia, pp. 17-34  BIENIAWSKI, Z.T., (1989), “Engineering rock mass classifications”, John Wiley & Sons, New York, p. 251  BIENIAWSKI, Z.T., (1993), “Classification of Rock Masses for Engineering: The RMR System and Future Trends”, In: Comprehensive Rock Engineering: Principles, Practice & Projects, Volume 3 Rock testing and site characterization, pp. 553-572  BRADY, B.H.G., and BROWN, E.T., (2004), “Rock mechanics for underground mining – Third edition”, Springer, Netherlands, 626 p  BRIDGES, M.C. (1976), “Presentation of Fracture Data for Rock Mechanics,” 2nd Australian-New Zealand Conference on Geomechanics, Brisbane, pp. 144-148 120  BROOK, N. and DHARMARATNE, P.G.R., (1985), “Simplified rock mass rating system for mine tunnel support”, Transactions of the Institution of Mining & Metallurgy, Section A, 94, pp. 148-154  BROWN, E.T., (1981), “Putting the NATM into perspective”, Tunnels and Tunneling, pp. 13-17  CAI M. et al., (2004), “Estimation of rock mass deformation modulus and strength of jointed hard rock masses using the GSI system” International Journal of Rock Mechanics & Mining Sciences, Vol. 41, pp. 3-19  CAI, M., KAISER, P.K., (2004) “Numerical simulation of the Brazilian test and the tensile strength of anisotropic rocks and rocks with pre-existing cracks”, Proceedings of the SINOROCK international symposium on rock mechanics: rock characterization, modelling and engineering design methods, 18–21 May, Three Gorges Project site, China. Paper 2B-03  CALL, R.D., SAVELY, J. and NICHOLAS, D.E. (1976), “Estimation of joint set characteristics from surface mapping data”, 17th U.S. Symposium on Rock Mechanics, pp. 2B2-1 - 2B2-9  CLARK, I.H. (2006), “Simulation of Rockmass Strength Using Ubiquitous Joints”, In: R. Hart & P. Varona (eds), Numerical Modeling in Geomechanics — 2006; Proc. 4th International FLAC Symposium, Madrid, May 2006, Paper No. 08-07, Minneapolis: Itasca  COATES, D. F. and PARSONS, R.C., (1966), “Experimental criteria for classification of rock substances”, Int. J. Rock Mech. Min., 3, pp. 181-189  COATES, D. F. and PATCHING, T.H., (1968), “A recommended rock classification for rock mechanics purposes”, CIM Bulletin, October, pp. 1195-1197  COATES, D.F., (1964), “Classification of rocks for rock mechanics”, Int. J. Rock Mech. Min., 1, pp. 421-429  CUMMINGS, R.A., KENDORSKI, F.S. and BIENIAWSKI, Z.T., (1982), “Caving rock mass classification and support estimation”, U.S. Bureau of Mines, Contract Report  #J0100103 [Referred in Hoek, 2000]  DEERE, D. U. and MILLER, R.P. (1966), “Engineering classification and index properties for intact rock”, Tech. Rept. No., AFWL-TR-65-116, Air Force Weapons Lab., Kirtland Air Force Base, New Mexico p. 300   DEERE, D. U., (1968), “Geological considerations”, In: Rock mechanics in engineering practice, (R. G. Stagg and D. C. Zienkiewicz.(Eds.)), Division of Civil Engineering, School of Engineering, University of Wales, Swansea, John Wiley & sons, New York, 1-20  DEERE, D.U., (1989), “Rock quality designation (RQD) after 20 years”, U.S. Army Corps Engrs, Contract Report GL-89-1, Vicksburg, MS: Waterways Experimental Station  121  DEERE, D.U., DEERE, D.W., (1988), “The rock quality designation index in practice”, Rock Engineering Systems for Engineering Purposes, ASTM STP 984, Louis Kirkaldie, Ed., American Society for Testing Materials, Philadelphia, pp. 17-34  DERSHOWITZ, W.S. (1984), “Rock joint systems”, PhD dissertation, M.I.T., Cambridge, Massachusetts.  DERSHOWITZ, W.S., (1979) “Jointed Rock Mass Deformability: A Probabilistic Approach” S.M. Thesis, Massachusetts Institute of Technology, Cambridge, MA  DERSHOWITZ, W.S., EINSTEIN, H.H., (1988), “Characterizing rock joint geometry with joint system models”, Rock Mechanics and Rock Engineering, Vol. 21, pp. 21-51  DERSHOWITZ, W.S., EINSTEIN, H.H., (1988), “Characterizing rock joint geometry with joint system models”, Rock Mechanics and Rock Engineering, Vol. 21, pp. 21-51  DERSHOWITZ, W.S., LA POINTE, P.R., CLADOUHOS, T. (1998), “Derivation of Fractures Spatial Pattern Parameters from Borehole Data”, International Journal of Rock Mechanics, Mining Sciences and Geomechanics Abstracts, Vol. 35, p. 508  DIEDERICHS, M.S., (2003), “Rock fracture and collapse under low confinement conditions”, Rock Mech. Rock Eng., Vol. 36, pp. 339-381  EBERHARDT, E (2009), “Weak rock Tunneling, PowerPoint Presentation”, EOSC-547 Lecture Material  EBERHARDT, E., STEAD, D., COGGAN, J.S., (2004) “Numerical analysis of initiation and progressive failure in natural rock slopes-the 1991 Randa rockslide” Int. J. Rock Mech. Min. Sci., Vol. 41, pp. 69–87  EDELBRO, C., (2003), Rock Mass Strength A Review, Technical Report, Lulea University of Technology, URL: http://epubl.luth.se/1402-1536/2003/16/LTU-TR-0316-SE.pdf, May 2010  EINSTEIN, H.H. and al. (1980), “Risk analysis for rock slopes in open pit mines -- final technical report”, Publication No. R80-17, Order No. 669, Department of Civil Engineering, M.I.T., Cambridge, Massachusetts  ELMO, D., (2006) “Evaluation of a hybrid FEM/DEM approach for determination of rock mass strength using a combination of discontinuity mapping and fracture mechanics modelling, with particular emphasis on modelling of jointed pillars” Ph.D. thesis, Camborne School of Mines, University of Exeter, UK  ELMO, D., and STEAD, D. (2010) “An integrated numerical modelling-discrete fracture network approach applied to the characterisation of rock mass strength of naturally fractured pillars”, Rock Mech. Rock Eng. 122  ELMO, D., COGGAN, J.S., PINE, R.J., (2005) “Characterisation of rock mass strength by combination of field mapping and numerical modelling” Proceedings of the 40th US rock mechanics symposium, Anchorage, Alaska  ELMO, D., VYAZMENSKY, A., STEAD, D., RANCE, J. (2008) “Numerical analysis of pit wall deformation induced by block caving mining” Proceedings of the 5th Conference and Exhibition on Mass Mining, Lulea, Sweden, June 2008, pp. 1073-1082  FANG, Z., HARRISON, J.P., (2001), “A mechanical degradation index for rock”, International Journal of Rock Mechanics & Mining Sciences, Vol. 38, pp. 1193-1199  FANG, Z., HARRISON, J.P., (2002a), “Numerical analysis of progressive fracture and associated behavior of mine pillars by use of a local degradation model”, Mining Technology, Volume 111, Number 1, April 2002 , pp. 59-72(14)  FANG, Z., HARRISON, J.P., (2002b), “Development of a local degradation approach to the modeling of brittle fracture in heterogeneous rocks”, International Journal of Rock Mechanics & Mining Sciences, Vol. 39, pp. 443-457  GOEL, R.K., JETHWA, J.L. and PAITHANKAR, A.G., (1996), “Correlation between Barton´s Q and Bieniawski's RMR – A new approach”, Technical Note, Int. J. Rock Mech. Min., 33 (2), pp. 179-181  Golder Associates, (2008), “An Introduction to Discrete Fracture Network Modelling using FracMan”, Course Manual  GOODMAN, R.E., SHI, G., (1985.) “Block theory and its application to rock engineering”, Englewood Cliffs, NJ: Prentice-Hall, p. 338  HOEK, E., CARRANZA-TORRES, C. & CORKUM B. (2002) “Hoek- Brown failure criterion” Proc. 5th NARMS, Toronto, July 2002, pp 267-273  HOEK, E. (2000). Practical rock engineering. URL: http://www.rocscience.com/hoek/pdf/3_Rock_mass_classification.pdf, May 2010"  HOEK, E., (2001), “Rock mass properties for underground mines”, In: Underground Mining  Methods: Engineering Fundamentals and International Case Studies, (Edited by W. A. Hustrulid and R. L. Bullock), Littleton, Society for Mining, Metallurgy, and Exploration, Inc. (SME)  HOEK, E., KAISER, P.K. and BAWDEN, W.F., (1995), “Support of underground excavations in hard rock” A.A. Balkema/Rotterdam/Brookfield  http://rfpa.cn/en/, MAY 2010  http://www.itascacg.com/pfc2d/, MAY 2010  123  http://www.rockfield.co.uk/elfen.htm, MAY 2010  HUDSON, J.A. and LA POINTE, P.R. (1980), “Printed Circuits for Studying Rock Mass Permeability,” Technical Note, International Journal of Rock Mechanics, Mining Sciences and Geomechanics Abstracts, Vol. 17, pp. 297-301  HUDSON, J.A. and PRIEST, S.D. (1979), “Discontinuities and rock mass geometry”, International Journal of Rock Mechanics, Mining Sciences and Geomechanics Abstracts, Vol. 16, pp. 339-362  International Society for Rock Mechanics (ISRM), (1981), “Commission on classification of rocks and rock masses, International Journal of Rock Mechanics and Mining Abstract, Vol. 18, pp. 85-110  Itasca Consulting Group, Inc. (2008a), FLAC – Fast Lagrangian Analysis of Continua, Ver. 6.0, User’s Manual  Itasca Consulting Group, Inc. (2008b), FLAC – Fast Lagrangian Analysis of Continua, Ver. 6.0, Verification Problems  Itasca Consulting Group, Inc. (2008c) FLAC – Fast Lagrangian Analysis of Continua, Ver. 6.0, Theory and Background  JING, L., and HUDSON J.A. (2002) “Numerical methods in rock mechanics” International Journal of Rock Mechanics, Mining Sciences and Geomechanics Abstracts, Vol. 39, pp. 409-427  JING, L., STEPHANSSON, O., (1994) “Topological identification of block assemblages for jointed rock masses” International Journal of Rock Mechanics, Mining Sciences and Geomechanics Abstracts, NO 2, pp. 163-172  KENDORSKI, F.S., CUMMINGS, R.A., BIENIAWSKI , Z.T. and SKINNER, E., (1983), “Rock mass classification for block caving mine drift support”, Proc. 5th congr. Int. Soc. Rock Mech., Melbourne, B51-B63 [Referred in Hoek, 2000]  LA POINTE, P.R., and HUDSON, J.A. 1985. “Characterization and Interpretation of Rock Joint Patterns”, GSA Special Paper 199, GSA, Denver  LAUBSCHER, D.H., (1975), “Geomechanics classification of jointed rock masses – Mining applications”, T. I. Min. Metal, A, 86, pp. A1-A8  LAUBSCHER, D.H., (1990), “A Geomechanics classification system for the rating of rock mass in mine design”, J. S. Afr. Inst. Min. Metall., Vol. 90, No. 10, pp. 257-273  LAUBSCHER, D.H., and TAYLOR, H. W., (1976), “The importance of Geomechanics classification of jointed rock masses in mining operations”, In: Proceedings of the Symposium 124  on Exploration for Rock Engineering, Johannesburg, A. A. Balkema, 1, pp. 119- 128 [Referred in Edelbro, 2003]  LAUFFER, H., (1958), “Classification for tunnel construction”, Geologie und Bauwesen, 4(1), 46-51, [In German] [Referred in Edelbro, 2003]  LORIG , L.J., (2007), “Using numbers from geology”, 11th congress of the international Society for Rock mechanics, pp. 1369-1377  MAS IVARS, D., PIERCE, M., DEGAGNÉ, D. & DARCEL, C. (2008), “Anisotropy and scale dependency in jointed rock mass strength – A synthetic rock mass study”, In R. Hart, C. Detournay & P. Cundall (eds), Proceedings of the First International FLAC/DEM Symposium on Numerical Modeling, August 25-27, 2008, Minneapolis, MN, USA, Minneapolis: Itasca  McMAHON, B.K. (1974), “Design of rock slopes against sliding on pre-existing fractures,” Proceedings, International Congress on Rock Mechanics, pp. 803-808  MEYER, T., EINSTEIN, E.E. (2002), “Geologic stochastic modeling and connectivity assessment of fracture systems in the Boston Area” Rock Mechanics and Rock Engineering, Vol. 35, pp. 23–44  MIN, K.-B., JING, L. (2004), “Stress Dependent Mechanical Properties and Bounds of Poisson’s Ratio For Fractured Rock Masses Investigated by DFN-DEM Technique”, International Journal of Rock Mechanics, Mining Sciences and Geomechanics Abstracts, Vol. 41, No. 3, CD-ROM  PALMSTRÖM, A., (1982), “The volumetric joint count - a useful and simple measure of the degree of rock jointing” Proc. 4th congress Int. Assn Engng Geol., Delhi 5, pp. 221-228  PALMSTRÖM, A., (1996), “RMi – A system for characterising rock mass strength for use in rock engineering” J. of Rock Mech. and Tunnelling Tech., India, 1 (2), pp. 69-108 [Referred in Edelbro, 2003]  PALMSTRÖM, A., and BROCH, E., (2006), “Ground behaviour and rock engineering design tools”, Tunnelling Underground Space Technology, pp. 575-593  PALMSTRÖM, A.,and  STILLE, H., (2006), “Ground behaviour and rock engineering design tools”, Tunnelling Underground Space Technology  PARK, E-S, MARTIN, C. D., CHRISTIANSSON, (2004), “Numerical simulation of the mechanical behaviour of discontinuous rock masses” In Y. Shimizu et al. (eds), Numerical modeling in micromechanics via particle methods-2004 (Kyptp, Japan, October 2004) pp. 85-91  PIERCE, M., CUNDALL, P., POTYONDY, D. and MAS IVARS, D., (2007) “A Synthetic Rock Mass Model for Jointed Rock”, Rock Mechanics Meeting Society’s Challenges and Demands – E. Eberhardt, D. Stead and T. Morrison (eds), Taylor and Francis Group, London, pp. 341–349  125  PINE, R.J. and al. (2007), “A new discrete modelling approach for rock masses”, Géotechnique, Vol. 57, No. 9, pp. 757-766  PINE, R.J., COGGAN, J.S., FLYNN, Z.N., ELMO, D., (2006), “The Development of a new Numerical Modelling Approach for Naturally Fractured Rock Masses”, Rock Mech. Rock Eng. Vol. 39, pp.395-419  PITEAU, D.R. (1970), “Geologicla Factors Significant to the Stability of Slopes Cut in Rock,” Proceedings, Symposium on the Theoretical Background to the Planning of Open Pit Mines with Special Reference to Slope Stability, pp. 33-53  PRIEST, S.D. and HUDSON, J. (1976), “Discontinuity spacings in roc”, International Journal of Rock Mechanics, Mineral Science and Geomechanics Abstracts, Vol. 13, pp. 135-148  RABCEWICZ, L.V., (1964/65), “The new Austrian tunnelling method”, Water Power, art 1 November 1964, 511-515, Part 2 January 1965, 19-24  RAMAMURTHY, T. and ARORA, V.K., (1993), “A classification for intact and jointed rocks”, In: Geotechnical engineering of hard soils - Soft Rocks, Anagnostopoulos et al., (Eds.), Balkema, Rotterdam [Referred in Edelbro, 2003]  RIEDMÜLLER, G. and SCHUBERT, W, (1999), “Critical comments on quantitative rock mass classifications”, Felsbau, 17(3)  ROBERTSON, A. (1970), “The Interpretation of Geological Factors for Use in Slope Theory,” Proceedings, Symposium on the Theoretical Background to the Planning of Open Pit Mines with Special Reference to Slope Stability, pp. 55-71  ROBERTSON, A., (1988), “Estimating weak rock strength”, AIME annual general meeting, Tucson, Arizona  ROCSCIENCE, (2010), Help Manual Phase² Version 7.0, Electronic Format  ROMANA, M., (1985), “New adjustment ratings for application of Bieniawski classification to slopes”, In: Proceedings of the International Symposium on the Role of Rock Mechanics, pp. 49- 53 [Referred in Edelbro, 2003]  SAINSBURY, B., PIERCE, M. and MAS IVARS, D. (2009), “Simulation of rock mass strength anisotropy and scale effects using a Ubiquitous Joint Rock Mass (UJRM) model”, Proceedings First International FLAC/DEM Symposium on Numerical Modelling, 25–27 August, Minneapolis, USA, Minneapolis, Itasca  SINGH, B. and GOEL, R.K., (1999), “Rock mass classification-A practical approach in civil engineering”, Elsevier, Netherlands  126  SJÖBERG, J. (1999), “ Analysis of large scale rock slopes” Doctoral thesis Luleå University of Technology, Department of Civil and Mining Engineering, Division of Rock Mechanics  SKINNER, E.H., (1988), “A ground support prediction concept: the rock structure rating (RSR) model”, Rock Engineering Systems for Engineering Purposes, ASTM STP 984, Louis Kirkaldie, Ed., American Society for Testing Materials, Philadelphia, pp. 33-51  SNOW, D.T. (1965), “A Parallel Plate Model of Fractured Permeable Media”, Ph.D. Dissertation, University of California, Berkeley, pp. 331  SNOW, D.T. (1968), “Rock fracture spacings, openings, and porosities”, ASCE Journal, Vol. 94, pp. 73-91  SNOW, D.T. (1970), “The frequency and apertures of fractures in rock”, International Journal of Rock Mechanics and Mineral Science, Vol. 7, pp. 23-40  SONG, J., LEE, C., SETO, M. (2001), “Stability analysis of rock blocks around a tunnel using a statistical joint modeling technique “, Tunnelling and Underground Space Technology, Vol. 16, pp. 341-351  STAUB, I., FREDRIKSSON, A., OUTERS, N. (2002), “Strategy for a rock mechanics site descriptive model, development and testing of the theoretical approach” Svensk Karnbranslehantering AB, Rapport R-02-02, May 2002  STEAD, D., COGGAN, J.S., EBERHARDT, E., (2004) “Realistic simulation of rock slope failure mechanisms: the need to incorporate principles of fracture mechanics”, Proceedings of SINOROCK international symposium on rock mechanics: rock characterization, modelling and engineering design methods, 18–21 May, Three Gorges Project site, China, Paper 2B-17  TANG, C.A. et al. (2000a) “Numerical studies of the influence of microstructure on rock failure in uniaxial compression - Part I: effect of heterogeneity”, International Journal of Rock Mechanics & Mining Sciences, Vol. 37, pp. 555-569  TANG, C.A. et al. (2000b) “Numerical studies of the influence of microstructure on rock failure in uniaxial compression - Part I: constraint, slenderness and size effect”, International Journal of Rock Mechanics & Mining Sciences, Vol. 37, pp. 571-583  TANG, C.A., (1995) “Numerical simulation of rock failure process”, Proceedings of the 2nd Youth Symposium on Rock Mechanics in China, Chengdu  TANG, C.A., (1997) “Numerical simulation on progressive failure leading to collapse and associated seismicity”, International Journal of Rock Mechanics & Mining Sciences, Vol. 34, pp. 249-261  127  TANG, C.A., KAISER, P.K., (1998) “Numerical simulation of cumulative damage and seismic energy release in unstable failure of brittle rock - Part I: Fundamentals”, International Journal of Rock Mechanics & Mining Sciences, Vol. 34, pp. 249-261  TERZAGHI, K. (1946), “Introduction to Tunnel Geology”, In Rock Tunneling with Steel Supports, pp. 47-99  VYAZMENSKY, A., ELMO, D., STEAD, D., (2009), “Role of Rock Mass Fabric and Faulting in the Development of Block Caving Induced Surface Subsidence”, Springer-Verlag, 24 p.   WARBURTON, P.M. (1980) “Stereological interpretation of joint trace data”, International Journal of Rock Mechanics, Mining Sciences and Geomechanics Abstracts, Vol. 17, pp. 181–190  WICKHAM, G.E., TIEDEMAN, H.R. and SKINNER, E.H., (1972), “Support determination based on geological predictions”, In: Proceedings of the North American Rapid Excavation & Tunnelling Conference, American Society of Mechanical Engineers, New York, pp. 43-64  WILLIAMSON, D.A., and KUHN, C.R., (1988), “The unified rock classification system”, Rock Engineering Systems for Engineering Purposes, ASTM STP 984, Louis Kirkaldie, Ed., American Society for Testing Materials, Philadelphia, pp. 7-16  WOODWORTH, J.B. (1897), “On the Fracture Systems of Joints, with Remarks on Certain Great Fractures,” Boston Society of Natural History Proceedings, Vol. 27, pp. 163-184     128  Appendix A: DFN Models Review  A1 Orthogonal Model The orthogonal model, as described by Snow, is characterised by two or three orthogonal sets of parallel unbounded joints, with constant spacing (Snow, 1965).  However, a Poisson distribution can be used to locate the joints in each set.  This assumption is the most commonly used for DFN (Dershowitz and Einstein, 1977).  In this case, the mean and the variance for the joint spacing can be described by the following equations:  2/1)(/1)(λλ==xVxE  Where, - x is a random variable; and - λ is the intensity. Furthermore, instead of being ubiquitous, joints can defined as bounded, either by assuming termination of joints at joint plane intersection or by assuming joint termination independently of joint planes intersections.  A2 Beacher Model In the Baecher model, the center points non-ubiquitous discs (joints) are randomly and independently distributed in space forming a Poisson field, thus generating Poisson lines in a 2D trace plan.  The radii of the joints follow a given distribution, often lognormal.   A3 Enhanced Beacher Model The Enhanced Baecher model is based on the Baecher model but, it accounts for fracture terminations at intersections with pre-existing fractures and for more general fracture shapes (Staub et al., 2002).  Fractures generated have generally three to six sides and termination probability is specified, i.e. the probability that a fracture of that set will end when it intersects a fracture from a pMasterrevious set.  The Enhanced Baecher offers a better simulation of the connectivity of natural fracture networks.  129  A4 Veneziano Model The Veneziano model uses Poisson lines on the Poisson planes to define polygonal areas.  A determined proportion of those polygons are then defined as open joints while the rest of the polygons are labelled as intact rock.  Bounded joints on joint planes are thus produced and can be used to model different rock mass conditions.  This model has been used in 2D for slope stability and hydrology in jointed rock.  A5 Dershowitz Model The Dershowitz model is based on the Veneziano model but instead of defining the polygons on the Poisson planes with Poisson lines, the polygons are defined by the intersections of the Poisson planes themselves.  As a result, all joint intersections occur at joint edges.  One of the advantages of the Dershowitz model over the Veneziano model is that it is simpler to represent and analyse.  A5 Stochastic Mosaic Block Tessellation Model In the stochastic mosaic block tessellation, points are generated in space by a Poisson process.  They can either be used to define the polyhedral centers, like in the Voronoi tessellation, or represent the block vertices, like in the Delaunay tessellation.  The faces of the different block thus generated can then either be defined as open features, joints, or intact rock.  The Orthogonal, Veneziano and Dershowitz models can all be considered as special cases of Mosaic tessellation.  A6 BART Model The BART model is also based on the Baecher model, but in this case termination is assigned by termination percentage.  This method represents more accurately field data (Staub et al., 2002).  In the BART model, locations of secondary fractures are controlled by the location of the primary fractures.  It results in spatially correlated and often clustered or chained fracture populations.    130  A7 3D Hierarchical Fracture Model The 3D Hierarchical Fracture model is based on Poisson plane and line processes.  Fractures are modelled as convex polygons that are randomly oriented and randomly located in space.  Three different stochastic processes are required for this model: the generation of Poisson planes to create fractures along planes of maximum shear and tension in rocks, a Poisson line tessellation and polygon marking to divide fracture planes in open features, joints, and intact rock and a random polygon translation and rotation to account for folding and faulting (Staub et al., 2002).  The 3D Hierarchical Fracture model enables the representation of relationship between fracture geometry and underlying mechanisms.  A8 Nearest NeighbourModel The Nearest Neighbour model is similar to the Enhanced Baecher model previously described except for its assumption about spatial distribution of fractures.  This model reflects the tendency of fractures to be clustered around major structures, for example faults.  Fractures are generated in three groups: primary, secondary and tertiary.  Location of joints from the primary group dictates the location of joints from the secondary and tertiary groups (Staub et al., 2002).  The following equation represents the probability of occurrence of a fracture at point x in space:  bx CLxP−=)(  Where,      - L is the distance from point x to the nearest fracture of a previous group;      - b is an empirical constant; and      - C is an empirical constant  A9 War Zone Model The “War Zone” model has been conceived to represent shear zones.  The areas comprised between two subparallel joints from a primary set, referred to as “war zones”, are identified and assigned a higher density function (Staub et al., 2002).  Apart from this assumption, everything else in the “War Zone” model is identical to the Enhanced Baecher model.   131  A10 Non-Planar Zone Model  The Non-Planar Zone model is a semi-stochastic model which generates fractures with location and orientation varying to a user defined non-planar surface (Staub et al., 2002).  A11 Levy-Lee Model The Levy-Lee model is a fractal model where joint centers are created sequentially by a Levy flight process and size of a joint is related to its distance from the previous fracture.     The Levy Flight process is a type of random walk.  The following probability equation gives the length L’ of each step (Staub et al., 2002).  DL LLLP−=> )'(   Where,      - D is the fractal dimension of the point field formed by the fracture centres The Levy-Lee model generates cluster of smaller fractures around bigger, more scattered, fractures.  A12 Fractal POCS Model The Fractal POCS model is a general procedure of generation of stochastic fields according to constraints such as spatial correlation which can be represented mathematically as convex sets (Staub et al., 2002).  A13 Poisson Rectangle Model The Poisson rectangle model is a special case of the Enhanced Baecher model where joints are represented by rectangles with prescribed length and width (Staub et al., 2002).  A14 MIT Model Meyer and Einstein introduced the MIT model in 2002, the latest step in the work on geometric-mechanical modeling started by Ivanova (Meyer and Einstein, 2002).  It allows the user to model joint systems in faulted geologies and makes it possible to represent geometric connectivity of the fractures.  In contrast to purely stochastic model, the MIT model imitates the natural 132  geomechanical process that generated the fractures.  Four stochastic process are necessary to generate the fractures in this model: the primary process models the orientation of the potential fracture planes according to an orientation distribution, with homogeneous, anisotropic Poisson process, the secondary process uses Poisson line tessellation to create polygon which are then divided into open features, joints, and intact rock, the tertiary process defines zone within the model and either retain or discard polygon located in those zones based on a certain probability that can vary from zone to zone and finally, the quaternary process translates and rotates remaining polygons to account for folding and faulting (Meyer and Einstein, 2002).        133  Appendix B: Rock Mass Rating Systems Review  B1 Rock load theory The rock load theory is the earliest system to use a rock mass classification for engineering purposes (Hoek, 2000).  It has been conceived by Terzaghi (1946) based on his experience with steel supported railroad tunnels in the Alps.  He divided rock masses into 7 groups: intact, stratified, moderately jointed, blocky and seamy, crushed, squeezing and swelling (Terzaghi, 1946).  A more thorough description of the 7 groups is presented in Terzaghi (1946)  According to Terzaghi, each one of these rock masses reacts differently to stress changes so they were consequently assigned different rock load.  They were estimated by studying the failure of wooden blocks of known strengths that were used for blocking the steel arches to the surrounding rock masses (Edelbro, 2003).  As this technique was developed for tunnels having steel arches as their only support element, it is not suitable for modern tunneling.  B2 Stand up time and NATM Lauffer (1958) suggested that the stand-up time of the unsupported span is related to rock mass quality in which the span is excavated (Hoek, 2000).  The unsupported span of the tunnel is defined as the distance between the excavation face and the nearest support if it is greater than the tunnel span. Lauffer’s original classification has been modified through the years and was the precursor of the New Austrian Tunneling Method (NATM) proposed by Rabcewicz (1964/65).    According to Lauffer, stand-up time significantly decreases with the increase of the tunnel span.  For example, a small pilot tunnel can be stable with minimal support whereas a larger tunnel would require the immediate installation of substantial support to achieve the same effect.  The concept at the base of the NATM is the conservation and mobilization of the strength of the soil or rock and the formation of largely self-supporting ring of soil or rock around the tunnel (Brown, 1981).  To accomplish that, primary support, usually shotcrete, is rapidly placed to help the rock mass preserve its load-carrying capacity. Deformation of the excavation and the build-up of load in the support or reinforcement are also monitored to check on the performance and the safety of the tunnel and to guide the provision of secondary and tertiary support.  Finally, the invert must be closed to form a load-bearing support ring to control deformation of the rock 134  mass.  The NATM method is also often associated with sequential excavation to maximize the stand-up time of the rock mass, hence allowing support to be installed.  B3 Rock Quality Designation (RQD) The Rock Quality Designation system was introduced by Deere (1964, 1968) to quantify discontinuity spacing.  RQD is a core recovery percentage indirectly based on the number of fractures and the amount of softening in the rock mass (Edelbro, 2003).  Only intact pieces longer than 10 cm are summed and divided by the total length of core run.  The core should be drilled as a double-tube core barrel and should be at least NW size.  Great care should be taken to ensure that drilling induced fractures are identified and ignored when determining the RQD.  Another way to determine the RQD, if core logs are not available but discontinuity traces are visible at the surface, is to estimate RQD from the number of discontinuities per unit volume (Palström, 1982).  The suggested relation for clay-free rock masses is:  g1844g1843g1830 = 115−3.3g1836g3049  where Jv (volumetric joint count) is the sum of the number of joints per unit length for all sets.  Drill core RQD is highly dependent on the direction.  The use of volumetric joint count can help reduce the dependence significantly.  B4 A recommended rock classification for rock mechanical purposes  Patching and Coates (1968) developed their classification system by modifying previous classifications by Coates and Coates and Parsons (1966).  The goal was to create a classification that would be easy to use but would still be detailed enough to divide rock masses into different categories.   B5 The Unified Classification of Soils and Rocks  The unified classification of soils and rocks by Deere and Miller (1969) is based on the testing of 257 NX-size rock cores from 27 localities.  Laboratory tests were conducted on these samples as follow: unit weight, Shore scleroscope hardness, Schmidt hammer hardness, abrasion hardness, 135  absorption, sonic-velocity stress-strain to failure and point-load tensile strength (Deere and Miller, 1969).  This classification is intended for intact rock.  As a result, it does not take into consideration any of the joint properties.  The strength classes are based on previous work by Coates and Patching.  Five charts have been produced to correlate strength or modulus properties for intact rock with either the Schmidt hammer, the Shore scleroscope or the sonic pulse velocity data.  B6 Rock Structure Rating (RSR) The RSR system, developed by Wickham et al. (1972), is the first classification to evaluate multiple components of the rock mass to arrive to a numerical value.  The 3 components are called A, B and C.  A accounts for the generic rock type with an index value and the general type of structure in the rock mass.  B considers the joint pattern in function of the direction of drive.  C combines the effect of groundwater inflow and joint condition.  The RSR value is a numerical value comprised between 0 and 100 and is calculated by adding the scores of A, B and C together.  The higher the RSR, the less support is needed for a same size tunnel.  The U.S. bureau of mine further developed the technique and identified 6 factors as being the most crucial for the prediction of support requirement (Edelbro, 2003).  Those 6 factors are:  1. Rock type with a strength index  A (maximum=30) 2. Geologic structure  3. Rock Joint spacing  B (maximum=45) 4. Orientation with respect to tunnel drive  5. Joint condition  C (maximum=25) 6. Groundwater inflow  The RSR system is based on 53 tunnel projects (322 km), mostly small to normal size (5.3-11m) and supported by steel sets, studies in eleven mines in western USA and case histories given by Proctor et al. (1946) of 183 tunnels and 64 mine examples (Edelbro, 2003).  The RSR system is 136  based on drill and blast tunnels, but can be adapted to excavation executed with a tunnel-boring machine by adding an added factor based on tunnel size.  B7 Rock Mass Rating (RMR) Bieniawski presented his first version of the RMR in 1973 based on his experience in the design of shallow tunnels in sedimentary rocks.  At first, only 49 unpublished case histories were used to conceive the system.  RMR has since been updated in 1974, 1975, 1976, 1979 and 1989 (Edelbro, 2003) more than trippling the number of case studies involved in the conception of the system.  As a result of the addition of those new case records, class boundaries have been changed and parameters have been reduced from 8 to 6 since the original version.  Those six parameters are:  1. Uniaxial Compressive Strength (UCS) of intact rock material 2. RQD 3. Spacing of discontinuities 4. Condition of discontinuities 5. Groundwater conditions 6. Orientation of discontinuities  RMR is obtained by adding the numerical value of the six parameters.  A second RMR value called RMRbasic can be obtained by adding only the 5 first parameters.  The sixth parameter is treated separately because it depends upon engineering application and not rock mass properties.     RMR has been applied to more than 351 projects such as tunnels, chambers, mines, slopes, foundations and rock cavern (Bieniawski, 1993).  The first step in characterizing the rock mass with the RMR system is to divide it into zones of similar geotechnical properties.  Such zones are often defined by major structures such as faults or change in rock type although significant change in discontinuity properties within the same unit can justify the division of the rock mass into a number of smaller structural sections.  137  Together with his 1989 version of the RMR, Bieniawski (1989) also published guidelines for the selection of the proper excavation sequence and support system.  Those guidelines are based on a 10 m span horseshoe-shaped tunnel excavated using drill and blast method and subjected to a vertical of less than 25 MPa (approximately equivalent to a depth of 900 m below surface).  It should also be noted that those guidelines have not been updated since 1973.  Therefore, new, more efficient, support method may exist now, like for example steel fiber reinforced shotcrete instead wire mesh and shotcrete.  B8 Q-System On the basis of the analysis of 200 tunnel case records, Barton et al. (1974), proposed a rock mass quality Q for the determination of tunnel support requirement.  The numerical values for Q range from 0.001 for exceptionally poor quality squeezing-ground up to 1000 for exceptionally good quality, practically unjointed, rock. (Barton et al., 1974).  Q is a function of six parameters:   1. RQD; 2. Jn, the number of joint sets; 3. Jr the roughness of the weakest joint; 4. Ja the degree of alteration or filling along the weakest joints; 5. Jw the joint water reduction factor; and 6. SRF the stress reduction factor.  and is defined as:  g1843 = g3436g1844g1843g1830g1836g3041g3440×g3436g1836g3045g1836g3028g3440×g3436 g1836g3028g1845g1844g1832g3440    According to Barton et al. the rock mass quality Q can be considered a function of the following three parameters:  1. Block size (RQD/Jn) 2. Inter-block shear strength (Jr/Ja) 138  3. Active stress (Jw/SRF)  In order to relate the index Q to the stability and support requirements, Barton et al. (1974), suggested the use of an additional parameter called the Equivalent Dimension (De).  De is defined as: g1830g3032 = g1831g1876g1855g1853g1874g1853g1872g1861g1867g1866 g1871g1868g1853g1866,g1856g1861g1853g1865g1857g1872g1857g1870 g1867g1870 planckg1857g1861g1859planckg1872 (g1865)g1831g1845g1844   The Excavation Support Ratio refelects the intended use of the excavation and the degree of safety and support demanded.  The Q system has its best applications in jointed rock masses where instability is caused by block falls (Palmström and Stille, 2006).   B9 Mining RMR (MRMR) The MRMR was introduced by Laubscher (1975) and modified by Laubscher and Taylor (1976) as a development of the CSIR Geomechanics classification system (Edelbro, 2003).  The fundamental difference between RMR and MRMR is the adjustment of the in-situ ratings (RMR) to account for the mining environment so the final ratings (MRMR) can be used for mine design (Laubscher, 1990).  The adjustment parameters are: weathering, mining-induced stresses, joint orientation, and blasting effects.  MRMR values range from 0 to 100 and are divided in five classes and ten sub-classes.  A set of support recommendations is suggested for the final and adjusted RMR value.  B10 The Unified Rock Classification System (URCS) The URCS was developed by Williamson in 1959 and 1960 (Williamson and Kuhn, 1988).  It has since been used on projects for the U.S. Army Corps of Engineers, USDA Soil Conservation Service, and the USDA Forest Service (Williamson and Kuhn, 1988).  According to Williamson and Kuhn (1988), the URCS has proven especially helpful to those involve in blasting design because it provides them with the necessary parameters for the design of a shot (Williamson and Kuhn, 1988).  The four fundamental physical properties at the used in the URCS are: 139   1. Weathering; 2. Strength; 3. Discontinuity; and 4. Density  Each of those properties is divided in five categories ranging from A to E, establishing the limiting values of each of the four basic elements.  A rock mass described as AAAAA would be requiring the least design, while a rock mass described as EEEEE would require the most.  The equipment used for the field tests is simple and available: your fingers, 10-power hand lens, a 1 lb ball peen hammer, and spring-loaded scales with 10 lb range.   B11 Basic Geotechnical Description (BGD) The BGD was established in 1981 by the International Society for Rock Mechanics (ISRM) to meet three major requirements:  1. To provide a language enabling the observer to transmit his general impression of a rock mass in a non-ambiguous fashion.  Different observers of a given rock mass should describe it in the same way; 2. To contain as much as possible quantitative data of interest to the solution of definite problems; and 3. Whenever possible, to use simple measurements, rather than visual observation alone. (ISRM)   140  Before applying the BGD to a rock mass, one should first divide it in geotechnical units having uniform characteristics with regard to the project requirements.  The BGD value is then assessed by: 1. The rock name and simplified geological description; 2. The layer thickness; 3. The fracture intercept; 4. The uniaxial compressive strength of the rock material; and 5. The angle of friction of the fractures.  B12 Rock Mass Strength (RMS) The RMS system was developed by Stille et al. (1982) based on only 8 different cases where the rock mass strength was measured or estimated from different field observations.  RMS is a strength classification derived from the RMR system as it includes the first five parameters of RMRbasic (loading conditions and initial stress are not considered) (Edelbro, 2003).    The RMS value is composed of the RMRbasic value plus a rating reduction factor to take the presence of joints and joint sets into consideration.  B13 Modified Basic RMR (MBR) The MBR system was introduced in 1982 by Cummings et al. (1982) for block caving operations in the Western United States.  It has been developed from data collected from horizontal drifts in block, panel or mass caving mines so its applicability to vertical structure is uncertain.  MBR is a modified RMR for mining applications; consequently it uses many of the same input parameters.  MBR involves different ratings for the original RMR parameters and subsequent adjustments to allow for blast damage, induced stresses, structural features, distance from the cave front and size of the caving block (Hoek, 2000).  Guidelines for support are provided by support charts and tables (see Kendroski et al., 1983).   141  B14 Simplified Rock Mass Raring (SRMR) The SRMRM is born from a study by Brook and Dharmarante (1985) of the application of Q, RMR and MRMR.  They found the later to be the most successful system and decided to develop a simpler version with only three major components:   1. The intact rock strength, 2. The joint spacing, and 3. The joint type  together with a consideration for groundwater.  In order to account for environmental conditions, Brok and Dharmarante suggest the application of adjustment factors.    The SRMR system is based on data from two graphite mines and one quartz-gneiss tunnel.   Brook and Dharmarante intentionally omitted the RQD with the purpose of simplifying their system and avoiding the difficulties associated with its measurement.   B15 Slope Mass Rating (SMR) The Slope Mass Rating system was introduced in 1985 by Romana in an effort to adjust the RMR system to slope design.  Romana suggests adding adjustment factors depending on the relation between the slope and joints and the excavation method to the RMRbasic value:  g1845g1839g1844 = g1844g1839g1844g3029g3028g3046g3036g3030 +(g1832g2869 ∙g1832g2870 ∙g1832g2871)+g1832g2872  where F1 is given by:   g1832g2869 = (1−g1871g1861g1866g1827)g2870  where A is the angle between the strike of the slope face and strike of the joint.  F2 depends on the joint dip angle of failure in the planar mode and F3 refers to the relationship between the slope face and joint dip.  F4 is the adjustment factor depending on the excavation method of the slope (Edelbro, 2003). 142   When introduced in 1985, SMR was first based on 28 slopes, both natural and excavated.  The system was then applied to 44 slopes in 1988 and later to slopes in a quarry (Edelbro, 2003).  B16 Slope Rock Mass Rating The Slope Rock Mass Rating system was developed by Robertson (1988) to make CSIR’s RMR system consistent for both weak and strong rock masses. Three major differences exist between RMR and the Slope Rock Mass Rating.  The groundwater parameter has been dropped because, according to Robertson, it is a destabilizing force and should be accounted for as such in stability analysis not in the rock mass strength. The RQD is replaced by the Handeled RQD (HRQD) which is measured in the same way as RQD but after the core has been firmly handled and “worried” in an attempt to break the core into smaller fragments (Robertson, 1988). And the rating for spacing of discontinuities is determined from the handled core.  B17 Ramamurthy and Arora Classification (RAC) The RAC for intact and jointed rock was introduced by Ramamurthy and Arora (1993) because they are of the opinion that neither strength nor modulus alone can represent the rock mass quality of the rock.  The classification is based on compressive strengths and modulus values in their unconfined state.  The modulus ratio Mrj, used to determine the class, is given by:  g1839g3045g3037 = g1831g3047g3037g2026g3030g3037= 1g2035g3033  where j stands for jointed rock (it can also be replaced by I for intact rock) and Et is the tangent modulus at 50% of the failure stress.  Values for σcj and Etj are obtained from the following equations:  g2026g3030g3037g2026g3030g3036 = expg3435−0.008∙g1836g3033g3439  g1831g3047g3037g1831g3047g3036 = exp g3435−0.0115∙g1836g3033g3439 143   where Jf is the joint factor and is given by:  g1836g3033 = g1836g3041g1866 ∙g1872g1853g1866   where Jn is the joint frequency, and the inclination parameter (n) depends on the inclination of the joint plane with respect to the major principal stress (Edelbro, 2003).    B18 Geological Strength Index (GSI) The Geological strength index was introduced by Hoek et al. (1995) to complement the Hoek Brown failure criterion, as a way to estimate s a and mb, and to overcome the fact that RMR is not suitable for very poor rock masses since the minimum value which RMR can assume is 18 (Hoek et al., 1995).  GSI values range from about 10 for very weak rock masses to 100 for intact rock.  Whenever using the GSI tables, it is very important to carry out the classification on an undisturbed rock mass (Hoek and Brown, 1997).   B19 Rock Mass Number (N) and Rock Condition Rating (RCR) Goel et al. (1996) proposed modifications to the RMR and Q systems in an attempt to find a better relation between both systems.  According to Goel et al. (1996) the major discrepancies between the two systems is that the RMR system does not consider the stress condition of the rock mass, while the Q-system does not consider joint orientation and intact rock strength as independent parameters (Goel et al., 1996).  In order to take those differences into consideration, a stress-free version of the Q system, the N-system, and a version of the RMR system without ratings for the compressive strength of the intact rock material and adjustments of joint orientation, the RCR system, have been developed.  g1840 = g1844g1843g1830g1836g3041∙ g1836g3045g1836g3028∙g1836g3050  g1844g1829g1844 = g1844g1839g1844 −(g1844g1853g1872g1861g1866g1859 g1858g1867g1870 g2026g3030 +g1827g1856g1862g1873g1871g1872g1865g1857g1866g1872 g1867g1858 g1862g1867g1861g1866g1872 g1867g1870g1861g1857g1866g1872g1853g1872g1861g1867g1866)  144  RCR and rock mass number N from 63 cases were used to obtain the following correlation between the two systems:  g1844g1829g1844 = 8∙ln (g1840)+30  Out of the 63 cases studied, 36 were from India, 4 from the Kielder experimental tunnel and 23 from NGI cases (Goel et al., 1996).  B20 Rock Mass Index (RMi) Palmström developed the RMi system in 1996 to characterize the strength of the rock mass for construction purposes (Palmström, 1996).  The system has been developed with data from triaxial laboratory tests on Panguna andesite, a large scale compressive laboratory test on granitic rock from Stripa, in-situ tests on mine pillars of sandstone in the Laisvall mine, strength data found from back-analysis of a slide in the Langsele mine and large-scale laboratory triaxial tests from Germany (Edelbro, 2003).    The RMi is essentially the reduce rock strength caused by jointing (Edelbro, 2003) and can be calculated with the following equation:  g1844g1839g1861 = g2026g3030 ∙g1836g1842  where JP is the jointing parameter, varying from 0 for crushed rock to 1 for intact rock, expressed by the following expression:  g1836g1842 = 0.2g3496g1862g1838 ∙g3436g1862g1844g1862g1827g3440∙g1848g3029g2868.g2871g2875∙g4672g3037g3013∙g3037g3019g3037g3002g4673g3127g3116.g3118  where jL is the size factor, jR the roughness factor, jA the alteration factor and Vb the block volume.  For more information on those parameters see Palmström, 1995 or Edelbro, 2003.  145  Appendix C: C++ Programming  To facilitate and accelerate the process of sampling the fracture properties to assign to the different regions of the models, a C++ program was created (see Appendix C).  The program is designed to find any fracture going through a cell of an imaginary grid superimposed over a DFN model.  It then selects, based on different criteria, the properties that will be inputted into each cell of a FLAC model.  Finally, the C++ program generates a FLAC .dat file with all the information necessary to run the model.  The input parameters for the program are the DFN fracture file, the surface file, the intact rock properties, the fracture properties and the number of cells desired for the grid.   /*-----------------------------------------------*/ /* FILE: FLAC3B.cpp                              */ /* AUTHOR:  Thierry Lavoie                       */ /* DATE: 06/08/2009                              */ /* DESCRIPTION: Generates FLAC file from data    */ /*              on the fractures and surfaces    */ /*              from a DFN model                 */ /*                                               */ /*-----------------------------------------------*/  #include <iostream>  //For cin et cout #include <iomanip>  //For cin et cout #include <fstream>   //For writing and reading files #include <ctime>   //For randomize() #include <cstring>   //For strcmp() #include <cstdlib>   //For srand() et rand() #include <cmath>   //For sin, cos, pow... const double PI=3.141592653589793238462643383279502884197169399375105820974944592; //PI using namespace std;  void main(void)  //Begining of the main function   {   /* Variables */  ifstream in, in2;  ofstream out, outh, outs;  int Big,H,i,t,q,s,x,z,y,w,Area,ColCell,RowCell;//i=nb fractures in file read, x count, ColCell nb of columns, RowCell number of rows  char InName[104],InName2[104], OutName[104],OutNameI[104], OutNameh[114], OutNameS[114];//InName InName2 name of files to read, name of files to create   double Frac [1000][13];// col 0: Frac Number, col 1: X1, col 2: Y1,col 3: X2, col 4: Y2, col 5: Xmax, col: 6 Ymax, Col 7: Xmin, col 8: Ymin, clo 9; y=Ax+b, col 10: y=ax+B, col 11: DIP, col 12: Frac Length 146   double buffer [1][15];// to read the file  double ***cell;//0:Young, 1:Poisson, 2:Cohesion, 3:Tension, 4:Friction, 5:Dilation,6:Jcohesion,7:Jfriction,8:Jdilation,9:Jtension,10:Jangle,11:Xmin,12Ymin,13Xmax,14Ymax,15:P21,16:P00,17:P20,18:P11*,19:P01,20:rand DIP, 21:rand Length, 22:Resultant DIP, 23:resultant length, 24: most represented DIP, 25:most represented length , 26: Critical DIP , 27: Critical length, 28: Critical*length Dip  char fill[50];// to read the file  char QGSI, Ubi, Grav, Damp;//GSI Yes or No, Ubiquitous joint model Y or N, Gravity Y or N, Combined Damping Y or N  double Young, Poisson, Cohesion, Tension, Friction, Dilation,Jcohesion,Jfriction,Jdilation,Jtension,Density; //Rock and Joint Properties  double BulkF, ShearF, CohesionF, TensionF, FrictionF, DilationF, DensityF;//joint properties (Mohr)  double Surface [6][6];//0:Num, 1:Xmax, 2:Ymax, 3:Xmin, 4:Ymin, 5:Area  double SizeCellx, SizeCelly;//size of cells  double fraclen;//Total Fracture length by cell  double random [5][30];//properties of the fractures in a cell 0 dip, 1 length, 2 nb of same orientation in cell, 3 length of same orientation, 4 total length  int V,P00,Choice,Most;//V for random number, P00 number of frac by cell, Choice between random, resultant, most represented or critical, Most number or length  double resulx, resuly;//for the resultant frac  double GSI,S,mb,UCS;//For the calculation of the rock mass properties  double conf, ini;//For the confinement pressure, initial velocity  int width, height;//Sample width and height  double criti;//For the critical fracture    cout<<"Name of FRACTURE file to read, no spaces or special caracters (max 100 car):";  cin>>InName;   strcat_s(InName,".txt");    in.open(InName);   if (in.fail())  {   cout<<"opening problem with: "<<InName<<endl<<endl;  }   else  {   cout<<"Name of SURFACE file to read, no spaces or special caracters (max 100 car):";   cin>>InName2;    strcat_s(InName2,".txt");      in2.open(InName2);    if (in2.fail())   {    cout<<"opening problem with: "<<InName2<<endl<<endl; 147    }    else   {    cout<<endl<<"Name of FLAC file to create, no sapces or special caracters (max 100 car):";    cin>>OutName;    strcpy_s(OutNameI,OutName);    strcat_s(OutName,".dat");    out.open(OutName);     if (out.fail())    {     cout<<"opening problem with: "<<OutName<<endl<<endl;     }    else    {     //initialization of parameters     Young=0;     Poisson=0;     Cohesion=0;     Tension=0;     Friction=0;     Dilation=0;     Jcohesion=0;     Jfriction=0;     Jdilation=0;     Jtension=0;      //initialization of the seed     srand((unsigned)time(NULL));     cout<<endl<<endl;     cout<<"Enter confinement pressure (Pa)";     cin>>conf;     cout<<endl;     cout<<"Enter initial velocity (mm/time step)";     cin>>ini;     cout<<endl;      //Rock and Joint properties     cout<<endl<<endl;     cout<<"Enter the properties of the different parameters"<<endl<<endl<<"INTACT ROCK"<<endl;     cout<<"Young's Modulus [Pa]:";     cin>>Young;     cout<<endl;     cout<<"Poisson's Ratio:";     cin>>Poisson;     cout<<endl;     cout<<"Cohesion [Pa]:";     cin>>Cohesion;     cout<<endl;     cout<<"Tensile Strength[Pa]:";     cin>>Tension;     cout<<endl;     cout<<"UCS [Pa]:"; 148      cin>>UCS;     cout<<endl;     cout<<"Friction angle:";     cin>>Friction;     cout<<endl;     cout<<"Dilation angle:";     cin>>Dilation;     cout<<endl;     cout<<"Density(kg/m3):";     cin>>Density;     cout<<endl<<endl;       cout<<"Do you want to use the Ubiquitous Joint Model? (Y/N)";     cin>>Ubi;     cout<<endl;     cout<<endl;             if (Ubi=='Y' || Ubi=='y')     {      //Choice for the way too get the angle of the fracture for the ubiquitous joint model      cout<<"1) Random;"<<endl<<"2) Resultant;"<<endl<<"3) Most Represented;"<<endl<<"4) Critical; or"<<endl<<"5) Critical*Length;";      cin>>Choice;      cout<<endl;      cout<<endl;       if (Choice==3)      {       cout<<"1) Number; or"<<endl<<"2) Length.";       cin>>Most;       cout<<endl;       cout<<endl;      }       else      {}       cout<<"JOINT PROPERTIES"<<endl;      cout<<"Joint Cohesion [Pa]:";      cin>>Jcohesion;      cout<<endl;      cout<<"Joint Tensile Strength[Pa]:";      cin>>Jtension;      cout<<endl;      cout<<"Joint Friction angle:";      cin>>Jfriction;      cout<<endl;      cout<<"Joint Dilation angle:";      cin>>Jdilation;  149       cout<<endl<<endl;     }      else     {      cout<<endl<<endl;      cout<<"Enter the properties of the different parameters"<<endl<<endl<<"FRACTURES"<<endl;      cout<<"Bulk Modulus [Pa]:";      cin>>BulkF;      cout<<endl;      cout<<"Shear Modulus:";      cin>>ShearF;      cout<<endl;      cout<<"Cohesion [Pa]:";      cin>>CohesionF;      cout<<endl;      cout<<"Tensile Strength[Pa]:";      cin>>TensionF;      cout<<endl;      cout<<"Friction angle:";      cin>>FrictionF;      cout<<endl;      cout<<"Dilation angle:";      cin>>DilationF;      cout<<endl;      cout<<"Density(kg/m3):";      cin>>DensityF;      cout<<endl<<endl;     }      cout<<"Do you want to use GSI? (Y/N)";     cin>>QGSI;     cout<<endl;     cout<<endl;      cout<<"Do you want to use Gravity? (Y/N)";     cin>>Grav;     cout<<endl;     cout<<endl;       cout<<"Do you want to use Combined Damping? (Y/N)";     cin>>Damp;     cout<<endl;     cout<<endl;      cout<<"Enter the dimensions of the model:"<<endl;     cout<<endl;     cout<<"Width:";     cin>>width;     cout<<endl;     cout<<"Height:";     cin>>height;     cout<<endl<<endl;  150      cout<<"Enter the number of cells for the model (positive integer divisible by (Height/Width)"<<endl;     cout<<endl;     cout<<"Rows:";     cin>>RowCell;     cout<<endl<<endl;          //reads surface file     in2>>fill;     in2>>fill;     in2>>fill;     in2>>fill;     in2>>fill;     in2>>fill;          //Reads surface file     y=0;     Area=0;     while( ! in2.eof() )     {      x=0;      while (x<6)      {       in2>>Surface[y][x];       x=x+1;      }      //Finds the biggest surface      if (Surface[y][5]>=Surface[Area][5])      {       Area=y;      }      else      {}      y=y+1;     }          //Number of columns     H=Surface[Area][2]/Surface[Area][1];     ColCell=RowCell/H;          //calculates size of cells     SizeCellx=Surface[Area][1]/ColCell;     SizeCelly=Surface[Area][2]/RowCell;      //dynamic allocation for the cells     cell=new double** [RowCell];     x=0;     while (x<RowCell)     {      cell[x]=new double* [ColCell];      z=0;      while (z<ColCell)      {       cell[x][z]= new double [29];       z=z+1;      }      x=x+1; 151      }      //coordinates Xmin, Ymin, Xmax, Ymax of the first cell     cell[0][0][11]=0;     cell[0][0][12]=0;     cell[0][0][13]=SizeCellx;     cell[0][0][14]=SizeCelly;          //coordinates Xmin, Ymin, Xmax, Ymax of all other cells     x=0;     while (x<RowCell)     {      if (x!=0)      {       cell[x][0][11]=0;       cell[x][0][12]=cell[x-1][0][14];       cell[x][0][13]=SizeCellx;       cell[x][0][14]=cell[x][0][12]+SizeCelly;      }      else      {}      z=1;      while (z<ColCell)      {       cell[x][z][11]=cell[x][z-1][13];       cell[x][z][12]=cell[x][z-1][12];       cell[x][z][13]=cell[x][z-1][13]+SizeCellx;       cell[x][z][14]=cell[x][z-1][14];       z=z+1;      }      x=x+1;     }          //reads the fracture file     i=0;     x=0;     while (x<17)     {      in>>fill;      x=x+1;     }     //reads the fracture file     while( ! in.eof() )     {      x=0;      while (x<15)      {       in>>buffer[0][x];       x=x+1;      }            //transfers fractures only, not surfaces lines, into Frac[][]      if (buffer [0][14]==0) 152       {       Frac[i][0]=i+1;       Frac[i][1]=buffer[0][1];       Frac[i][2]=buffer[0][2];       Frac[i][3]=buffer[0][4];       Frac[i][4]=buffer[0][5];       Frac[i][5]=buffer[0][8];       Frac[i][6]=buffer[0][9];       Frac[i][7]=buffer[0][11];       Frac[i][8]=buffer[0][12];       Frac[i][12]=buffer[0][7];        //Calculates A       //Checks if x1<x2       if (Frac[i][1]<Frac[i][3])       {        Frac[i][9]=(Frac[i][4]-Frac[i][2])/(Frac[i][3]-Frac[i][1]);       }       else       {        Frac[i][9]=(Frac[i][2]-Frac[i][4])/(Frac[i][1]-Frac[i][3]);       }         Frac[i][11]=(atan(Frac[i][9]))*180/PI;         //Calculates B and DIP       Frac[i][10]=Frac[i][2]-Frac[i][9]*Frac[i][1];        i=i+1;             }      else {}          }           // Calculates P's (P21, P20...), also determinates random and resultant frac     x=0;     while (x<RowCell)     {      z=0;      while (z<ColCell)      {       w=0;       fraclen=0;       P00=0;       while (w<i)       {                //Checks if x1 y1, and x2,y2 are in the cell 153         if (((Frac[w][1]<=cell[x][z][13]+0.0001) && (Frac[w][1]>=cell[x][z][11]-0.0001) && (Frac[w][2]<=cell[x][z][14]+0.0001) && (Frac[w][6]>=cell[x][z][12]-0.0001)) && ((Frac[w][3]<=cell[x][z][13]+0.0001) && (Frac[w][3]>=cell[x][z][11]-0.0001) && (Frac[w][4]<=cell[x][z][14]+0.0001) && (Frac[w][4]>=cell[x][z][12]-0.0001)))        {         P00=P00+1;         fraclen=fraclen+sqrt(pow(Frac[w][4]-Frac[w][2],2)+pow(Frac[w][3]-Frac[w][1],2));         random[0][P00-1]=Frac[w][11];         random[1][P00-1]=sqrt(pow(Frac[w][4]-Frac[w][2],2)+pow(Frac[w][3]-Frac[w][1],2));         random[4][P00-1]=Frac[w][12];        }        //checks if X1,Y1 only is in the cell        else if ((Frac[w][1]<cell[x][z][13]) && (Frac[w][1]>cell[x][z][11]) && (Frac[w][2]<cell[x][z][14]) && (Frac[w][2]>cell[x][z][12]))        {         P00=P00+1;         //checks if x2<=x1 and y2<=y1         if((Frac[w][3]<=Frac[w][1]+0.0001) && (Frac[w][4]<=Frac[w][2]+0.0001))         {          //checks if fracture goes through the side of the cell          if (((Frac[w][9]*cell[x][z][11]+Frac[w][10])>=cell[x][z][12]-0.0001 && (Frac[w][9]*cell[x][z][11]+Frac[w][10])<=cell[x][z][14]+0.0001))          {           fraclen=fraclen+sqrt(pow((Frac[w][2]-(Frac[w][9]*cell[x][z][11]+Frac[w][10])),2)+pow(Frac[w][1]-cell[x][z][11],2));           random[0][P00-1]=Frac[w][11];           random[1][P00-1]=sqrt(pow((Frac[w][2]-(Frac[w][9]*cell[x][z][11]+Frac[w][10])),2)+pow(Frac[w][1]-cell[x][z][11],2));           random[4][P00-1]=Frac[w][12];          }           //Fracture goes through the bottom of the cell          else          {           fraclen=fraclen+sqrt(pow(Frac[w][1]-((cell[x][z][12]-Frac[w][10])/Frac[w][9]),2)+pow(Frac[w][2]-cell[x][z][12],2));           random[0][P00-1]=Frac[w][11]; 154            random[1][P00-1]=sqrt(pow(Frac[w][1]-((cell[x][z][12]-Frac[w][10])/Frac[w][9]),2)+pow(Frac[w][2]-cell[x][z][12],2));           random[4][P00-1]=Frac[w][12];          }         }          //checks if x2>x1 and y2<y1         else if((Frac[w][3]>=Frac[w][1]-0.0001) && (Frac[w][4]<=Frac[w][2]+0.0001))         {          //checks if fracture goes through the side of the cell          if (((Frac[w][9]*cell[x][z][13]+Frac[w][10])>=cell[x][z][12]-0.0001 && (Frac[w][9]*cell[x][z][13]+Frac[w][10])<=cell[x][z][14]+0.0001))          {           fraclen=fraclen+sqrt(pow((Frac[w][2]-(Frac[w][9]*cell[x][z][13]+Frac[w][10])),2)+pow(Frac[w][1]-cell[x][z][13],2));           random[0][P00-1]=Frac[w][11];           random[1][P00-1]=sqrt(pow((Frac[w][2]-(Frac[w][9]*cell[x][z][13]+Frac[w][10])),2)+pow(Frac[w][1]-cell[x][z][13],2));           random[4][P00-1]=Frac[w][12];          }           //Fracture goes through the bottom of the cell          else          {           fraclen=fraclen+sqrt(pow(Frac[w][1]-((cell[x][z][12]-Frac[w][10])/Frac[w][9]),2)+pow(Frac[w][2]-cell[x][z][12],2));           random[0][P00-1]=Frac[w][11];           random[1][P00-1]=sqrt(pow(Frac[w][1]-((cell[x][z][12]-Frac[w][10])/Frac[w][9]),2)+pow(Frac[w][2]-cell[x][z][12],2));           random[4][P00-1]=Frac[w][12];          }         }          //checks if x2>x1 and y2>y1         else if((Frac[w][3]>=Frac[w][1]-0.0001) && (Frac[w][4]>=Frac[w][2]-0.0001))         {          //checks if fracture goes through the side of the cell 155           if (((Frac[w][9]*cell[x][z][13]+Frac[w][10])>=cell[x][z][12]-0.0001 && (Frac[w][9]*cell[x][z][13]+Frac[w][10])<=cell[x][z][14]+0.0001))          {           fraclen=fraclen+sqrt(pow((Frac[w][2]-(Frac[w][9]*cell[x][z][13]+Frac[w][10])),2)+pow(Frac[w][1]-cell[x][z][13],2));           random[0][P00-1]=Frac[w][11];           random[1][P00-1]=sqrt(pow((Frac[w][2]-(Frac[w][9]*cell[x][z][13]+Frac[w][10])),2)+pow(Frac[w][1]-cell[x][z][13],2));           random[4][P00-1]=Frac[w][12];          }           //Fracture goes through the top of the cell          else          {           fraclen=fraclen+sqrt(pow(Frac[w][1]-((cell[x][z][14]-Frac[w][10])/Frac[w][9]),2)+pow(Frac[w][2]-cell[x][z][14],2));           random[0][P00-1]=Frac[w][11];           random[1][P00-1]=sqrt(pow(Frac[w][1]-((cell[x][z][14]-Frac[w][10])/Frac[w][9]),2)+pow(Frac[w][2]-cell[x][z][14],2));           random[4][P00-1]=Frac[w][12];          }         }          //checks if x2<x1 and y2>y1         else         {          //checks if fracture goes through the side of the cell          if (((Frac[w][9]*cell[x][z][11]+Frac[w][10])>=cell[x][z][12]-0.0001 && (Frac[w][9]*cell[x][z][11]+Frac[w][10])<=cell[x][z][14]+0.0001))          {           fraclen=fraclen+sqrt(pow((Frac[w][2]-(Frac[w][9]*cell[x][z][11]+Frac[w][10])),2)+pow(Frac[w][1]-cell[x][z][11],2));           random[0][P00-1]=Frac[w][11];           random[1][P00-1]=sqrt(pow((Frac[w][2]-(Frac[w][9]*cell[x][z][11]+Frac[w][10])),2)+pow(Frac[w][1]-cell[x][z][11],2));           random[4][P00-1]=Frac[w][12];          } 156            //Fracture goes through the top of the cell          else          {           fraclen=fraclen+sqrt(pow(Frac[w][1]-((cell[x][z][14]-Frac[w][10])/Frac[w][9]),2)+pow(Frac[w][2]-cell[x][z][14],2));           random[0][P00-1]=Frac[w][11];           random[1][P00-1]=sqrt(pow(Frac[w][1]-((cell[x][z][14]-Frac[w][10])/Frac[w][9]),2)+pow(Frac[w][2]-cell[x][z][14],2));           random[4][P00-1]=Frac[w][12];          }         }        }               //checks if x2,y2 only is in the cell        else if ((Frac[w][3]<cell[x][z][13]) && (Frac[w][3]>cell[x][z][11]) && (Frac[w][4]<cell[x][z][14]) && (Frac[w][4]>cell[x][z][12]))        {         P00=P00+1;         //checks if x1<=x2 and y1<=y2         if((Frac[w][3]>=Frac[w][1]-0.0001) && (Frac[w][4]>=Frac[w][2]-0.0001))         {          //checks if fracture goes through the side of the cell          if (((Frac[w][9]*cell[x][z][11]+Frac[w][10])>=cell[x][z][12]-0.0001 && (Frac[w][9]*cell[x][z][11]+Frac[w][10])<=cell[x][z][14]+0.0001))          {           fraclen=fraclen+sqrt(pow((Frac[w][4]-(Frac[w][9]*cell[x][z][11]+Frac[w][10])),2)+pow(Frac[w][3]-cell[x][z][11],2));           random[0][P00-1]=Frac[w][11];           random[1][P00-1]=sqrt(pow((Frac[w][4]-(Frac[w][9]*cell[x][z][11]+Frac[w][10])),2)+pow(Frac[w][3]-cell[x][z][11],2));           random[4][P00-1]=Frac[w][12];          }           //Fracture goes through the bottom of the cell          else          {           fraclen=fraclen+sqrt(pow(Frac[w][3]-((cell[x][z][12]-Frac[w][10])/Frac[w][9]),2)+pow(Frac[w][4]-cell[x][z][12],2)); 157            random[0][P00-1]=Frac[w][11];           random[1][P00-1]=sqrt(pow(Frac[w][3]-((cell[x][z][12]-Frac[w][10])/Frac[w][9]),2)+pow(Frac[w][4]-cell[x][z][12],2));           random[4][P00-1]=Frac[w][12];          }         }          //checks if x2<x1 and y2>y1         else if((Frac[w][3]<=Frac[w][1]+0.0001) && (Frac[w][4]>=Frac[w][2]-0.0001))         {          //checks if fracture goes through the side of the cell          if (((Frac[w][9]*cell[x][z][13]+Frac[w][10])>=cell[x][z][12]-0.0001 && (Frac[w][9]*cell[x][z][13]+Frac[w][10])<=cell[x][z][14]+0.0001))          {           fraclen=fraclen+sqrt(pow((Frac[w][4]-(Frac[w][9]*cell[x][z][13]+Frac[w][10])),2)+pow(Frac[w][3]-cell[x][z][13],2));           random[0][P00-1]=Frac[w][11];           random[1][P00-1]=sqrt(pow((Frac[w][4]-(Frac[w][9]*cell[x][z][13]+Frac[w][10])),2)+pow(Frac[w][3]-cell[x][z][13],2));           random[4][P00-1]=Frac[w][12];          }           //Fracture goes through the bottom of the cell          else          {           fraclen=fraclen+sqrt(pow(Frac[w][3]-((cell[x][z][12]-Frac[w][10])/Frac[w][9]),2)+pow(Frac[w][4]-cell[x][z][12],2));           random[0][P00-1]=Frac[w][11];           random[1][P00-1]=sqrt(pow(Frac[w][3]-((cell[x][z][12]-Frac[w][10])/Frac[w][9]),2)+pow(Frac[w][4]-cell[x][z][12],2));           random[4][P00-1]=Frac[w][12];          }         }          //checks if x2<x1 and y2<y1         else if((Frac[w][3]<=Frac[w][1]+0.0001) && (Frac[w][4]<=Frac[w][2]+0.0001))         {          //checks if fracture goes through the side of the cell 158           if (((Frac[w][9]*cell[x][z][13]+Frac[w][10])>=cell[x][z][12]-0.0001 && (Frac[w][9]*cell[x][z][13]+Frac[w][10])<=cell[x][z][14]+0.0001))          {           fraclen=fraclen+sqrt(pow((Frac[w][4]-(Frac[w][9]*cell[x][z][13]+Frac[w][10])),2)+pow(Frac[w][3]-cell[x][z][13],2));           random[0][P00-1]=Frac[w][11];           random[1][P00-1]=sqrt(pow((Frac[w][4]-(Frac[w][9]*cell[x][z][13]+Frac[w][10])),2)+pow(Frac[w][3]-cell[x][z][13],2));           random[4][P00-1]=Frac[w][12];          }           //Fracture goes through the top of the cell          else          {           fraclen=fraclen+sqrt(pow(Frac[w][3]-((cell[x][z][14]-Frac[w][10])/Frac[w][9]),2)+pow(Frac[w][4]-cell[x][z][14],2));           random[0][P00-1]=Frac[w][11];           random[1][P00-1]=sqrt(pow(Frac[w][3]-((cell[x][z][14]-Frac[w][10])/Frac[w][9]),2)+pow(Frac[w][4]-cell[x][z][14],2));           random[4][P00-1]=Frac[w][12];          }         }          //checks if x2>x1 and y2<y1         else         {          //checks if fracture goes through the side of the cell          if (((Frac[w][9]*cell[x][z][11]+Frac[w][10])>=cell[x][z][12]-0.0001 && (Frac[w][9]*cell[x][z][11]+Frac[w][10])<=cell[x][z][14]+0.0001))          {           fraclen=fraclen+sqrt(pow((Frac[w][4]-(Frac[w][9]*cell[x][z][11]+Frac[w][10])),2)+pow(Frac[w][3]-cell[x][z][11],2));           random[0][P00-1]=Frac[w][11];           random[1][P00-1]=sqrt(pow((Frac[w][4]-(Frac[w][9]*cell[x][z][11]+Frac[w][10])),2)+pow(Frac[w][3]-cell[x][z][11],2));           random[4][P00-1]=Frac[w][12];          } 159            //Fracture goes through the top of the cell          else          {           fraclen=fraclen+sqrt(pow(Frac[w][3]-((cell[x][z][14]-Frac[w][10])/Frac[w][9]),2)+pow(Frac[w][4]-cell[x][z][14],2));           random[0][P00-1]=Frac[w][11];           random[1][P00-1]=sqrt(pow(Frac[w][3]-((cell[x][z][14]-Frac[w][10])/Frac[w][9]),2)+pow(Frac[w][4]-cell[x][z][14],2));           random[4][P00-1]=Frac[w][12];          }         }        }         //Check if fracture goes through the sides of the cell        else if (((((Frac[w][9]*cell[x][z][11]+Frac[w][10])>=cell[x][z][12]-0.0001) &&((Frac[w][9]*cell[x][z][11]+Frac[w][10])<=cell[x][z][14]+0.0001)) &&((((Frac[w][9]*cell[x][z][11]+Frac[w][10])>=Frac[w][2]-0.0001) &&((Frac[w][9]*cell[x][z][11]+Frac[w][10])<=Frac[w][4]+0.0001)) ||(((Frac[w][9]*cell[x][z][11]+Frac[w][10])<=Frac[w][2]+0.0001) &&((Frac[w][9]*cell[x][z][11]+Frac[w][10])>=Frac[w][4]-0.0001)))) ||((((Frac[w][9]*cell[x][z][13]+Frac[w][10])>=cell[x][z][12]-0.0001) &&((Frac[w][9]*cell[x][z][13]+Frac[w][10])<=cell[x][z][14]+0.0001)) &&((((Frac[w][9]*cell[x][z][13]+Frac[w][10])>=Frac[w][2]-0.0001) &&((Frac[w][9]*cell[x][z][13]+Frac[w][10])<=Frac[w][4]+0.0001)) ||(((Frac[w][9]*cell[x][z][13]+Frac[w][10])<=Frac[w][2]+0.0001) &&((Frac[w][9]*cell[x][z][13]+Frac[w][10])>=Frac[w][4]-0.0001)))))        {                          //if fracture goes through both sides of the cell         if (((Frac[w][9]*cell[x][z][11]+Frac[w][10])>=cell[x][z][12]-0.0001 && (Frac[w][9]*cell[x][z][11]+Frac[w][10])<=cell[x][z][14]+0.0001) && ((Frac[w][9]*cell[x][z][13]+Frac[w][10])>=cell[x][z][12]-0.0001 && (Frac[w][9]*cell[x][z][13]+Frac[w][10])<=cell[x][z][14]+0.0001) &&((Frac[w][1]<=cell[x][z][11]+0.0001 && Frac[w][3]>=cell[x][z][13]-0.0001)||(Frac[w][3]<=cell[x][z][11]+0.0001 && Frac[w][1]>=cell[x][z][13]-0.0001)))         {          P00=P00+1;          fraclen=fraclen+sqrt(pow(((Frac[w][9]*cell[x][z][11]+Frac[w][10])-(Frac[w][9]*cell[x][z][13]+Frac[w][10])),2)+pow(cell[x][z][13]-cell[x][z][11],2));          random[0][P00-1]=Frac[w][11];          random[1][P00-1]=sqrt(pow(((Frac[w][9]*cell[x][z][11]+Frac[w][10])-160  (Frac[w][9]*cell[x][z][13]+Frac[w][10])),2)+pow(cell[x][z][13]-cell[x][z][11],2));          random[4][P00-1]=Frac[w][12];          }          //if fracture goes only through the xmin side of the cell         else if (((Frac[w][9]*cell[x][z][11]+Frac[w][10])>=cell[x][z][12]-0.0001 && (Frac[w][9]*cell[x][z][11]+Frac[w][10])<=cell[x][z][14]+0.0001) && ((Frac[w][9]*cell[x][z][13]+Frac[w][10])<=cell[x][z][12]+0.0001 || (Frac[w][9]*cell[x][z][13]+Frac[w][10])>=cell[x][z][14]-0.0001))         {                           //if positive slope          if (Frac[w][9]>=0-0.0001 && (((cell[x][z][14]-Frac[w][10])/Frac[w][9]<=Frac[w][1]+0.0001 && (cell[x][z][14]-Frac[w][10])/Frac[w][9]>=Frac[w][3]-0.0001)||((cell[x][z][14]-Frac[w][10])/Frac[w][9]>=Frac[w][1]-0.0001 && (cell[x][z][14]-Frac[w][10])/Frac[w][9]<=Frac[w][3]+0.0001)))          {           P00=P00+1;           fraclen=fraclen+sqrt(pow((cell[x][z][14]-(Frac[w][9]*cell[x][z][11]+Frac[w][10])),2)+pow((1/Frac[w][9])*(cell[x][z][14]-(Frac[w][9]*cell[x][z][11]+Frac[w][10])),2));           random[0][P00-1]=Frac[w][11];           random[1][P00-1]=sqrt(pow((cell[x][z][14]-(Frac[w][9]*cell[x][z][11]+Frac[w][10])),2)+pow((1/Frac[w][9])*(cell[x][z][14]-(Frac[w][9]*cell[x][z][11]+Frac[w][10])),2));           random[4][P00-1]=Frac[w][12];          }          //if negative slope          else if (Frac[w][9]<0 && (((cell[x][z][12]-Frac[w][10])/Frac[w][9]<=Frac[w][1]+0.0001 && (cell[x][z][12]-Frac[w][10])/Frac[w][9]>=Frac[w][3]-0.0001)||((cell[x][z][12]-Frac[w][10])/Frac[w][9]>=Frac[w][1]-0.0001 && (cell[x][z][12]-Frac[w][10])/Frac[w][9]<=Frac[w][3]+0.0001)))          {           P00=P00+1;           fraclen=fraclen+sqrt(pow(((Frac[w][9]*cell[x][z][11]+Frac[w][10])-cell[x][z][12]),2)+pow((1/Frac[w][9])*((Frac[w][9]*cell[x][z][11]+Frac[w][10])-cell[x][z][12]),2));           random[0][P00-1]=Frac[w][11];           random[1][P00-1]=sqrt(pow(((Frac[w][9]*cell[x][z][11]+Frac[w][10])-cell[x][z][12]),2)+pow((1/Frac[w][9])*((Frac[w][9]*cell[x][z][11]+Frac[w][10])-cell[x][z][12]),2)); 161            random[4][P00-1]=Frac[w][12];          }          else          {}                  }          //if fracture goes only through the xmax side of the cell         else if (((Frac[w][9]*cell[x][z][13]+Frac[w][10])>=cell[x][z][12]-0.0001 && (Frac[w][9]*cell[x][z][13]+Frac[w][10])<=cell[x][z][14]+0.0001) && ((Frac[w][9]*cell[x][z][11]+Frac[w][10])<=cell[x][z][12]+0.0001 || (Frac[w][9]*cell[x][z][11]+Frac[w][10])>=cell[x][z][14]-0.0001))         {          //if positive slope          if (Frac[w][9]>=0 && (((cell[x][z][12]-Frac[w][10])/Frac[w][9]<=Frac[w][1]+0.0001 && (cell[x][z][12]-Frac[w][10])/Frac[w][9]>=Frac[w][3]-0.0001)||((cell[x][z][12]-Frac[w][10])/Frac[w][9]>=Frac[w][1]-0.0001 && (cell[x][z][12]-Frac[w][10])/Frac[w][9]<=Frac[w][3]+0.0001)))          {           P00=P00+1;           fraclen=fraclen+sqrt(pow((cell[x][z][12]-(Frac[w][9]*cell[x][z][13]+Frac[w][10])),2)+pow((1/Frac[w][9])*(cell[x][z][12]-(Frac[w][9]*cell[x][z][13]+Frac[w][10])),2));           random[0][P00-1]=Frac[w][11];           random[1][P00-1]=sqrt(pow((cell[x][z][12]-(Frac[w][9]*cell[x][z][13]+Frac[w][10])),2)+pow((1/Frac[w][9])*(cell[x][z][12]-(Frac[w][9]*cell[x][z][13]+Frac[w][10])),2));           random[4][P00-1]=Frac[w][12];          }          //if negative slope          else if (Frac[w][9]<0 && (((((cell[x][z][14]-Frac[w][10])/Frac[w][9])<=Frac[w][1]+0.0001) && (cell[x][z][14]-Frac[w][10])/Frac[w][9]>=Frac[w][3]-0.0001)||((((cell[x][z][14]-Frac[w][10])/Frac[w][9])>=Frac[w][1]-0.0001) && (cell[x][z][14]-Frac[w][10])/Frac[w][9]<=Frac[w][3]+0.0001)))          {           P00=P00+1;           fraclen=fraclen+sqrt(pow(((Frac[w][9]*cell[x][z][13]+Frac[w][10])-cell[x][z][14]),2)+pow((1/Frac[w][9])*((Frac[w][9]*cell[x][z][13]+Frac[w][10])-cell[x][z][14]),2));           random[0][P00-1]=Frac[w][11];           random[1][P00-1]=sqrt(pow(((Frac[w][9]*cell[x][z][13]+Frac[w][10])-cell[x][z][14]),2)+pow((1/Frac[w][9])*((Frac[w][9]*cell[x][z][13]+Frac[w][10])-cell[x][z][14]),2));           random[4][P00-1]=Frac[w][12]; 162           }          else          {}          }         }         // check if the fracture goes through the top and bottom of the cell        else if((((cell[x][z][12]-Frac[w][10])/Frac[w][9])>=cell[x][z][11]-0.0001 && ((cell[x][z][12]-Frac[w][10])/Frac[w][9])<=cell[x][z][13]+0.0001) && (((cell[x][z][14]-Frac[w][10])/Frac[w][9])>=cell[x][z][11]-0.0001 && ((cell[x][z][12]-Frac[w][10])/Frac[w][9])<=cell[x][z][13]+0.0001)&& ((((cell[x][z][12]-Frac[w][10])/Frac[w][9])<=Frac[w][1]+0.0001 && ((cell[x][z][12]-Frac[w][10])/Frac[w][9])>=Frac[w][3]-0.0001)|| (((cell[x][z][12]-Frac[w][10])/Frac[w][9])>=Frac[w][1]-0.0001 && ((cell[x][z][12]-Frac[w][10])/Frac[w][9])<=Frac[w][3]+0.0001))&&((((cell[x][z][14]-Frac[w][10])/Frac[w][9])<=Frac[w][1]+0.0001 && ((cell[x][z][14]-Frac[w][10])/Frac[w][9])>=Frac[w][3]-0.0001)|| (((cell[x][z][14]-Frac[w][10])/Frac[w][9])>=Frac[w][1]-0.0001 && ((cell[x][z][14]-Frac[w][10])/Frac[w][9])<=Frac[w][3]+0.0001)))        {         P00=P00+1;         fraclen=fraclen+sqrt(pow((((cell[x][z][12]-Frac[w][10])/Frac[w][9])-((cell[x][z][14]-Frac[w][10])/Frac[w][9])),2)+pow(cell[x][z][14]-cell[x][z][12],2));         random[0][P00-1]=Frac[w][11];         random[1][P00-1]=sqrt(pow((((cell[x][z][12]-Frac[w][10])/Frac[w][9])-((cell[x][z][14]-Frac[w][10])/Frac[w][9])),2)+pow(cell[x][z][14]-cell[x][z][12],2));         random[4][P00-1]=Frac[w][12];        }         else        {}         w=w+1;       }       //P21       cell[x][z][15]=fraclen/(SizeCellx*SizeCelly);       //P00       cell[x][z][16]=P00;       //P20       cell[x][z][17]=P00/(SizeCellx*SizeCelly);       //P11*       cell[x][z][18]=fraclen/((SizeCellx+SizeCelly)*2);       //P01       cell[x][z][19]=fraclen;        //RandDIP and RandLength       if (P00!=0)       { 163         V=((rand()%P00));        cell[x][z][20]=random[0][V];        cell[x][z][21]=random[1][V];       }       else       {        cell[x][z][20]=0;        cell[x][z][21]=0;       }        //For the calculation of the resultant Frac and the nb of similar oriented fracture in the cell       //if there is a fracture in the cell and the size of the frac is not null       if (P00!=0 && cell[x][z][19]!=0)       {        s=0;        resulx=0;        resuly=0;        Big=0;        cell[x][z][26]=random[0][0];        cell[x][z][27]=random[1][0];        cell [x][z][28]=random[0][0];        if (random[0][0]<35)        {         criti=(random[4][0])/(-0.0274*fabs(random[0][0])+5.5544);        }        else        {         criti=(random[4][0])/(0.00478225*pow(fabs(random[0][0]),2)-0.58035*fabs(random[0][0])+19.0497);        }        while (s<P00)        {         //Calculations for the resultant fracture         resulx=resulx+sqrt(pow(random[1][s]*cos(random[0][s]*(PI/180)),2));         resuly=resuly+(random[1][s]*sin(random[0][s]*(PI/180)));          //Calculations for the most represented fracture         q=0;         random[2][s]=0;         random[3][s]=0;         while (q<P00)         {                    //if angle comprise between random-15 et random+15          if ((random[0][q]>=(random[0][s]-15)) && (random[0][q]<=(random[0][s]+15)))          { 164            random[2][s]=random[2][s]+1;           random[3][s]=random[3][s]+random[1][q];          }           else          {}          q=q+1;         }          //if number option chosen for most represented         if (Most==1)         {           //Finds the most represented (number) fracture or most critical if equally represented          if (random[2][Big]<=random[2][s])          {           //finds most critical if equally represented           if ((random[2][Big]==random[2][s])&&(Big!=s)&&(random[0][Big]!=random[0][s]))           {            //if DIP Big more critical than DIP s            if (fabs(45-fabs(random[0][Big]))<fabs(45-fabs(random[0][s])))            {             Big=Big;            }             //if DIP Big less critical than DIP s            else            {             Big=s;            }           }           else           {            Big=s;           }          }          else          {}         }          else         {                    //Finds the most represented (length) fracture or most critical if equally represented 165           if (random[3][Big]<=random[3][s])          {           //finds most critical if equally represented           if ((random[3][Big]==random[3][s])&&(Big!=s)&&(random[0][Big]!=random[0][s]))           {            //if DIP Big more critical than DIP s            if (fabs(45-fabs(random[0][Big]))<fabs(45-fabs(random[0][s])))            {             Big=Big;            }             //if DIP Big less critical than DIP s            else            {             Big=s;            }           }           else           {            Big=s;           }          }          else          {}                   }          //to find the most critical fracture (fracture must have a certain length to be considered)         if (fabs(60-fabs(cell[x][z][26]))>fabs(60-fabs(random[0][s]))&& (random[1][s]>=(SizeCelly/10)))         {          cell[x][z][26]=random[0][s];          cell[x][z][27]=random[1][s];         }         else         {         }          //to find most critical fracture with length component          if (random[0][s]<35)         {          if ((criti)>((random[4][s])/(-0.0274*fabs(random[0][s])+5.5544))) 166           {          }           else          {            criti=((random[4][s])/(-0.0274*fabs(random[0][s])+5.5544));                     cell[x][z][28]=random[0][s];          }                   }         else         {          if ((criti)>((random[4][s])/(0.00478225*pow(fabs(random[0][s]),2)-0.58035*fabs(random[0][s])+19.0497)))          {          }           else          {            criti=((random[4][s])/(0.00478225*pow(fabs(random[0][s]),2)-0.58035*fabs(random[0][s])+19.0497));                     cell[x][z][28]=random[0][s];          }          }          s=s+1;        }        cell[x][z][22]=(atan(resuly/resulx))*180/PI;        cell[x][z][23]=sqrt(pow(resulx,2)+pow(resuly,2));        cell[x][z][24]=random[0][Big];        cell[x][z][25]=random[3][Big];       }        //if no fracture in the cell or fracture length is null       else       {        cell[x][z][22]=0;        cell[x][z][23]=0;       }        //cell Rock Properties        //if use GSI       if ((QGSI=='Y' || QGSI=='y') && (cell[x][z][16]!=0))       { 167         if(((sqrt(285/(SizeCellx*SizeCelly))*cell[x][z][19])/285)>=0.09)        {         GSI=-10.29*log((sqrt(285/(SizeCellx*SizeCelly))*cell[x][z][19])/285)+68.86;        }         else        {         GSI=-77.778*((sqrt(285/(SizeCellx*SizeCelly))*cell[x][z][19])/285)+100;        }                S=exp((GSI-100)/(9-(3*0)));//Hoek Brown 2002        mb=12*exp((GSI-100)/(28-(14*0)));//Hoek Brown 2002        cell[x][z][0]=Young*(0.02+(1-(0/2))/(1+exp((60+15*0-GSI)/11)));//Hoek and Diederichs 2006 (Rock Mass Young Modulus)         cell[x][z][1]=Poisson;        cell[x][z][2]=Cohesion;        cell[x][z][3]=((S*UCS)/mb);//Hoek Brown 2002 (Rock Mass Tension)        cell[x][z][4]=Friction;        cell[x][z][5]=Dilation;        cell[x][z][6]=Jcohesion;        cell[x][z][7]=Jfriction;        cell[x][z][8]=Jdilation;        cell[x][z][9]=Jtension;       }               //if don't use GSI or no joint in cell       else       {        cell[x][z][0]=Young;        cell[x][z][1]=Poisson;        cell[x][z][2]=Cohesion;        cell[x][z][3]=Tension;        cell[x][z][4]=Friction;        cell[x][z][5]=Dilation;        cell[x][z][6]=Jcohesion;        cell[x][z][7]=Jfriction;        cell[x][z][8]=Jdilation;        cell[x][z][9]=Jtension;       }              //Joint Angle       //If Ubiquitous joint model and random frac chosen       if ((Ubi=='Y' || Ubi=='y') && Choice==1)       {        //if angle positive        if (cell[x][z][20]>0)        { 168          cell[x][z][10]=cell[x][z][20];        }        //if angle negative        else        {         cell[x][z][10]=180+cell[x][z][20];        }       }           //If Ubiquitous joint model and most represented chosen       else if ((Ubi=='Y' || Ubi=='y') && Choice==3)       {        //if angle positive        if (cell[x][z][24]>0)        {         cell[x][z][10]=cell[x][z][24];        }        //if angle negative        else        {         cell[x][z][10]=180+cell[x][z][24];        }       }        //If Ubiquitous joint model and Critical chosen       else if ((Ubi=='Y' || Ubi=='y') && Choice==4)       {         //if angle positive        if (cell[x][z][26]>0)        {         cell[x][z][10]=cell[x][z][26];        }         //if angle negative        else        {         cell[x][z][10]=180+cell[x][z][26];        }         }              //If Ubiquitous joint model and Overall Critical chosen       else if ((Ubi=='Y' || Ubi=='y') && Choice==5)       {         //if angle positive 169         if (cell[x][z][28]>0)        {         cell[x][z][10]=cell[x][z][28];        }         //if angle negative        else        {         cell[x][z][10]=180+cell[x][z][28];        }         }       //if Resultant Frac       else       {        //if angle positive        if (cell[x][z][22]>0)        {         cell[x][z][10]=cell[x][z][22];        }        //if angle negative        else        {         cell[x][z][10]=180+cell[x][z][22];        }       }        z=z+1;      }      x=x+1;     }        //Generates the grid     out<<"new"<<endl;     out<<"grid "<<ColCell<<" "<<RowCell<<endl;      //If Ubiquitous joint model     if (Ubi=='Y' || Ubi=='y')     {      x=0;      t=1;      while (x<RowCell)      {       z=0;       while (z<ColCell)       {         //if no fractures in cell        if (cell[x][z][19]==0)        { 170          out<<"group 'User:new"<<t<<"' "<<"i "<<(z+1)<<" j "<<(x+1)<<endl;         out<<"model mohr group 'User:new"<<t<<"'"<<endl;         out<<"prop b "<<((cell[x][z][0]/3)/(1-2*cell[x][z][1]))<<" c "<<cell[x][z][2]<<" d "<<Density<<" di "<<cell[x][z][5]<<" f "<<cell[x][z][4]<<" s "<<((cell[x][z][0]/2)/(1+cell[x][z][1]))<<" t "<<cell[x][z][3]<<" group 'User:new"<<t<<"'"<<endl;         }         //if fracture in cell        else        {         out<<"group 'User:new"<<t<<"' "<<"i "<<(z+1)<<" j "<<(x+1)<<endl;         out<<"model ubiquitous group 'User:new"<<t<<"'"<<endl;         out<<"prop b "<<((cell[x][z][0]/3)/(1-2*cell[x][z][1]))<<" c "<<cell[x][z][2]<<" d "<<Density<<" di "<<cell[x][z][5]<<" f "<<cell[x][z][4]<<" s "<<((cell[x][z][0]/2)/(1+cell[x][z][1]))<<" t "<<cell[x][z][3]<<" ja "<<cell[x][z][10]<<" jc "<<cell[x][z][6]<<" jd "<<cell[x][z][8]<<" jf "<<cell[x][z][7]<<" jt "<<cell[x][z][9]<<" group 'User:new"<<t<<"'"<<endl;        }                        z=z+1;        t=t+1;       }       x=x+1;      }           }          //Mohr model     else     {      x=0;      t=2;      while (x<RowCell)      {       z=0;       while (z<ColCell)       {        //if no fractures in cell        if (cell[x][z][19]==0)        {         out<<"group 'User:new"<<t<<"' "<<"i "<<(z+1)<<" j "<<(x+1)<<endl;         out<<"model mohr group 'User:new"<<t<<"'"<<endl;         out<<"prop b "<<((cell[x][z][0]/3)/(1-2*cell[x][z][1]))<<" c "<<cell[x][z][2]<<" d "<<Density<<" di "<<cell[x][z][5]<<" f "<<cell[x][z][4]<<" s 171  "<<((cell[x][z][0]/2)/(1+cell[x][z][1]))<<" t "<<cell[x][z][3]<<" group 'User:new"<<t<<"'"<<endl;        }         //if fracture in cell        else        {         out<<"group 'User:new"<<t<<"' "<<"i "<<(z+1)<<" j "<<(x+1)<<endl;         out<<"model mohr group 'User:new"<<t<<"'"<<endl;         out<<"prop b "<<BulkF<<" c "<<CohesionF<<" d "<<DensityF<<" di "<<DilationF<<" f "<<FrictionF<<" s "<<ShearF<<" t "<<TensionF<<" group 'User:new"<<t<<"'"<<endl;        }        z=z+1;        t=t+1;       }       x=x+1;      }      }      //Sets the sample height and width     out<<"gen 0,0 0,"<<height<<" "<<width<<","<<height<<" "<<width<<",0 i=1,"<<ColCell+1<<" j=1,"<<RowCell+1<<endl;      //To calculate Strain     out<<"def ve"<<endl;     out<<" ve=(ydisp("<<ColCell/2+1<<",1)-ydisp("<<ColCell/2+1<<","<<RowCell+1<<"))/(y("<<ColCell/2+1<<","<<RowCell+1<<")-y("<<ColCell/2+1<<",1))"<<endl;     out<<"end"<<endl;      //To calculate Sigma Vertical     out<<"def sigmav"<<endl;     out<<" sum=0.0"<<endl;     out<<" loop i (1,igp)"<<endl;     out<<"  sum=sum+yforce(i,jgp)"<<endl;     out<<" end_loop"<<endl;     out<<" sigmav=sum/(x(igp,jgp)-x(1,jgp))"<<endl;     out<<"end"<<endl;          //averaging major and minor principal stress in pillar     out<<"def pillar1"<<endl;     out<<" sum1=0.0"<<endl;     out<<" sum3=0.0"<<endl;     out<<" loop i (1,"<<ColCell<<")"<<endl;     out<<"  loop j (1,"<<RowCell<<")"<<endl;     out<<"   temp1=-0.5*(sxx(i,j)+syy(i,j))"<<endl;     out<<"   temp2=sqrt(sxy(i,j)^2+0.25*(sxx(i,j)-syy(i,j))^2)"<<endl;     out<<"   s1=max(temp1+temp2,-szz(i,j))"<<endl; 172      out<<"   s3=min(temp1-temp2,-szz(i,j))"<<endl;     out<<"   sum1=sum1+s1"<<endl;     out<<"   sum3=sum3+s3"<<endl;     out<<"  end_loop"<<endl;     out<<" end_loop"<<endl;     out<<" pillar1=sum1/"<<ColCell*RowCell<<endl;     out<<" pillar3=sum3/"<<ColCell*RowCell<<endl;     out<<"end"<<endl;      //sets the Damping     //If Combined Damping mode chosen     if (Damp=='Y' || Damp=='y')     {      out<<"set st_damp comb"<<endl;     }      else     {}      //sets the gravity     //If gravity mode chosen     if (Grav=='Y' || Grav=='y')     {      out<<"set grav=9.81"<<endl;      out<<"solve"<<endl;     }      else     {}      out<<"save sample"<<OutNameI<<".sav"<<endl;     out<<"call Hydrostat"<<OutNameI<<".dat"<<endl;      strcpy_s(OutNameh,"Hydrostat");     strcat_s(OutNameh,OutName);      outh.open(OutNameh);      if (outh.fail())     {      cout<<"opening problem with: "<<OutNameh<<endl<<endl;      }      outh<<"new"<<endl;     outh<<"restore sample"<<OutNameI<<".sav"<<endl;      //Apply hydrostat pressure     outh<<"ap pr "<<conf<<" j=1"<<endl;     outh<<"ap pr "<<conf<<" j="<<RowCell+1<<endl;     outh<<"ap pr "<<conf<<" i=1"<<endl;     outh<<"ap pr "<<conf<<" i="<<ColCell+1<<endl;     outh<<"solve"<<endl;      outh<<"save Hydrostat"<<OutNameI<<".sav"<<endl; 173      outh<<"call Servo"<<OutNameI<<".dat"<<endl;      strcpy_s(OutNameS,"Servo");     strcat_s(OutNameS,OutName);      outs.open(OutNameS);      if (outs.fail())     {      cout<<"opening problem with: "<<OutNameS<<endl<<endl;      }      outs<<"new"<<endl;     outs<<"restore hydrostat"<<OutNameI<<".sav"<<endl;      //Apply initial velocity     outs<<"ap yvel="<<(-ini)<<" j="<<(RowCell+1)<<endl;     outs<<"ap yvel="<<ini<<" j=1"<<endl;      //Servocontrol to apply velocity     outs<<"def servo"<<endl;     outs<<"while_stepping"<<endl;     outs<<" if unbal > 1e6 then"<<endl;     outs<<"  loop i (1,"<<(ColCell+1)<<")"<<endl;     outs<<"   yvel(i,"<<RowCell+1<<")=yvel(i,"<<RowCell+1<<")*.975"<<endl;     outs<<"   yvel(i,1)=yvel(i,1)*.975"<<endl;     outs<<"  end_loop"<<endl;     outs<<" end_if"<<endl;     outs<<" if unbal > 1e5 then"<<endl;     outs<<"  loop i (1,"<<(ColCell+1)<<")"<<endl;     outs<<"   yvel(i,"<<RowCell+1<<")=yvel(i,"<<RowCell+1<<")*1.025"<<endl;     outs<<"   yvel(i,1)=yvel(i,1)*1.025"<<endl;     outs<<"  end_loop"<<endl;     outs<<" end_if"<<endl;     outs<<"end"<<endl;      outs<<"history sigmav"<<endl;     outs<<"history ve"<<endl;     outs<<"history pillar1"<<endl;     outs<<"history pillar3"<<endl;       // deletes memory space used by cell[][][]     x=0;     while (x<RowCell)     {            z=0;      while (z<ColCell) 174       {       delete []cell[x][z];       z=z+1;      }      delete []cell[x];      x=x+1;     }     delete []cell;           }    outh.close();    outs.close();    out.close();      }   in2.close();  }  in.close(); }      Appendix D: Slope Models   Figure D1-1 Slope2 1 m mesh x-displacement contours  175   Figure D1-2 Slope2 1 m mesh xx-stress contours   Figure D1-3 Slope2 1 m mesh yy-stress contours  176   Figure D1-4 Slope2 2 m mesh x-displacement contours   Figure D1-5 Slope2 2 m mesh xx-stress contours  177   Figure D1-6 Slope2 2 m mesh yy-stress contours    

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0052524/manifest

Comment

Related Items