You may notice some images loading slow across the Open Collections website. Thank you for your patience as we rebuild the cache to make images load faster.

Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Derivation of an equivalent boundary method for ground-support interaction problems Mitelman, Amichai 2020

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.

Item Metadata

Download

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

Full Text

i  DERIVATION OF AN EQUIVALENT BOUNDARY METHOD FOR GROUND-SUPPORT INTERACTION PROBLEMS by Amichai Mitelman    A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE AND POSTDOCTORAL STUDIES (Mining Engineering)  THE UNIVERISTY OF BRITISH COLUMBIA (Vancouver)   April 2020  © Amichai Mitelman, 2020     ii  The following individuals certify that they have read, and recommend to the Faculty of Graduate and Postdoctoral Studies for acceptance, the dissertation entitled:  Derivation of an equivalent boundary method for ground-support interaction problems                 submitted by    Amichai   Mitelman                          in partial fulfillment of the requirements for the degree of    Doctor of Philosophy  in               Mining Engineering  Examining Committee: Davide Elmo, Professor, Mining Engineering, UBC Supervisor  Scott Dunbar, Mining Engineering, UBC Supervisory Committee Member  Erik Eberhardt, Professor, Earth, Ocean and Atmospheric Sciences, UBC Supervisory Committee Member Marek Pawel, Professor, Mining Engineering, UBC University Examiner Carlos Ventura, Professor, Civil Engineering, UBC  University Examiner   Additional Supervisory Committee Members: Doug Stead, Professor,  Earth Sciences, SFU Supervisory Committee Member Yossef Hatzor, Professor, Geological and Environmental Sciences, BGU External Examiner iii  Abstract This thesis presents a novel approach termed ‘Equivalent Boundary’ (EB) for the analysis of ground-support interaction problems. The basic idea of the proposed method is to simplify these problems by representing the ground as analogous structural entities. Similar to the convergence-confinement method, increased efficiency is attained by focusing on the ground-support boundary, rather than simulating a great portion of the surrounding ground, as is required in finite-element models. Within this thesis, a number of different ground-support problems are addressed: room and pillar mines, vertical circular shafts, and tunnels. Different structural analogues are chosen for each problem according to the nature of each problem: the pillar is represented by a spring, the shaft by a ring, and the tunnel by a series of beam elements. Expressions for the stiffness of the structural analogues are derived. Subsequently, the unsupported ground displacements are used as an input; based on these, the displacements, internal forces, and factor of safety of the supported ground can readily be computed. All results have been validated against numerical models. For the tunnel problem, a methodology for the analysis of circular tunnels in elastic ground is derived. The underlying assumption of the traditional convergence-confinement approach is that the tunnel is subjected to a hydrostatic stress field, a simplification which poses considerable practical limitations. Within the proposed method, the case of a tunnel subjected to a non-uniform stress field can be addressed. Subsequently, the methodology is further modified to address two more complex conditions: 1) a circular tunnel in plastic ground, and 2) a non-circular horseshoe shaped tunnel in elastic ground.  Due to the efficient computational process, the method developed in this thesis is well-suited for probabilistic analyses that require a large number of iterations. A methodology for the cost estimation of tunnel support based on the construction method is presented. A practical example is used to demonstrate the advantages of the EB method developed in this thesis. Additionally, this methodology constitutes a useful stand-alone concept, which can be implemented using other available tunnel analysis methods.   iv  Lay Summary Ground-support problems are a type of problems in engineering where a support system interacts with the surrounding ground. Projects that include ground-support interaction are dealt with by practicing civil, mining and geotechnical engineers. Tunnels, vertical shafts, and room and pillar mines are all examples of such problems, all addressed in this thesis. For the design of ground support there are a number of methods available, including analytical, empirical, and numerical methods. In this thesis, a novel approach termed ‘Equivalent Boundary’ (EB) for the analysis of ground-support problems is introduced. In this approach, the problems are simplified by converting the surrounding ground into equivalent analogous structures. In terms of simplicity and computational efficiency, this method can be viewed as an intermediate method; compared to analytical closed-form solutions, the proposed method requires less simplifying assumptions; compared to available numerical programs, the EB method requires significantly less computational resources.   v  Preface This thesis represents original and independent work. A number of papers have been published based on the research presented in this thesis. Amichai Mitleman is the first author on two peer-reviewed published journal papers and a conference paper, in addition to being a co-author on a conference paper on work related to the thesis.   The journal paper “Development of a spring analogue approach for the study of pillars and shafts” was published in the International Journal of Mining Science and Technology (2018) and is based on Chapters  4 and  5.   The journal paper “Analysis of tunnel-support interaction using an equivalent boundary beam” was published in the journal of Tunnelling and Underground Space Technology (2019) and is based on Chapter  6.  The conference paper “A proposed probabilistic analysis methodology for tunnel support cost estimation depending on the construction method” was published in the Proceedings of the 52nd U.S. Rock Mechanics/Geomechanics Symposium (2018) and is based on Chapter  7.   The conference paper “Design and construction of a dense array of small diameter holes in a large cavern sidewall” was published in the Proceedings of the 53rd U.S. Rock Mechanics/Geomechanics Symposium (2019) and is related to the practical example discussed in Section 6.5. In this paper the thesis supervisor was the lead author and contributed FEM/DEM simulation and analysis. The project tunnel designers, Arnon Rozen and Dan Tzuker, the tunnel engineers responsible for the design of the tunneling project, are credited as additional co-authors. During the work on this thesis, the author contributed to an addition paper titled “A methodology for the analysis of tunnel intersections using two-dimensional numerical modeling”, which has recently been submitted as a conference paper. This paper summarizes the Master’s thesis by Sahar Cohen from Graz University of Technology, who is the lead author. The author of this thesis helped develop the idea and work for the paper, and is co-author along with Dr. Elmo.   vi   Table of Contents Abstract ...................................................................................................................................... ii Lay Summary ............................................................................................................................ iv Preface ........................................................................................................................................ v Table of Contents ...................................................................................................................... vi List of Tables........................................................................................................................... viii List of Figures ........................................................................................................................... ix List of Abbreviations................................................................................................................ xii Acknowledgements ................................................................................................................. xiii Dedication ............................................................................................................................... xiv 1. CHAPTER 1  INTRODUCTION ....................................................................................... 1 1.1. Problem Statement ...................................................................................................... 1 1.2. Research Objectives .................................................................................................... 3 1.3. Thesis Organization ..................................................................................................... 4 2. CHAPTER 2  LITERATURE REVIEW ............................................................................ 6 2.1. Analytical and Analogue Methods .............................................................................. 7 2.2. Numerical Methods ................................................................................................... 22 2.3. Relationship Between Different Analysis Methods .................................................. 23 3. CHAPTER 3 METHODOLOGY ..................................................................................... 34 3.1. General Methodology ................................................................................................ 34 3.2. Basic Formulation ..................................................................................................... 35 4. CHAPTER 4 APPLICATION OF THE EB METHOD FOR PILLAR PROBLEMS ..... 43 4.1. Introduction ............................................................................................................... 43 4.2. Derivation .................................................................................................................. 44 4.3. Comparison with Empirical Pillar Formulae ............................................................. 51 4.4. Influence of Support on the Idealized Mine-Pillar System ....................................... 54 4.5. Discussion of Post-Peak Pillar Analysis ................................................................... 58 4.5.1. Dynamic Phase ...................................................................................................... 59    4.5.2. Residual State ..................................................................................................... 60 vii  5. CHAPTER 5  APPLICATION OF THE EB METHOD FOR VERTICAL                                                                                        SHAFT PROBLEMS ............................................................................................................... 63 5.1. Introduction ............................................................................................................... 63 5.2. Derivation and Results .............................................................................................. 65 6. CHAPTER 6  APPLICATION OF THE EB METHOD FOR TUNNEL PROBLEMS .. 70 6.1. Introduction ............................................................................................................... 70 6.2. Derivation for Elastic Ground ................................................................................... 72 6.3. Validation for Tunnel Liner Analysis ....................................................................... 94 6.4. Rock Bolt Analysis .................................................................................................... 99 6.4.1. Pre-tensioned Bolts .......................................................................................... 100 6.4.2. Mechanically Anchored Bolts .......................................................................... 101 6.4.3. Grouted Bolts ................................................................................................... 102 6.5. Derivation for Plastic Ground ................................................................................. 108 6.6. Horseshoe Shaped Tunnel ....................................................................................... 114 6.6.1. Results and Discussion ..................................................................................... 122 7. CHAPTER 7  METHODOLOGY FOR TUNNEL SUPPORT PROBABILISTIC ANALYSIS ............................................................................................................................ 126 7.1. Introduction ............................................................................................................. 126 7.2. Tunneling Method Influence on Probabilistic Analysis .......................................... 129 7.3. Brittle Failure Criterion ........................................................................................... 135 7.4. Practical Tunnel Analysis Example......................................................................... 141 8. CHAPTER 8  CONCLUSIONS AND RECOMMENDATIONS .................................. 150 8.1. Summary of key contributions ................................................................................ 150 8.2. Recommendation for Further Studies ...................................................................... 152 References .............................................................................................................................. 156 Appendix 1 ............................................................................................................................. 165    viii  List of Tables Table 1 Structural analogues for different ground-support problems ...................................... 36 Table 2 Starting point of known's and unknowns in typical structural problems compared to the EB method .......................................................................................................................... 40 Table 3 Properties for mine-pillar parametric study ................................................................ 49 Table 4 Comparison of pillar factor-of-safety using the Lunder and Pakalnis and the proposed spring analogue ........................................................................................................................ 53 Table 5 Properties and geometry for parametric study ............................................................ 57 Table 6 Comparison of bolting mine roof with bolting of pillar on FoS at different depths ... 57 Table 7 Properties for shaft parametric study .......................................................................... 68 Table 8 Notable analytical solutions for circular and deep tunnels and their limitations ........ 71 Table 9 Comparison of results of maximum major principal stresses for lined tunnels .......... 91 Table 10 Result comparison for lined tunnels. E&S refers to results obtained via the analytical solution given by Einstein and Schwartz (1979)...................................................................... 96 Table 11 Parameters for model comparison ............................................................................ 97 Table 12 Result comparison for bolted tunnels...................................................................... 102 Table 13 Liner and bolt parameters used for combined liner-bolt support comparison ........ 102 Table 14 Result comparison for bolted tunnels...................................................................... 107 Table 15 Liner and bolt parameters used for combined liner-bolt support comparison ........ 107 Table 16 Result comparison of combined liner-bolt support system ..................................... 107 Table 17 Tunnel parameters ................................................................................................... 111 Table 18 Tunnel crown and wall displacements per stage ..................................................... 111 Table 19 Results for liner axial forces, bending moments, and displacements, obtained from numerical models in RS2 and the proposed EBB method ..................................................... 111 Table 20 Rock input parameters ............................................................................................ 118 Table 21 Comparison of results using different methods, for maximum displacements at the tunnel crown ........................................................................................................................... 124 Table 22 Comparison of results using different methods, for maximum displacements at the tunnel wall .............................................................................................................................. 124 Table 23 Parameters and standard deviations for analysis..................................................... 143 Table 24 Support classes for SEM modeling ......................................................................... 145 Table 25 Cost estimation result comparison .......................................................................... 148 Table 26 Comparison of the required computational resources for Monte Carlo analysis of EB vs. FEM (RS2) ....................................................................................................................... 149    ix  List of Figures Figure 1 Thesis structure, main research chapters ..................................................................... 4 Figure 2 Example of a convergence-confinement plot .............................................................. 9 Figure 3 Ground-support problem analogue representation by: a) two springs in series,  and b) alternative graphical representation ......................................................................................... 10 Figure 4 Comparison of stress results from FE models for a) hydrostatic stress field vs. b) non-hydrostatic stress field. The tunnel deformation deformed shape is exaggerated for the sake of illustration and appears in gray .................................................................................... 12 Figure 5 Convergence-confinement curves (blue and purple, respectively) given by the program RocSupport, by RocScience (2015) .......................................................................... 13 Figure 6 Longitudinal displacement profile according to Chern et al. (1998). The ratio of radial displacement to the maximum longitudinal displacement is plotted against the ratio of the distance from tunnel face to tunnel radius. ........................................................................ 15 Figure 7 Example of a simplification of a truss problem using a spring analogue .................. 18 Figure 8 Foundation slab analysis using "Winkler" springs .................................................... 19 Figure 9 HRM model (modified from Oreste, 2006) ............................................................... 21 Figure 10 Evolutionary processes in a) biology and b) technology c) scientific analysis methods. Images modified from public websites ..................................................................... 25 Figure 11 Diagram of forces on a failing soil block used for Coulomb's theory for retaining wall analysis. ............................................................................................................................ 27 Figure 12 Newmark's sliding block analogy for dynamic analysis of slopes .......................... 27 Figure 13 General trend of scientific progress vs. the minimalistic approach known as “Occam's razor” ....................................................................................................................... 28 Figure 14 Viewing alternative analysis methods as complementary methods ........................ 30 Figure 15 The role of uncertainty, from data collection to complex numerical analysis methods .................................................................................................................................... 31 Figure 16 Example of the major stress contours from a 2D FE tunnel model, which can be divided into 3 different parts: 1) the ground-support interface, 2) the zone of influence, and 3) far field region. ......................................................................................................................... 39 Figure 17 Summary of the EB calculation stages. Ffic is the vector of fictitious forces, Kg represents the stiffness matrix of the equivalent ground, Ks is the stiffness matrix of the support, Dus is the unsupported displacement vector, and Ds is the displacement of the supported ground. ..................................................................................................................... 42 Figure 18  Rectangular excavation with pillar at its centre ...................................................... 45 Figure 19 Equivalent force 𝐹𝐸𝑄 .............................................................................................. 45 Figure 20 Mine-pillar system idealization to a spring model .................................................. 47 Figure 21 Proposed mine-pillar interaction curve .................................................................... 48 Figure 22 Mesh for numerical model built in RS2 for L=26 m. .............................................. 49 Figure 23 Comparison of pillar deformation computed by the tributary method, FEM analysis, and the proposed spring analogue ............................................................................................ 50 x  Figure 24 Percentage of error for displacement calculation using the tributary method vs. ratio of pillar to mine rock mass Young’s modulus ......................................................................... 51 Figure 25 Lunder and Pakalnis (1997) pillar formula and expected FoS for i) 3.5m wide, 7m high; and ii) 7m wide, 7m high pillar. ...................................................................................... 53 Figure 26 Artificial support effect on mine-pillar interaction for a) passive roof support b) active roof support and c) pillar confinement. The black 'X' denotes pillar failure. ................ 56 Figure 27 Illustration of spalling impact model ....................................................................... 60 Figure 28 SDOF Mine-pillar representation of post-peak dynamic phase. 𝐹𝑃 − 𝑅𝑒𝑠 and 𝐹𝑀 − 𝑅𝑒𝑠 denote the portions of the load imposed on the pillar and mine, respectively. ∆𝑃 − 𝑅𝑒𝑠 denotes the residual deformation of the pillar......................................................... 61 Figure 29 Shaft idealized into a double ring under uniform compression analogy ................. 64 Figure 30 Unsupported shaft axisymmetric model in RS2. Hs is the shaft height, and Ds is the shaft diameter ........................................................................................................................... 66 Figure 31 Shaft strain vs. shaft normalized depth .................................................................... 67 Figure 32 Results showing comparison of (a) liner hoop axial force, and (b) displacements, vs. liner thickness using FE model and ring analogy ............................................................... 69 Figure 33 Conversion of the plane strain problem to a double beam problem ........................ 74 Figure 34 Summary of the two main EBB computation steps. Kg represents the global stiffness matrix of the equivalent beam, Dus is the unsupported tunnel displacement vector, Ffic is the vector of fictitious forces, Ks ad Kst are the stiffness matrices of the support and supported tunnel, and Ds is the supported tunnel displacement vector. .................................. 75 Figure 35 Boundary beam elements. 𝐹𝑛 and 𝑀𝑛are the fictitious nodal forces and moments acting on: a) the equivalent boundary beam at the first stage, and b) the lined tunnel at the second stage. ............................................................................................................................ 76 Figure 36 Degrees of freedom for a beam element .................................................................. 77 Figure 37 Example of two beam elements ............................................................................... 80 Figure 38 Application of symmetry and constraints ................................................................ 81 Figure 39 Example of coordinate and element formulation for the assembly of two equivalent beam elements .......................................................................................................................... 82 Figure 40 Result convergence vs. number of boundary-beam elements ................................. 83 Figure 41 Examples of plots of displacement contours generated through MATLAB using the proposed methodology, where a) the K ratio is smaller than 1, and b) the K ratio is larger than 1. ............................................................................................................................................... 85 Figure 42 a) Example of two beam cross sections which share the same neutral axis.        b) Example of two beams with different neutral axes. ................................................................. 88 Figure 43 Examples of the superposition of axial and bending stresses, when a) the entire cross section is in compression, and b) the bottom fiber is in tension. .................................... 89 Figure 44 Summary of design process with EBB elastic analysis ........................................... 93 Figure 45 a) RS2 numerical model, showing the angle measured from horizontal (wall) to vertical (crown) b) Result comparison between RS2 and EBB of axial forces (using the xi  compression-negative sign convention)  c) Result comparison between RS2 and EBB of bending moments ..................................................................................................................... 98 Figure 46 Illustration of a) grouted bolts, and b) mechanically anchored bolts .................... 100 Figure 47 Schematic representation of the equivalent “rock-bar” and bolt system, and the displacements and fictitious forces along a single bar. .......................................................... 105 Figure 48 Comparison of axial forces along bolts computed in EBB and RS2. .................... 106 Figure 49 Computational process for plastic ground-support interaction analysis ................ 110 Figure 50 Comparison of liner maximum stresses computed using finite-element models in RS2, the proposed EBB method, and the Convergence-Confinement method        (denoted C-C in the legend) ...................................................................................................................... 113 Figure 51 Example of different two-lane road tunnel cross sections (modified from Thewes et al., 2013). ............................................................................................................................... 115 Figure 52 “Menuhot” cavern cross section. ........................................................................... 117 Figure 53 Drilling of small diameter holes in the cavern sidewalls ....................................... 118 Figure 54 Typical geological profile cross-section. ............................................................... 119 Figure 55 RS2 cavern model .................................................................................................. 121 Figure 56 Portion of tunnel contour computed with the EB method, and queried in RS2 for comparison. The boundary constraints for EB are shown as well. ........................................ 122 Figure 57 Comparison of displacement results along the tunnel boundary (as defined in Figure 56) using RS2 and EBB, for different liner thickness of: a) 150 mm, b) 250 mm, and c) 350 mm ......................................................................................................................................... 125 Figure 58 Example of probability of failure histogram ......................................................... 128 Figure 59 Flowchart diagram for: a) Tunnels with constant support (TBM tunnels) ............ 133 Figure 60 Example of convergence-confinement plot for SEM computation scheme .......... 134 Figure 61 Depth of Failure (DoF) for a tunnel with a radius r .............................................. 137 Figure 62 Increase in liner stress vs. impact factor ................................................................ 140 Figure 63 Tunnel required inner dimensions (in mm units) for SEM (solid line). The dashed line represents the circular TBM cross section of an 8 m diameter. ...................................... 142 Figure 64 Support class distribution per 1 km of tunnel for the SEM alternative ................. 146 Figure 65 Probability of Failure of support system for the TBM alternative ........................ 147 Figure 66 Flowchart for generalized application of equivalent boundaries ........................... 153   xii  List of Abbreviations  𝜈- Poisson ratio 𝑡𝑒𝑞- Equivalent thickness C-C- Convergence-confinement method DOF- Degree of Freedom EB- Equivalent Boundary EBB- Equivalent Boundary Beam (application of EB for tunnels). Erm- Rock mass Young’s modulus Ffic - Fictitious forces FE, FEM- Finite-element Method. FoS- Factor of Safety GRC- Ground Reaction Curve kij – a single stiffness matrix component, which corresponds to the force reaction in the ith degree-of-freedom, to a unit displacement at the jth degree-of-freedom. K- coefficient of lateral pressure, equal to the ratio of horizontal to vertical in-situ stresses. In other instances, depending on the context, K denotes the stiffness of an element or structure. LDP- Longitudinal Displacement Profile 𝑝 – (Vertical) ground pressure PoF- Probability of Failure u- Nodal displacement VoS- Volume of Support xiii  Acknowledgements I would like to express my sincere gratitude to my supervisor, Dr. Davide Elmo for his continuous support and guidance, and for sharing his broad knowledge in rock mechanics and vast experience in numerical modeling. Dr. Elmo is more than a guide and teacher; he has been a life-changing role model to me.  I would also like to thank my thesis committee members, Profs. Dunbar, Eberhardt and Stead, for taking the time to read my thesis and add their invaluable comments. Special thanks to tunnel engineer Arnon Rozen for his help regarding the study case of the Menuhot tunneling project (briefly reviewed in this thesis), and for sharing his important insights about tunnel engineering in general. Finally, I would like to thank my family, for being an everlasting source of meaning and inspiration in my life. All else pales in comparison.    xiv  Dedication To Shiri, the love of my life               1  1. CHAPTER 1  INTRODUCTION 1.1. Problem Statement Projects that include ground-support interaction are dealt with by practicing civil, mining and geotechnical engineers. Tunnels, vertical shafts, and room and pillar mines are all examples of such problems.  In the case of above-ground structures, which are more commonly dealt with, the design process includes a number of separate stages: first, the loads upon the structure (such as the self-weight of the structural components, people, equipment, etc.) are estimated. Subsequently, a stress and strain analysis of the loaded structure is carried out. Finally, the structural response is compared to its capacity. Based on the results, and according to engineering judgment and/or different local building code requirements, the structure may require modifications in order to increase its capacity. In contrast, ground-support problems cannot be addressed using a similar design approach, i.e. by assessing the ground load and then conducting a separate structural analysis of the loaded support system. This is because interaction problems consist of a more complex cause-and-effect process. In this process, the gravitational loads are distributed between both components, ground and support, depending on their relative structural stiffness. One consequence of this interaction process is that the stiffness of the support affects the magnitude of load shed from the ground to the support.  For the design of ground support there are a number of methods available, including analytical, empirical, and numerical methods. There are no universally accepted guidelines, and design practice tends to vary greatly. Hence, designers must apply engineering judgment in order to determine if their design is acceptable (Hoek et al., 2008). Additionally, there are often no clear boundaries between the roles of the different engineers working on geotechnical projects. This may lead to erroneous analysis, due to lack of fundamental understanding of the ground-support problem. For example, a civil engineer focusing on the design of a support element may wish to increase the stiffness of the element without accounting for the additional load imposed by the ground due to the effect of this 2  stiffening. Such design errors could lead to a negative impact on project performance (i.e. cost and schedule), and in the more severe case, to catastrophic outcomes. Finite-element (FE, or FEM) codes are powerful tools for structural analysis, and different problems consisting of complex geometries and material failure models can be computed and studied. However, there is still a need to develop solving methods which are more efficient in terms of user expertise, time consumption, and computational resources that they require. The practical applications where these matters have a significant impact from the user’s standpoint are discussed in this thesis.   Typically, ground-support problems are solved for the purpose of attempting to predict whether failure will occur, and what type of measurements are required for assuring an acceptable Factor of Safety (FoS). For this matter, in some problems it is only necessary to examine the domain of maximum stresses and/or displacements. For example, in a room and pillar problem, the maximum stresses that develop within the pillar dictate whether failure will occur, whereas stresses that develop in far field regions of the “room” do not impact the overall stability. Nevertheless, in order to properly solve a FE model, it is required to build a model with boundaries extending well beyond the critical region of the problem, so that the boundary constraints do not impact results. Inevitably, the model is built and computed to solve results for a great number of nodes that of are no practical significance to the designer.  There have been various techniques that have been developed to mediate this issue. For example, many codes include mesh generators that automatically construct a coarser mesh for the regions more distant than the zone of interest. Although, these techniques serve to reduce the problem, they do not eliminate it altogether. An alternative and more efficient approach that can be undertaken is treating only the problem boundary, (i.e. the ground-support interface), rather than creating a mesh that represents the entire surrounding ground. One example of such a method which is commonly used in practice is the Winkler method (Winkler, 1867), in which the ground medium is represented by a series of independent springs. Common applications that rely on the Winkler method include foundation slabs and piles. One drawback of this approach is that the springs lack continuity, thus in some instances results may be highly inaccurate. 3  The convergence-confinement method for tunnel design has been developed based on the Winkler approach. This method serves both as an important academic tool for the conceptual understanding of tunnel-support interaction, and also as a practical and efficient design tool. As will be reviewed and discussed in Chapter ‎2, this method is limited to a series of simplifying assumptions. This thesis introduces a new and general "Equivalent Boundary" approach, developed to solve a number of ground-support interaction problems using an approach. In this process, each type of problem is simplified so that the interaction of the ground and the support are fully accounted for, without the need of a more complex numerical analysis. This is achieved by deriving expressions for structural entities with an equivalent stiffness to the ground boundary.  For the simpler ground-support interaction problems, the method can be viewed as an extension of the convergence-confinement method. For the more complex tunnel problem, the method can be viewed as a modification of the convergence-confinement method. Practical applications of the method, including non-circular tunnel analysis, and cost probabilistic analysis for tunnel support cost estimation, are also presented.  1.2. Research Objectives The primary research objective is to develop a simplified analogous method for solving ground-support interaction problems. This is achieved by considering the following tasks: 1. To derive equations for assessing the room-pillar interaction problem. 2. To develop equations for assessing the shaft-liner interaction problem. 3. To develop a novel and simplified method for tunnel-support interaction analysis, with fewer limitations than the currently available methods.  Within the discussion regarding the primary objective, a number of secondary objectives are met as well: 1. To provide a conceptual framework for the static and dynamic analysis of pillars in room and pillar mines. 4  2. To develop a methodology for tunnel support cost estimation at the pre-feasibility project stage. This methodology serves to demonstrate the advantage of the method developed for this thesis, but also constitutes a stand-alone concept. 3. To lay the foundations for the development of a more generalized methodology to solve more complex problems.  1.3. Thesis Organization The thesis is divided into 8 chapters; following this introductory chapter, are a chapter summarizing the literature review and a chapter discussing the general methodology. The final chapter summarizes the conclusions and recommendations. The 3 research chapters in between, follow the structure displayed in Figure 1.  Figure 1 Thesis structure, main research chapters  First, the Equivalent Boundary, EB method is demonstrated for solving the simpler cases of pillars and circular shafts that can be reduced to single-degree-of-freedom (SDOF) problems. Next, the case of tunnels under non-uniform stresses, which involves a multi-degree-of-freedom (MDOF) solving process, is dealt with. (The explanation of the terms SDOF and MDOF is provided in Section 3.1). Finally, an extensive analysis is presented in Chapter ‎7. There, through a practical example, the implementation of the proposed method for tunnel-support analysis is given, in 5  conjunction with a proposed probabilistic analysis for tunnel support cost estimation. This analysis highlights the advantage of the EB method, which allows quick and simple analysis of varying ground conditions and alternative support systems.  6  2. CHAPTER 2  LITERATURE REVIEW From ancient times, civilizations have built underground tunnels and caverns for various purposes. Presumably, the process of planning and construction relied mainly on basic rules of thumb gained by trial and error processes. The earliest recorded formal empirical classification system developed for tunnel support design is attributed to Ritter ("Die statik der tunnelgewölbe", 1879). Since this time , a constant effort has been made to develop more improved and accurate design methods, including empirical, analytical and numerical methods. In addition to the obvious goals of guaranteeing the immediate safety of the workers, and the long-term stability of the underground facility, support is a major component of the costs associated with underground construction. Support costs can vary greatly, depending on the geological conditions (Hoek et al., 2006). Currently, studies have shown that cost over-runs in tunnel projects are a global phenomenon. For example, Efron and Read (2012) examined 158 tunnel projects from 35 countries, and concluded that in every case the final cost was higher than the initial estimation. Considering the high cost and risks of underground projects, it is essential that reliable analysis tools are made available. Amongst the empirical methods, the Q (Barton et al., 1974) and RMR (Bieniawski, 1989) classification systems are commonly used in practice for tunnels and slopes. In these methods the rock is rated using both quantitative parameters (such as the uniaxial rock strength), semi-quantitative indexes (e.g. RQD, Rock Quality Designation, Deere et al., 1967), and more qualitative parameters, such as the joint conditions. According to the estimated rating, the classification systems provide guidelines and recommendations for support selection and implementation. These methods rely on extensive databases and have been used successfully in many projects.  However, empirical methods are limited by their integration of their contributing factors and do not account directly for the stresses and deformations of the ground and support system. Hence, the empirical methods cannot be used for different applications, such as comparing 7  and optimizing the support system for given ground conditions  or running probabilistic analyses that account for varying ground conditions.  In the context of rock mass rating, Hoek et al. (2013) argued that although empirical methods have played a critical role in the development of practical design tools for rock engineering, researchers should turn their attention to actual physical properties of the rock, and to  numerical modeling of rock fracture networks, in order to develop and apply a better understanding of rock mass behavior. This statement is apparently valid for all applied sciences; if applied correctly, the advancements from qualitative to quantitative approaches increase objectivity and reduce errors.   In the following sections, a brief review of the current analytical and numerical methods is given, and their advantages and shortcomings highlighted.  2.1. Analytical and Analogue Methods Analytical methods can be defined as precise mathematical formulations for application to problems in mechanics. During most of its historical development, the science of mechanics of materials relied principally on mathematical theorists. Much of their work represents mathematical skills and intuitions of a very high order, challenging even for the most advanced researchers of today.  In the context of soil and rock mechanics, Coulomb wrote important and fundamental works as early as 1773. A notable analytical method used for tunnel analysis is the solution derived by Kirsch (Kirsch, 1898). Kirsch’s closed-form solutions can be used to calculate stresses and displacements around a circular opening.  The solution is obtained using a complex-variable method based on Airy's stress function (G.B. Airy, 1801-1892). The stress function is derived so that it satisfies the conditions of static equilibrium, stress and strain compatibility, and problem boundary conditions. Analytical methods are usually limited to simple and ideal conditions due to the complexity of the mathematics required for treating complicated geometries and constitutive relations for nonlinear materials. Similarly, the Kirsch equations are limited to simplified ground 8  conditions: homogenous, elastic, and isotropic material behavior. Moreover, the Kirsch solution does not account for the ground-support interaction, and therefore can be used only for initial assessment purposes (Brady and Brown, 2004).   Various analytical solutions have been derived for the analysis of lined tunnels in elastic ground. One such notable method, which has been widely used by practical engineers in the past due to its simplicity, is Einstein and Schwartz's solution (1979). Einstein and Schwartz derived closed form solutions of tunnel-support interaction analysis for the estimation of liner stresses. The Einstein-Schwartz solution is limited to circular tunnels, in elastic ground, supported only with uniform liner applied at excavation. Carranza-Torres et al. (2013) modified the solution by Einstein and Schwartz in to a two-staged problem in order to account for the effect of delay in installation support.  A list of other available analytical solutions for lined tunnels can be found in the Technical Manual for Design and Construction of Road Tunnels (AASHTO, 2009). These solutions have the advantage of enabling quick and simple implementations. However, the mathematics for deriving these solutions is relatively complex.  Another approach which can be categorized as analytical, and is widely studied and used, is the convergence-confinement method. The convergence-confinement method is based on three main independent components: 1. The tunnel convergence curve. 2. The support curve. 3. The Longitudinal Displacement Profile. Numerous solutions for tunnel convergence curves, commonly termed as the Ground Response Curve (GRC), as well as for the support response, have been published (amongst them: Brown et al., 1983; Duncan Fama, 1993; Panet, 1995; Oreste, 2003; Vrakas and Anagnostou, 2014).  Once both the tunnel convergence and support confining curves are plotted, their point of intersection shows the final displacement of the supported tunnel. In addition, this point shows the relative contribution of load between the ground and the support. Thus, by 9  knowledge of the stiffness of both components (ground and support), the ground-support interaction process can be inferred. Figure 2 shows a simple example of a convergence-confinement plot, where the point of curve intersection defines the system-equilibrium displacement by its reflection on the horizontal axis, and the relative load bearing of the ground and support can be seen in the vertical axis.  In addition to finding loads and system displacement, the graphical representation provides an intuitive means of understanding the influence of different conditions, such as stiff vs. weak ground and support, and the delay of support installation, as later discussed in this Section.  Figure 2 Example of a convergence-confinement plot  An alternative representation of the ground-support problem is viewing the ground and support as an analogue model with two springs subjected to an external force, as shown in Figure 3(a). The stiffness of the equivalent spring that represents the combined ground-support system is obtained by addition of the two springs (𝑘𝑒𝑞 = 𝑘1 + 𝑘2). As shown in 10  Figure 4(b), the previous graphical representation is modified to include a single curve that has a slope equal to the sum of the ground and support slopes.   Figure 3 Ground-support problem analogue representation by: a) two springs in series,  and b) alternative graphical representation of combined ground-support curve (marked in green)    Although commercial programs that rely on numerical finite-element methods are not restricted to many of the simplifications assumed for convergence-confinement solutions, the latter are still widely used in practice as they are simpler to carry out and require reduced computational resources. Brown et al. (1983) presented a summary of such solutions with 11  respect to the material strength criterion, along with an additional solution of their own. These solutions include mainly elastoplastic Mohr-Coulomb behavior, and vary primarily by their post-peak behavior. Many additional solutions have been developed in the past years to account for more complex conditions, such as different failure criteria, and even 3D effects (Svoboda & Mašín, 2010; Cui et al., 2015).  Traditionally, for obtaining solutions of GRCs, axisymmetric conditions are assumed. The axisymmetric assumption includes two conditions:  1. The shape of the tunnel is circular. 2. The loading is hydrostatic (i.e. in-situ vertical stresses are equal to the horizontal stresses).  Violation of axisymmetry significantly increases the complexity of the problem. With respect to its practicality, assuming a circular tunnel can be viewed as reasonable, as many tunnels are circular or near circular. Conversely, the condition of hydrostatic stresses rarely exists.  Two important phenomena that occur only under non-uniform loading are arching and bending. Figure 4 shows stress contours and a tunnel displacement contour plot, for two tunnels; one under a hydrostatic stress field, and the other under a non-hydrostatic stress field. It can be seen that under a hydrostatic stress field the stress contours that form around the circular tunnel are parallel to the tunnel and uniform for a given radial distance from the tunnel center. Similarly, the tunnel displacement contour is constant and parallel to the original tunnel contour. 12   Figure 4 Comparison of stress results from FE models for a) hydrostatic stress field vs. b) non-hydrostatic stress field. The tunnel deformation deformed shape is exaggerated for the sake of illustration and appears in gray  In contrast, under non-uniform loading stress concentrations appear around the tunnel walls. The magnitude of the stress at the points of concentration are a function of the ratio of the initial horizontal to vertical stresses, which in geotechnical engineering is referred to as the K ratio. Note that if the horizontal stresses are greater than the vertical stresses, the location of stress concentration shift so that they occur at the tunnel roof and floor. The phenomenon of stress concentration around tunnels is sometimes referred to as 'ground arching'. Arching around tunnels has an advantageous effect in terms of support design, as the dead weight of the ground above the tunnel is transferred to the sides of the tunnel (Terzaghi, 1946). However, in deep tunnels in rock, where there is an extremely high stress environment, stress concentrations can result in failure of the rock. In severe cases, this failure can be manifested by spalling and rockbursting, and special considerations must be incorporated into the support system design and implementation process (Kaiser, 2011). In addition to the non-uniformity in stress distribution, tunnel displacement contours also show a non-circular curved pattern. When applying a tunnel liner the non-uniform displacement will induce bending moments within the liner. The implication of both these 13  phenomena, arching and bending, is that the stresses that develop in the tunnel support system may vary greatly once the stress is non-uniform. On the basis of the convergence-confinement method, the commercial program RocSupport was developed (RocScience, 2015). RocSupport is a software tool for estimating deformation in circular excavations in weak rock, based on the tunnel interaction with various support systems. The program allows selecting between various GRC solutions by different authors. Figure 5 shows an example of a convergence-confinement plot given by RocSupport.   Figure 5 Convergence-confinement curves (blue and purple, respectively) given by the program RocSupport, by RocScience (2015)  As can be seen in Figure 5, the confinement curve does not start from the origin; this is due to the consideration of the fact that support is regularly installed at some distance from the tunnel face, as illustrated at the upper right corner of the figure. This condition implies that some relaxation occurs prior to support installation, thus the final stresses and displacements are affected by this occurrence.  14  The convergence-confinement plot demonstrates how this effect is advantageous in terms of the loading on the tunnel support: the displacement and load imposed upon the support is reduced as the confinement curve is shifted to the right. Hence, lighter support systems are required. This point is one of the fundamental principles of the Sequential Excavation Method (SEM), or as it often is referred to as New Austrian Tunneling Method (NATM), where the displacements are recorded throughout tunnel excavation and the optimal point of support installation is found accordingly.  The other implication of support delay is that as the delay is increased, the final displacement increases as well. It is noted that in shallow tunnels under urban areas, settlement of the surface can cause damage to roads and structures, thus the support is installed as soon as possible to avoid the detrimental effect of support delay. In these cases, displacements, rather than support pressure, dictate the limiting criterion (Fang et al., 2012). In order to quantify the effect of support delay on the tunnel convergence process, the Longitudinal Displacement Profile (LDP) concept has been introduced and studied. LDP is a term used for describing the deformation curve from the tunnel face to the point of support installation. Knowledge of the LDP can assist in assessing the magnitude of initial displacements induced prior to support installation. During construction, monitoring displacements of the tunnel wall can help validate and/or calibrate the preliminary assessment, in order to find the optimal support installation timing. Perhaps the best method for accounting for the LDP effects is by 3D numerical simulation. However, 3D simulations are complex and time consuming; therefore different researchers have suggested formulae for determination of the LDP (Panet 1993; Carranza-Torres and Fairhurst, 2000; Vlachopoulos and Diederichs, 2009). Note that these methods all assume hydrostatic stress conditions. Chern et al. (1998) proposed an empirical solution based on curve fitting analysis to LDPs measured in-situ. Figure 6 shows a plot generated using Chern's equation. From this plot it can be inferred that the tunnel face ceases to affect the tunnel convergence from a distance of 4-5 tunnel radii. Note that Chern's equation is a function of the tunnel geometry alone, i.e. the ground properties do not influence the LDP. 15   Figure 6 Longitudinal displacement profile according to Chern et al. (1998). The ratio of radial displacement to the maximum longitudinal displacement is plotted against the ratio of the distance from tunnel face to tunnel radius.    According to the underlying assumptions of the convergence-confinement method, limitations are imposed upon practical applications for actual tunneling projects. This point is emphasized in the RocSupport program manual (Rocscience, 2015), which states: “The main assumptions in the RocSupport analysis method are as follows:  The tunnel is circular.  In-situ stress field is hydrostatic (i.e. equal stress in all directions)  Rock mass is isotropic and homogeneous. Failure is not controlled by major structural discontinuities.  Support response is elastic-perfectly plastic. 16   Support is modeled as an equivalent uniform internal pressure around the entire circumference of the circular tunnel. This last assumption in particular, should be carefully considered, when comparing actual tunnel behavior, and calculated results using RocSupport. The assumption of uniform support pressure implies that:  Shotcrete and concrete linings are closed rings  Steel sets are complete circles  Mechanically anchored rockbolts are installed in a regular pattern which completely surrounds the tunnel. Because this will not usually be the case, actual support capacities will be lower, and deformations larger, than those assumed in RocSupport. In spite of the limitations discussed above, rock-support interaction analysis has many attractions and, when used in conjunction with numerical analysis, it can provide valuable insights into the mechanics of rock support and reasonable guidelines for the design of this support.” Schurch and Anagnostou (2012) investigated the effects of deviations from the underlying axisymmetric GRC assumption. They found the range of conditions for which the axisymmetric assumption returns acceptable results, considering the extent of uncertainty associated with geological materials. However, their findings are of limited practical use, as the support response cannot be found in a straightforward process from the GRC. This is because the traditional expressions for support pressure are also based on the axisymmetric assumption. Some attempts have been made to provide solutions that are not restricted to axisymmetric conditions. Exadaktylos and Stavropoulou (2002) presented a closed-form solution for stresses and displacements around D-shaped tunnels. Their solution relies on complex potential functions and conformal mapping representation. The solution is restricted to linear-elastic material behavior.  For deep circular tunnels under non-uniform loading conditions, Detournay (1986) derived an approximate solution for the excavation of a circular tunnel in elastoplastic Mohr-Coulomb 17  material. Note, that this topic of research is still of interest; in a recent paper, Lu et al. (2018) published a similar solution for the widely used Hoek-Brown failure criterion.  Based on Detournay’s solution, Detournay and St. John (1988) developed a series of design charts that provide estimates of the average plastic radius and the tunnel convergence at the major axes (i.e. at the tunnel wall and roof). Their solution can be categorized as semi-analytical and relies on a number of a-priori assumptions. However, this solution cannot provide results which account for the actual tunnel-support interaction. Instead, they assume a uniform support pressure, and estimate the reaction accordingly. Consequently, results of the actual support pressures can be highly inaccurate, especially under high horizontal in-situ stresses (K>1.75).  Within the domain of analytical methods, a distinction can be made between direct solutions, and solutions based on analogues. Both the Kirsch equations and the convergence-confinement are examples of solutions which have been derived by directly addressing a given problem, assuming a set of conditions.  Contrastingly, analogue methods are derived by converting the original problem to a simpler equivalent system. While direct analytical methods may be simpler than FEM from the perspective of the user, the derivation process commonly requires a high degree of mathematical skill. By using proper analogues, the need for advanced mathematics can be circumvented, as a simpler and equivalent system may be considerably easier to derive and analyze.   Analogue systems can be precise or approximate. In precise equivalent analogues, the solutions provided by the equivalent system are exactly equal to solutions obtained by analytical solutions. A simple example is shown in Figure 7, where a truss structure is represented by an equivalent spring. In this example, the load to deformation relationship is precisely identical; therefore the system can be regarded as precise. The truss structure consists of many elements, and solving the deformation to a given load requires a lengthy calculation process. Contrastingly, the equivalent spring is a single element, and any resultant deformation can be found by simply dividing the force by the stiffness of the spring. The conversion process, from truss to spring, requires first solving the full truss problem, and then 18  calculating the equivalent stiffness, Keq, of the spring. This analogue may be extremely useful for the case where the truss is an archetype that varies proportionally by different variables (e.g. size, material, strength, external load, etc.). The equivalent stiffness Keq of the spring is a function of these variables. This would obviate the need for recalculation of the truss sub-elements for each variation.     Figure 7 Example of a simplification of a truss problem using a spring analogue  Another approach is to derive equivalent systems that yield approximated solutions. In these cases, it is necessary to define the conditions for which the analogue is applicable. A notable example of an approximate approach is known as Winkler springs. Similar to the concept of the convergence-confinement method for tunnels, Winkler (1867) developed a modeling approach that uses a series of discrete and elastic springs to represent the ground at the soil-structure interface (see Figure 7).   19   Figure 8 Foundation slab analysis using "Winkler" springs    According to this idealization, deformation of the foundation due to the applied load is limited to the loaded regions only. Furthermore, as the springs are discrete, the response of a single spring does not affect the neighboring springs, as opposed to reality where the response in one point in the ground has an impact on its near vicinity. In some cases this can develop into a state of structural failure. Due to these reasons, under some instances this simplification can result in unrealistic results. In such cases, it is necessary to model the problem with FEM, which assumes continuum behavior for the ground. The Winkler approach has been modified and adjusted by many researchers. In addition to the application of slab foundation analysis, the method has been extended to account for problems such as foundation piles subjected to lateral loads, and retaining wall analysis. As an example of the latter problem, the commercial program WALLAP (2007), which is widely 20  used by engineers is based on the Winkler approach, and is used for analysis of a variety of types of retaining walls, such as cantilever walls, anchored walls, and strutted excavations. Much effort has been invested in finding appropriate values for spring stiffness, or, as commonly referred to in the context of foundation slabs, the coefficient of subgrade reaction. The Filanenko-Borodich model (1946) for slab foundations allows connecting of the springs with a thin elastic membrane under tension in order to achieve continuity between the springs. This approach was further modified by Pasternak (1954), Kerr (1964), and others.  The Winkler approach was used for solving dynamic problems as well, such as earthquake loading. For this subject, Gazetas (1985, 1991) has published many papers for different types of structures such as embedded foundations , piles , and caisson walls (Gerolymos and Gazetas, 2006).   In order to account for plastic failure of the soil, different authors proposed methods for derivation of inelastic springs with non-linear stress and strain curves. By this means, a better representation of the actual elastoplastic ground response is achieved. For the study of lateral loads on foundation piles, extensive studies have been carried out in order to find proper values for the stress and strain curves, referred to in this context as p-y curves, for different soil formations (Peck, 1974). P-y curves have been incorporated into different building codes for foundation analysis. On the basis of the Winkler approach, Oreste (2006) developed a numerical method for tunnel support analysis termed the Hyperstatic Reaction Method (HRM). In this method, the tunnel boundary is represented by a series of independent Winkler springs (see Figure 9).  Equations are derived for the stiffness of the springs and their yielding limit. The HRM method is very suitable for conducting probabilistic analyses of the rock-support interaction, and does not require the assumption of hydrostatic loading. Yet, there are some drawbacks associated with this method:  1. The vertical and horizontal loads cannot be found in a straightforward manner, and empirical or back analysis methods must be applied for estimating these loads. 2. The nature of independent springs which lack continuity differs from that of the actual ground, or rock mass. For example, in the case of a point load acting on one point, the 21  response is affected only by the loaded spring, whereas in the case of an actual rock mass the neighboring area contributes to the resistance and to the load and affects the final strains and stresses. This point becomes critical when dealing with conditions that resemble the action of point loads. In tunneling, active bolts and loose wedges serve as such examples.   Figure 9 HRM model (modified from Oreste, 2006)  The attractiveness of using an equivalent beam for engineering analogues, rather than a series of equivalent and independent Winkler springs, stems from the continuous nature of the beam. Following this rational, Burland (1995) developed a method based on an equivalent beam for the assessment of soil-structure interaction when tunneling underneath an existing structure. Burland’s method can be categorized as an approximated analog system. Nevertheless, depending on the derivation process, methods that rely on a beam analogue can be developed to be precise.   22  2.2.  Numerical Methods The finite-element (FE) method was developed due to the practical need to solve complex structural problems. It is difficult to trace the exact origin of the method, and there are different publications from the 1940s and 50s that introduced the basis for this method. Alongside with the vast advancements in computer science, in the last few decades numerical methods have become an extremely powerful and useful tool for structural analysis. Nowadays, FE methods are commonly used in both practice and academic research. Many of the limitations that apply to analytical methods do not apply to numerical modeling. Numerical methods can model various 2D and 3D geometries, effects of staged excavations, dynamic problems, and additional conditions that are of such a degree of complexity that deriving an analytical solution would most likely be unachievable.  In the field of geotechnical design there are different numerical methods used for modeling. For the modeling of soils and other continuum materials, a detailed review of numerical methods and material failure models is given in a book by Anandarajah (2010). An extensive review will not be given here. In order to avoid misconceptions related to terminology in later references in this thesis, a number of notable methods are mentioned here.  Finite-element method (FEM): in the FEM the problem is subdivided into smaller parts, called finite elements. This is in contrast to analytical methods, where the mathematical equations assume infinitely small elements. FEM is an approximate method, which uses variational methods from the calculus of variations to approximate a solution by minimizing an associated error function. The accuracy of FEM solutions depends on the size of the mesh elements. There are other methods that rely on different mathematical formulations (such as the finite difference method, FDM), but the general concept described is similar.  Direct Stiffness Method (DSM): This method, also known as the Matrix Stiffness Method, is a one that makes use of the members' stiffness relations for computing member forces and displacements in structures. The direct stiffness method can be seen as a particular implementation of FEM. Historically, the DSM is a precursor to the more complex FEM method. FEM codes are derived using formulations based on the principle of minimum potential energy, whereas in the DSM only the force-displacement relationship is used as a 23  basis for solving problems. In this thesis, the DSM is used for solving the equivalent beam in the tunnel problem (Chapter ‎6).    Boundary-element method (BEM): in this method, partial linear differential equations, which have been formulated as integral equations, are used to solve the problem. Unlike FEM, only the boundaries of the problem are discretized, thus a significantly lesser number of elements are required and the method is more computationally efficient. However, this method is more limited, and bounded by many simplifying conditions. Note that although BEM and the Equivalent Boundary method developed for this thesis both share the term 'boundary', the methods are unrelated, and should not be confused with each other.  Distinct-element method (DEM): is a numerical method used to simulate the behavior of discrete bodies that make up the system. In this method the contacts or interfaces of the different bodies must be defined so that the interaction between the bodies can be computed. There are various types of DEM codes, including the Discrete Elements Method, Discontinuous Deformation Analysis, Particle Flow Code, and more. (UDEC theory manual, 2005). In DEM methods, dynamic solving algorithms are used to track relative movement between bodies. In soil mechanics, DEM methods can be used for simulating the ground as a series of discrete elements, thus replicating actual interaction phenomena at the soil particle level. In rock mechanics, a rock mass is represented as an assembly of discrete blocks, whereas in FEM the rock mass is treated as an equivalent continuum. There are different hybrid methods (FEM/DEM, or FDEM) which incorporate both techniques to allow a realistic simulation of brittle fracture driven processes (Elmo et al., 2013).  2.3. Relationship between Different Analysis Methods This section provides a general discussion regarding the relationship between the different approaches reviewed in the previous sections. In this context, the role of simplified methods, such as the method developed as part of this research, is discussed. Scientific knowledge is comprised of systematic thinking and generalizations. One means of both gaining and presenting scientific knowledge is through conceptual models. Conceptual models can be defined as the reduction of information into an efficient and ideal description, or framework. Conceptual models are often abstractions of things in the real world, whether 24  physical or social. A conceptual model, when implemented properly, should satisfy a number of objectives, including enhancing individual understanding, providing a basis of communication and collaboration, and serving as a basis for future modification. In addition to these goals, physical models must provide proper predictions of measured physical behavior.  Physical and mathematical models are usually composed of relationships and variables. Relationships can be described by operators, such as algebraic operators, functions, differential operators, etc. Variables are abstractions of system parameters of interest that can be quantified. Through the course of time, scientific theory and models tend to become more complex, and require increasing efforts to study and compute. The purpose of adding increasing model complexity is to treat a wider range of problems, as well as obtaining higher levels of accuracy. Figure 10 illustrates three types of evolutionary processes: a) biological, b) technological, and c) scientific analysis methods.        25  a)  b)  c)  Figure 10 Evolutionary processes in a) biology and b) technology c) scientific analysis methods. Images modified from public websites  Each of these processes demonstrates that as time progresses, simple states evolve into more complex systems. Yet, the role of the simpler systems is not absolute; in some cases, the simpler systems become superfluous, and the complicated systems drive the simpler ones to extinction; while in other cases, the simple systems prevail and even flourish.  Two contrasting examples from the technological world are cellular phones and hand calculators: usually, there is no incentive in manufacturing older versions of cellular phones 26  once they have no advantage over the newer devices. Thus, these older versions gradually disappear. In contrast, hand calculators are still produced and used in massive numbers. Moreover, technological efforts are still being invested in developing more advanced hand calculator devices and apps, although personal computers provide significantly more advanced capabilities. This is because there are many instances where it is more convenient to use a hand calculator than a personal computer. In biology, the human evolution process demonstrates both contrasting examples: while archaic human species (such as homo-habilis and homo-erectus) no longer exist, various species of apes (such as gorillas and chimpanzees) have survived with relatively minor biological changes. Similar to the evolutionary processes in the technological world, the human prototypes pose no advantages over the current human “version” (I.e. homo-sapiens), while existing ape species are more fit for living in their natural environment.   To summarize, it is evident generally that as long as older and simpler systems have a clear advantage in terms of durability and/or efficiency, even if this advantage is restricted to a limited purpose, they continue to exist and evolve.  This principle applies to scientific and engineering analysis methods as well. One such notable example in geotechnical engineering is Coulomb's model for retaining walls. Coulomb studied the problem of lateral earth pressures on retaining structures as early as 1776. Coulomb used the failing soil block as a free body in order to determine the limiting horizontal earth pressure (see Figure 11). Although geotechnical engineering has advanced incredibly, the equations based on Coulomb's simple theory are still in use today for actual retaining wall design, not only in educational textbooks, but in official building codes. Moreover, Coulomb's basic model is still the subject of research and further development, even after the advent of more sophisticated methods such as FEM codes.  27   Figure 11 Diagram of forces on a failing soil block used for Coulomb's theory for retaining wall analysis.   Similarly, Newmark's sliding block model (1965) is an analogy of a single discrete block on a slope, used for the prediction of the permanent displacement of a slope subjected to seismic induced accelerations (see Figure 12). In this example as well, FEM codes with dynamic solvers can provide solutions slopes subjected to seismic loads, with fewer limitations compared to the sliding block method. Yet, Newmark's simple analogy is still a subject of research, and has been studied, advanced and modified by many researchers. In addition, Newmark’s method is still used in current software, e.g. the program SLAMMER (Jibson and Randall, 2011).   Figure 12 Newmark's sliding block analogy for dynamic analysis of slopes  28  The principle of Occam's razor (attributed to William of Occam, 1287–1347), sometimes referred to as the principle of parsimony, states that simpler theories are better than theories that are more complex (Sober, 2015). This principle plays an important role in current physics, biology, psychology, and particularly relevant to modeling. The essential idea being that among models with roughly equal predictive power, the simplest one is the most desirable. Added complexity usually makes model difficult to understand and analyze, and can also pose computational problems, including numerical instability. As opposed to the general trend of scientific progress, the principle of Occam’s razor incentives researchers and developers to promote simple and user friendly methods and programs, so long that results are acceptable (see Figure 13).   Figure 13 General trend of scientific progress vs. the minimalistic approach known as “Occam's razor”  Due to the rapid development of electronic processors and the spread of numerical methods, the role of the analytical and empirical based methods has changed, but nevertheless remains important in geotechnical engineering. These different approaches should not be considered 29  as rivalry to each other, or even as alternatives, but rather as being complementary (see Figure 14).  In this regard, advances in one approach can be utilized to advance other approaches. For example, empirical relationships can be incorporated into FEM and analytical methods in order to assist with parameter selection. For example, the Hoek-Brown failure criterion (Hoek et al., 2002) has been integrated into commercial numerical codes and analytical convergence-confinement solutions. Another example of method integration can be found in common practice, where preliminary simplified analyses based on empirical and/or analytical methods, are used for determining assumptions and conditions for more detailed and complex numerical analysis. Oreste (2009) discusses the current practical roles of the convergence-confinement method in geotechnical engineering. These roles include:  A quick and efficient tool for preliminary assessment of required support: This may be very useful for feasibility and cost assessment, when there is not enough data available to justify the execution of sophisticated numerical models.   Verification of the validity of complex numerical models: the more the numerical models are complicated, the greater the chance there is for obtaining unrealistic results. It is therefore good practice to compare results with simpler calculations, in order to assess the validity of the model.  Sensitivity analysis of the effect of different parameters: geotechnical characterization of a rock mass is always approximate, as well as costly. It is therefore very useful to know for a specific problem which parameters have the greatest impact on the rock and support system, so as to be able to best direct further site investigations and/or laboratory test towards the most influential parameters.  Estimation of tunnel convergence: it is necessary to have a good estimation early in the design process of the convergence that the tunnel will display. This is useful to define the type and precision of the measurement instruments, and monitoring procedures. 30   Back analysis of the monitoring measurements: monitoring measurements of a tunnel can play an important role on the calibration of geomechanical parameters. Back analysis techniques require numerous models to be carried out, usually under a short time frame, and therefore simplified methods are useful for this purpose.  Optimization of the support system: simple methods that are easy to setup and quick to run allow quick optimization of the support classes from the economic point of view.     Figure 14 Viewing alternative analysis methods as complementary methods  With respect to project progression, the role of simplified solutions, including the method proposed in this thesis is presented in Figure 15.  Simple analysis methods enable an efficient means of coping with geological uncertainty, prior to execution of more sophisticated models.  In this context, shortage of data at the beginning of a project (i.e. data uncertainty) makes it difficult to describe the observed variability of key geological parameters, and in the presence of data uncertainty even the most advanced FEM models would yield inaccurate results. As a project moves from prefeasibility through to detailed design, the amount of data collected will increase as efforts are made to minimize uncertainty and reduce risk. Design initially has to be completed based on a limited amount of core samples. Whereas these methods do not 31  overcome the issue of data uncertainty, they allow for multiple scenarios to be simulated in a shorter period of time, thus it is possible to identify key support design scenarios, and at the same time use the results to drive the data collection process required by sophisticated analysis methods.    Figure 15 The role of uncertainty, from data collection to complex numerical analysis methods    Much of the contemporary work carried out in the various fields of applied science focuses on developing and utilizing the most sophisticated and advanced computational tools to study highly complicated engineering problems. In terms of the evolutionary processes shown in Figure 10, this would imply that progress is viewed as advancement to the right side, rather than "moving backwards", to the left.  Simler and Hanson (2017) discuss the nature of academic progress and compare it to the course of a typical conversation. During conversations, each side responds with an associative matter relevant to the topic of discussion. Hence, the conversation progresses so that the focus of interest gradually shifts in a quite uncontrolled process, without the discussed topics being fully exhausted. Similarly, the topics of academic publishing focus on “hot topics”, rendering certain fields of research from being fully studied. Therefore, a 32  researcher that does not follow this trend, but publishes work on “old” topics, as important that his work may be, he will likely draw less attention and gain fewer citations.  This observation highlights the difference between the evolutionary processes that occur in the biological world and scientific progress. Biological processes have a mechanistic nature and develop strictly according to natural selection. Conversely, scientific progress is impacted by human psychology and sociology. Thus, the course of scientific development is prone to human biases, which are not necessarily compliant with actual needs and demands, and with the maximum benefit to society.   The implication of this observation is that it is important to contemplate on the course of scientific progress, so that research efforts are allocated correctly. In this regard, it is the author's opinion that concurrent to the conventional trend in science, it is important not to neglect fundamental theories, and to continuously develop and advance methods that are simpler to derive and implement, and consequently also easier to comprehend.  In terms of derivation, the process of deriving a code for FEM modeling is complex, both in terms of the mathematics involved, and coding of the computational processing and post-processing stages. It is important to provide students with simpler methods that can be fully developed on their own, with a minimal number of "black boxes". With respect to comprehension, sophisticated tools are extremely difficult to grasp. Per contra, analogous methods in engineering provide an ideal academic tool, as they are based on intuitive thinking. The term "academic tool" may connote to some as impractical, or hypothetical. A notable quote made by Kurt Lewin, the founding father of social psychology, is that "nothing is so practical as a good theory" (Lewin, 1945).  While this statement is valid for all fields of practice, it is argued that “good theory" is vitally important for managing geotechnical projects. In many projects, unanticipated conditions are encountered and the engineer is forced to make quick and crucial decisions. In these moments, it is the basic concepts that the engineer has grasped during his training and practice that guide him, or her, towards the practical solution. At this point, the subtle conclusions drawn from highly sophisticated analyses, which may be very useful throughout 33  the early design stages, cannot be part of the decision process, which has to be executed in a short period of time, and often under high pressure. In the context of this work, the proposed methodology, being an analogue method, provides not only an efficient means for solving ground-support problems, but also assists in deepening the understanding of the concept of ground-support interaction by using conceptual framework and analogues.   34  3. CHAPTER 3 METHODOLOGY 3.1. General Methodology The Equivalent Boundary (EB) method developed in this thesis deals with two types of problems: single-degree-of-freedom (SDOF) problems, and multi-degree-of-freedom (MDOF) problems. The number of degrees of freedom in a numerical problem is determined according to the number of different types of motion permitted for each node (vertical/horizontal/rotational). The transition from a SDOF to a MDOF problem greatly complicates the nature of the problem. For example, the convergence-confinement method assumes axisymmetric conditions, which is synonymous with a SDOF problem, and the problem can be solved using a series of simple equations.  Conversely, for solving MDOF problems, a stiffness matrix must be assembled, and linear algebra mathematics must be employed in their solution. In addition, the graphical representation of the convergence-confinement curves, as shown in Figure 5, are only suitable for the SDOF condition, as they describe the relaxation and loading of a single point.   The software MATLAB (Mathworks, 2012) is used in this research to develop an EB method suited for MDOF problems, Examples of the implementation of MATLAB for finite element analysis can be found in Ferreira (2008).  Soil and rock are both natural materials; therefore, a significant extent of the geotechnical engineering research and literature has to address topics such as site data collection, laboratory tests, and analysis of failure models. These topics are not explicitly considered in this thesis, as the main objective is the mathematical analysis of the ground-support interaction process. The proposed method assumes that proper characterization and analysis of the material (rock or soil) has been made prior to using the method. Note that this thesis adopts the word “ground” to indicate a generic material, whether rock or soil. However, for specific applications references are given to either soil or rock materials.   The proposed method assumes is valid only where the behavior of the rock mass can be treated as an equivalent continuum. When the rock mass consists of dominant structural 35  features (rock discontinuities) that violate isotropic conditions, an equivalent continuum approach should not be used (Hoek et al., 2002). The same assumption applies for soil materials with fissures or other types of discontinuities.  The FEM is the most common and simple modeling approach for equivalent continuum problems. With respect to the numerical analysis, it is important to recognize the difference between model verification and model validation, which can be defined as:   Model verification involves checking if the mathematical approach that is used is correct and consistent with the model, usually performed by comparing with other closed-form solutions and/or different available methods of calculation.   Model validation involves checking if the model correctly represents the physical reality, by comparing with actual structure measured response from laboratory and/or in-situ tests.  The goal defined for this thesis is developing an alternative solution method. Therefore, the proposed model must be verified by comparing the results to other available and scientifically established methods. Model validation would not change the outcome of the verification process. Agreement with field data would be a matter of calibration of the input properties due the inherent variability of natural materials. In other words, a verified model could be validated by running a sensitivity analysis on the material inputs, but a validated model is not necessarily an index of verification of a mathematical formulation. Accordingly, for method verification, numerical models are carried out using the FE program RS2 (RocScience, 2015).  3.2. Basic Formulation The basic idea of the proposed method is to simplify ground-support interaction problems by converting them into equivalent systems using appropriate structural entities to represent the ground at the ground-support boundary.  The structural analogues chosen for each problem are listed in Table 1. These choices are made according to the nature of the support in each problem:  36   A pillar typically consists of a single load and deformation in the vertical direction; hence it is analogous to the action of a spring.  A shaft liner is typically subjected to a uniform radial displacement; hence it is analogous to the action of a ring, or cylinder.  A tunnel liner is typically subjected to non-uniform pressures, which induce a coupled axial-bending response, hence the liner is best modeled as a beam (or beam elements in FEM).   Table 1 Structural analogues for different ground-support problems Structural analogue for the "Ground" "Support" "Ground"  Spring Pillar Room  Pillar problem Ring (cylinder) Liner Circular vertical excavation Shaft problem Circular beam Liner Circular horizontal excavation Tunnel problem  The key mechanism that dictates the ground-support interaction process is the distribution of gravitational loads between the ground and the support. By converting the ground to an idealized entity of the same nature of the support, the calculation process becomes significantly simpler; consequently, the combined ground-support stiffness is obtained by simple summation; and, the load distribution can readily be found according to the stiffness ratio of each of the components.  In other words, it is much more complicated to formulate and calculate the distribution of stresses between two structural entities of very different nature. In geotechnical problems the ground is commonly idealized to behave as an infinite (or semi-infinite) medium, while the support is localized at the ground-support interface.  Thus, it is difficult to assess the relative stiffness, or relative load bearing contribution, of each of these very different components.  The analogous problem also obviates treating and solving the entire zone of influence affected by the excavation of the surrounding ground. This is achieved by assuming an equivalent structural element that represents the ground boundary. The equivalent stiffness must be 37  found in order to solve the problem. The fictitious forces acting upon the fictitious “ground” elements are found based on the displacements of the assumed unsupported ground. Finally, the fictitious forces are divided by the stiffness of the supported ground to give the final displacements of the supported ground.  One shortcoming which results from the problem reduction associated with the proposed method is that the results provided are limited to the ground-support interface. Conversely, in FEM a full solution of the state of stresses and displacements across the entire zone of influence is given. These results are most commonly displayed by contour plots. However, the most important information of interest usually lies within the ground-support interface. For most cases of pillar, shaft, and tunnel problems, the displacements and stresses at the ground-support interface dictate the overall factor of safety (FoS).  To illustrate this point, Figure 16 shows an example from a 2D FE tunnel model, subjected to hydrostatic loading. Three distinct domains are recognized, with varying importance from the user’s perspective:  1. The far field region- in this part of the model the initial in-situ stresses remain unaltered by the excavation. This zone has no practical implications from the user’s perspective. This domain is only necessary to prevent the boundary constraints from erroneously affecting the model. It can be noticed that in contrast to the degree of importance, this portion of the model is largest. 2. The zone of influence- in this region the resultant stresses are magnified compared to their initial condition. For the most part, this domain is of secondary importance to the user, as the maximum stresses and displacements occur at the ground-support interface. Assuming elastic conditions, this zone is usually equal to 3-5 tunnel diameters. Under plastic conditions, the depth of plastic yielding around the tunnel helps assess the required length of rockbolts. 3. The ground-support interface- the tunnel-support interaction is manifested at this interface. The maximal and/or minimal results occur at this interface and dictate the tunnel FoS. 38  Considering this distinction, it can be argued that over 99% of nodes computed in a simple FE model simulating ground-support problem, are superfluous. Nonetheless, the FE method cannot be solved without constructing a model that includes all three domains. This argument applies not only to tunnels, but to other ground-support problems well.  Obviously, for more complex problems, it is necessary to conduct a full simulation of the problem in order to identify more complex failure mechanisms that occur at some distance from the ground-support interface. (e.g. a pillar or tunnel which are in proximity to other structures and/or geological features). Nevertheless, it is necessary to distance the FE model boundaries for these problems as well.   39   Figure 16 Example of the major stress contours from a 2D FE tunnel model, which can be divided into 3 different parts: 1) the ground-support interface, 2) the zone of influence, and 3) far field region.  The following is a detailed overview of the formulation of the proposed method. A detailed derivation for each problem discussed in this thesis will be given separately in Chapters 4, 5 and 6.  The fundamental equation defining the relationship between the external force imposed upon a structure and the resultant displacement is:  𝑭 = 𝑲 ∗ 𝑫    (1) 40   where F denotes forces, K stiffness, and D displacement. It is noted that the transition from a SDOF to a MDOF problem complicates the nature of the problem: in the case of a SDOF problem the expression show in Equation 1 is a simple equation, whereas in MDOF problems F and D become nodal vectors and K is the global stiffness matrix. Obviously, to solve any problem, it is required that two out of three of the components from Equation 1 are known, so that the third unknown can be found. For common structural analysis problems, the starting point is that the external forces are known, or assessed; the stiffness K can be assessed based on the geometry and material properties of the structure, or a number of sub-structures; hence, the unknown displacements can be computed.  Conversely, the starting point for an EB solution is exactly opposite: F and K are unknown, whereas the displacements are known, as shown in Table 2. Table 2 Starting point of known's and unknowns in typical structural problems compared to the EB method D K F   Unknown Known Known Typical structural problems 1 Known Requires derivation Unknown EB Method 2  More specifically, the forces due to gravity acting on the ground can generally be considered as known, assuming that the weight of the ground, the depth, and the coefficient of lateral pressure K are known, or can be properly estimated. However, for the proposed EB method, the forces that need to be found are the forces that would act on the equivalent ground structure and cause a response identical to that of the actual tunnel boundary. Those forces are fictitious, and therefore unknown; thus, they cannot be found directly with the knowledge of the gravitational loads alone.  Similarly, the stiffness of the boundary is unknown, as any structural element requires knowledge of the geometrical properties of elements (e.g. the beam cross section dimensions 41  for a beam element). In the case of the EB method, the elements that comprise the boundary are fictitious (as same are the forces); thus, the equivalent cross section is considered unknown. The displacement of the unsupported ground can be considered as known, based on either existing solutions, or solutions developed specifically for this purpose within this thesis. As shown in Table 2, by deriving a proper expression for the stiffness of the equivalent ground, the problem becomes mathematically solvable.  The methodology used for this thesis is to find the stiffness of an equivalent structural entity at the ground-support boundary. This initial process is required only for the derivation of the method and is not part of the computational process. Once established, the method can be used to solve any problem that conforms to the simplifying assumptions that underlie each specific problem.  It is important to make a distinction between the stiffness of a material, and the stiffness of a structure. The elastic modulus of a material is not the same as the stiffness of a component made from that material. Elastic modulus is a property of the constituent material; stiffness is a property of a structure (or sub-structure). Hence, it is dependent upon various physical dimensions that describe that component, which include the elastic moduli. That is, the modulus is an intensive property of the material; stiffness, on the other hand, is an extensive property of the solid body that is dependent on the material, its shape, and its boundary conditions (Baumgart, 2000).  Once the derivation is completed and the method has been established, solving each problem includes two main computational steps:  1. The fictitious forces acting upon the fictitious “ground” elements are found, based on the stiffness of the equivalent ground and the displacements of the unsupported ground.  2. The fictitious forces are then divided by the stiffness of the supported ground. This yields the displacement of the supported ground.  42   Figure 17 Summary of the EB calculation stages. Ffic is the vector of fictitious forces, Kg represents the stiffness matrix of the equivalent ground, Ks is the stiffness matrix of the support, Dus is the unsupported displacement vector, and Ds is the displacement of the supported ground.  These two steps are shown in Figure 17 in equation form. It is noted that for the multiple degree-of-freedom problem of tunnels, the equations must be arranged in matrix form (see Figure 34). With knowledge of the displacements of the supported ground, the stresses acting within the support components can be found based on the relationships between stresses and displacements, known from basic mechanical theory.  At this point, the selected support system can be assessed according to the project design criteria. The internal stresses that develop within the support must be acceptable, i.e. they must not exceed the maximum allowable stress. In some projects, maximum displacements may also be defined as limiting factors. As will be demonstrated in Chapter ‎7 through a practical example, this method allows quick assessment of various scenarios (i.e. varying ground conditions, and different support systems).     43  4. CHAPTER 4 APPLICATION OF THE EB METHOD FOR PILLAR PROBLEMS 4.1. Introduction In the design of room-and-pillar systems, the loading capacity of a pillar and the assessment of pillar’s factor of safety (FoS) have a significant economic impact, as they relate to the size of the opening and the extraction ratio. Pillar FoS is defined as the ultimate strength of the pillar divided by the stress acting on the pillar. The most generally accepted techniques for estimating pillar strength use empirical formulae based on survey data from actual mining conditions. Several different formulas for pillar strength can be found in the literature including  Hedley and Grant, 1972; Krauland and Soder, 1987; Sjoberg, 1992; and Martin and Maybee, 2000. In contrast to pillar strength, little effort has been devoted to investigation of the loading environment (Roberts et al., 2002). For initial analysis of the stresses acting on the pillar, the tributary area method is most commonly used (Brady and Brown, 2006). This method simply assumes that the vertical gravitational pressure acting on half the span of the mine is imposed on the pillar. The pillar and mine act together to resist gravitational deformations. This interaction is not directly considered by the tributary area method. Moreover, due to the heterogeneous nature of the rock mass, the mine and pillar strength depends on both the intact rock properties and the shear strength of the natural discontinuities. The concentration of stresses directed from the mine to the pillar and the transient stresses induced upon the pillar through the excavation process likely cause the pillar strength to degrade due to progressive brittle damage in the form of crack growth and shearing of pre-existing joints. In turn, the weakening of the pillar influences the distribution of loads between the pillar and mine. Oravecz (1997) developed an elastic analogue for determining stresses and displacements for pillars in coal mines, assuming that they were dictated by the compressibility of the seam.  The spring-analogue presented in this thesis examines the mine-pillar interaction in terms of the deformation of the mine-pillar system and is not limited to the case of stratigraphic units and seams. 44  Esterhuizen et al. (2010) used the concept of ground response curves for presenting mine-pillar interaction. In their work they used numerical analyses with the code FLAC (Itasca, 2016) to obtain pillar deformation results. The results were plotted in the form of mine-pillar interaction curves where the initial vertical pressure acting on the pillar was used as a reference point for plotting the vertical axis. In a similar context, Barczak et al.. (2005) used numerical modeling to obtain ground response curves specifically developed for longwall tailgate standing support design. The approach proposed in this thesis is to develop simple equations that would obviate the need to build a full numerical model for pillar analysis. Elastic conditions are assumed, as it is argued that for elastic-brittle rock masses, where the peak elastic stress is followed by a rapid drop in material strength, an elastic analysis would be sufficient for estimating the ultimate pillar strength and pillar FoS. However, for deformable plastic rock mass the yielding will likely initiate simultaneously in both the mine and pillar, therefore it would be recommended to conduct fully elasto-plastic numerical analyses in order to obtain the characteristic curves for the mine and pillar loading system rather than to rely solely on the equations proposed in this thesis.  A number of additional simplifying assumptions were made: 1. Plane strain conditions apply (i.e. the model is infinite in the out-of-plane direction. 2. The in-situ stresses are constant. 3. No high horizontal stresses are present (i.e. the horizontal-to-vertical stress ratio, k, is less than one).  4.2. Derivation Given a rectangular excavation, (see Figure 18), at a depth D below ground surface, and a rock material with unit weight 𝛾, the vertical pre-mining pressure 𝑃𝑂 acting on the elevation of the excavated roof is: 𝑷𝑶 = 𝜸 ∗ 𝑫   (2) Equation 1  45   Figure 18  Rectangular excavation with pillar at its centre  The deformation of the mine-pillar system is affected by both the stiffness of the pillar, 𝐾𝑃, and the stiffness of the mine, 𝐾𝑀. In order to simplify the mine-pillar system to a spring system with a single degree of freedom, the distributed pressure 𝑃𝑂 is substituted by an equivalent force, 𝐹𝐸𝑄. The equivalent force is defined as the force that would yield the same maximum displacement ∆𝑀 at the centre of the roof of the mine system alone (i.e. assuming the absence of the pillar), as shown in Figure 19.    Figure 19 Equivalent force 𝑭𝑬𝑸  46   In order to find  𝐹𝐸𝑄, a series of FEM models with varied parameters of the given system were constructed using the program RS2, and by a trial and error process  𝐹𝐸𝑄 was varied until the midspan displacement was equal to ∆M. Accordingly, the equivalent force 𝐹𝐸𝑄 was found, and can be approximated as: 𝑭𝑬𝑸 = 𝟎. 𝟐𝟒 ∗ 𝑳 ∗ 𝑷𝑶  (3) Equation 2  where L is the mine span. The mine span is equal to the width of the pillar plus the width of two rooms, and can be related to the extraction ratio, which is the ratio of extracted material to the total dimensions of the pillar and rooms. Note that 𝐹𝐸𝑄 is in units of N/m, as plane-strain conditions were assumed. The mine displacements are dependent and linearly proportional to the length of the unsupported span, the vertical pressure and the rock mass Young’s modulus. An equal displacement occurs at both the roof and floor level of the mine (according to the aforementioned constant stress assumption); hence the mine convergence is twice the displacement ∆M. The mine displacement ∆M is equal: ∆𝐌=𝑳𝑷𝒐𝑬𝑴  (4) Equation 3  where L is the mine span and 𝐸𝑀 is the rock mass Young’s modulus of the mine. Note that the maximum displacement of the mine roof and floor ∆M is calculated based on elastic analysis. The mine stiffness 𝐾𝑀 is therefore: 𝑲𝑴 =𝑭𝑬𝑸∆𝐌= 𝟎. 𝟐𝟒𝑬𝑴 (5) Equation 4  Equation 5 shows that the mine stiffness is not dependent on the excavation geometry and is approximately one fourth of the rock mass deformation modulus.   The stiffness of the pillar 𝐾𝑃 is a function of the rock mass Young’s modulus of the pillar, 𝐸𝑃, the pillar height H and the pillar’s cross section. In a plane strain analysis the pillar cross 47  section is affected only by the pillar width, W. According to elastic theory, the stiffness of the pillar is: 𝑲𝑷 = 𝑬𝑷𝑾𝑯  (6) Equation 5  The idealization of the mine-pillar system to a discretized spring system is illustrated in Figure 20. Initially, the system is converted into a system with two degrees of freedom, but owing to symmetry the problem can be further simplified and treated as a SDOF problem. Hence, for the stiffness of the single spring the pillar height is taken as half the original height, and the stiffness is therefore doubled. As the mine-pillar system is represented by two springs with deformation compatibility, the deformation of the pillar ∆𝑃 is affected by the combined stiffness of the mine and the pillar, thus: ∆𝐏= 𝑭𝑬𝑸 (⁄ 𝑲𝑴 + 𝟐𝑲𝑷) = 𝟎. 𝟐𝟒𝑳𝑷𝑶 (𝟎. 𝟐𝟒𝑬𝑴 + 𝟐𝑬𝑷𝑾𝑯)⁄   (7) Equation 6    Figure 20 Mine-pillar system idealization to a spring model  48  Equation 7 shows that the pillar displacement can be estimated using the rock mass Young’s modulus of the pillar and of the mine, the pillar dimensions, mine span (or extraction ratio) and the excavation depth. A graphical representation of this process in the form of a mine-pillar interaction curve is presented in Figure 21, which can readily be constructed based on the outputs from the above equations.  Figure 21 Proposed mine-pillar interaction curve  In order to compare results obtained by different methods, the pillar maximum vertical displacement was computed using i) the proposed spring analogue; ii) FEM analysis using RS2,  and iii) the tributary area method. The properties and geometry listed in Table 3 were used for this comparison, with all properties kept constant, with a constant width-to-height ratio W:H=1, and the span L varied in increments of 2 m. The initial stress field was taken as Po=9 MPa, corresponding to a depth in the range of 330-345 m, according to the relationship given in Equation 2, for typical rock unit weights (𝛾 =0.026-0.027 MN/m3).  49  Table 3 Properties for mine-pillar parametric study Po [MPa] 9 W [m] 6 H [m] 6 Erm Pillar [GPa] 5, 10, 15, 20, 25 Erm Mine [GPa] 25 L [m] 26, 28, 30, 32, 34  The FEM models were built in RS2. Figure 22 shows the mesh for the model with a constant span of L=26 m. The external boundaries of the model were set to approximately 5-6 times the span length in order to avoid boundary effects. A horizontal to vertical stress ratio of 0.33 was assumed for all models. For these models, the rock mass Young’s modulus Erm of both the pillar and mine was assigned an equal value of 25 GPa.    Figure 22 Mesh for numerical model built in RS2 for L=26 m.   The results are using Equation 7 were plotted graphically in Figure 23, showing that the proposed method agrees well with the results from the FEM analysis, the slight variation being interpreted as due to mesh effects.  50   Figure 23 Comparison of pillar deformation computed by the tributary method, FEM analysis, and the proposed spring analogue   An additional parametric study was conducted using the same properties and assumptions, assuming a constant span L of 32 m, and varying the Young’s modulus of the pillar, 𝐸𝑃, from 5-25 GPa. Calculations were made for width-to-height ratios of W:H=1 and 2. Figure 24 shows the percentage of error in calculated displacement using the tributary method for pillar width-to-height ratios of 1 and 2 respectively. The tributary area method is shown to be highly over-conservative, as it neglects the favorable effects of the stiffness of the mine which resists displacement together with the pillar. This is shown to be true even for the case where the Young’s modulus of the pillar and mine are taken as equal (Ep= Em); as the Young’s modulus of the pillar compared to the mine decreases, (Ep/Em<1), the error in the calculated displacements using the tributary method increases exponentially.  51   Figure 24 Percentage of error for displacement calculation using the tributary method vs. ratio of pillar to mine rock mass Young’s modulus    4.3.  Comparison with Empirical Pillar Formulae In order to compare the proposed method with empirical formulae, a room-and-pillar mine study case is used. The Middleton mine is a classic square room-and-pillar mine located in Derbyshire, UK, mostly under a cover of about 100m (Elmo, 2006). The excavation is within the Hoptonwood limestone. The pillar Factor of Safety (FoS) is calculated using the empirical formulae proposed by Lunder & Pakalnis (1997) and compared to FoS results calculated using the proposed spring analogue equations. 52  In order to find the pillar FoS with the spring analogue equations, the allowable pillar stress 𝜎𝑃 can be related to the allowable pillar displacement ∆PA using the elastic relationship between the stress and deformation of a bar: ∆𝐏𝐀= 𝝈𝑷𝑯𝟐𝑬𝒑     (8)  Equation 7  Lunder and Pakalnis (1997) examined the stress distribution in hard-rock pillars in Canadian mines and proposed that the average confinement in a pillar could be expressed in terms of the width-to-height ratio (W:H): 𝐂𝐩𝐚𝐯 = 𝟎. 𝟒𝟔 [𝐋𝐨𝐠 (𝐖𝐇+ 𝟎. 𝟕𝟓)]𝟏.𝟒𝐖𝐇  (9) Equation 8  where Cpav is the average minor/major stress ratio within the pillar core. Numerical models have shown that pillar confinement depends on the ratio between the far-field horizontal stress; however, for pillar ratios less than 1 the effect of k can be ignored. The confinement formula introduces a mine pillar friction term (k) calculated from the average minor/major stress ratio within the pillar core according to:  𝛔𝐏 = 𝟎. 𝟒𝟒𝐔𝐂𝐒(𝟎. 𝟔𝟖 + 𝟎. 𝟓𝟐𝐤)  (10) Equation 9  𝐤 = 𝐭𝐚𝐧 [𝐜𝐨𝐬−𝟏 (𝟏−𝐂𝐩𝐚𝐯𝟏+𝐂𝐩𝐚𝐯)]   (11) Equation 10  Equation 10 could be used as an input to Equation 8 to establish the limiting pillar displacement, which would correspond to the equilibrium point indicated in Figure 21. The average intact rock strength in the mine is 43MPa, and the corresponding values of cohesion and internal angle of friction are 9 and 45 degrees respectively (Mohr-Coulomb constitutive criterion).  53  Assuming a unit weight of 26kN/m3, Figure 25 below shows the corresponding FoS calculated using Equations 8-11 above for two pillars: i) 3.5m wide, 7m high (W:H=0.5); and ii) 7m wide, 7m high pillar (W:H=1). The room width is 7m for both scenarios.   Figure 25 Lunder and Pakalnis (1997) pillar formula and expected FoS for i) 3.5m wide, 7m high; and ii) 7m wide, 7m high pillar.   Table 4 Comparison of pillar factor-of-safety using the Lunder and Pakalnis and the proposed spring analogue  Pillar FoS Lunder & Pakalnis Spring analogue W:H 0.5 3.2 3.5 1 5 5.4   54  Table 4 shows a comparison of the FoS for the two pillars (W:H=0.5 and 1) calculated with the Lunder and Pakalnis pillar formula and  results using the proposed spring analogue equations. The results are generally in good agreement, with the FoS given by the spring analogue equations approximately 10% higher.  This difference can be explained considering that the Lunder and Pakalnis formula was originally developed for actual mine conditions including multiple pillar layouts, whereas the proposed method was derived using a single pillar layout. To extend the current approach to consider the effects of a multiple pillar layout would require a more complex derivation that involves a multiple-degree-of-freedom solution scheme for the interaction of a slab (representing the mine) and a series of springs, representing the pillars. Whereas the single pillar concept may be a limiting factor, it is argued that in the literature the numerical analysis of pillar behavior is typically limited to the case of single pillars under uniaxial loading conditions (Jaeger, 1969). The objective of this thesis is to introduce a simple concept first and study its validity when compared to common numerical and empirical applications used to characterize pillar stress and pillar strength.  4.4. Influence of Support on the Idealized Mine-Pillar System In underground mines, in order to maximize production, it is often required to consider adding support, most commonly in the form of rock bolts. Active support systems are implemented by pre-tensioning mechanical bolts, whereas grouted bolts and/or shotcrete act as passive support.  The effect of the passive support (shotcrete and/or bolts) installed in the mine roof and/or floor is analogous to the addition of a spring to the mine-pillar spring analogue. The slope of the mine stiffness curve increases in response to the support resistance that is added to the pillar and mine resistance, as shown in Figure 26a; Equation 6 can thus be modified accordingly: ∆𝐏= 𝑭𝑬𝑸 (⁄ 𝑲𝑴 + 𝑲𝑷 + 𝑲𝑺)  (12) Equation 11  where 𝐾𝑆 is the stiffness of the support system. In order to check the influence of bolts, numerical models using densely spaced and high strength bolt properties showed that passive bolting has a negligible influence on the mine stiffness and the overall pillar FoS. 55  When adding active support to the mine roof the slope of the mine stiffness curve remains the same, but is shifted parallel and downwards, as shown in Figure 26b. The magnitude of the parallel shift is 𝑁𝑒𝑞, where 𝑁𝑒𝑞 is the equivalent force for pre-tensioning pressure induced by the active support. Using the relationship between a distributed pressure and an equivalent force (see Equation 3), for a bolting pattern with a spacing s, and pre-tensioning force 𝐹𝐵, the equivalent active force 𝑁𝑒𝑞 is: 𝑵𝒆𝒒 = 𝟎. 𝟐𝟒𝑳𝑭𝑩𝑺     (13) Equation 12  Consideration of pillar confinement can be achieved by installing horizontal bolts and/or shotcrete on the pillar wall. The effect of pillar confinement on the mine-pillar interaction is shown in Figure 26c. According to the Mohr-Coulomb strength criterion the maximum pillar stress 𝜎𝑃 increases with an increase in the horizontal confining pressure 𝜎3. The slope ψ of this strength increase can be calculated using the friction angle ɸ𝑟𝑚 of the rock mass, and is: 𝒕𝒂𝒏𝛙 = (𝟏 + 𝒔𝒊𝒏ɸ𝒓𝒎 ) (𝟏 − 𝒔𝒊𝒏ɸ𝒓𝒎 )⁄  (14) Equation 13   56   Figure 26 Artificial support effect on mine-pillar interaction for a) passive roof support b) active roof support and c) pillar confinement. The black 'X' denotes pillar failure.  The spring analogue equations presented in the previous section can be used to evaluate the effectiveness of pillar confinement vs. mine roof support. Table 5 lists the properties used for the calculations and Table 6 presents the results. In addition to an overall increase in FoS, Table 6 also indicates the increase in FoS per bolt allowing the effectiveness of the support method to be evaluated.    57  Table 5 Properties and geometry for parametric study L [m] 30 Em [GPa] 15 Ep [GPa] 10 W  [m] 6 H [m] 6 Pre-tensioning force [MN] 0.5 Spacing [m] 1 Friction angle 40 Allowable stress [MPa] 10e6  Table 6 Comparison of bolting mine roof with bolting of pillar on FoS at different depths  Depth [m] 50 100 200 Bolting mine roof FoS increase [%] 23.81 1.37 0.68 FoS increase per bolt [%] 0.40 0.02 0.01 Bolting pillar FoS increase [%] 21.66 FoS increase per bolt [%] 1.80   Given that with increasing depth the stresses on the mine roof increases and the pre-tensioning force remains constant, the effectiveness of the roof bolting would be depth dependent, whereas the effect of confinement is not. The results would have important implications, as for deep mines the contribution of the roof bolts to the pillar’s FoS appears to be negligible, whereas for shallow mines the roof bolting was found to significantly contribute to the increase of the calculated FoS for the pillar. The analysis would suggest that bolting of the roof is not cost-effective when compared to the bolting of the pillar in the context of increasing the pillar’s FoS. However, bolting of the roof is frequently carried out for the purpose of preventing the falling of loose wedges; In this case it would be useful to examine and consider the contribution of the roof bolts to the overall pillar FoS. 58  4.5. Discussion of Post-Peak Pillar Analysis The support provided by bolting of the pillar is limited by the bolt tensile capacity and by maintaining reasonable bolt spacing. In general, there is a financial incentive to conduct mining operations with high extraction ratios and with a pillar FoS close to one. In this case, degradation in pillar strength caused by cracking and shearing could trigger a sudden drop of FoS below unity that would be followed by a dynamic response. The dynamic motion will rapidly attenuate and subsequently a new equilibrium state is reached. Hereinafter these two states will be referred to as the dynamic phase and residual state, each necessitating separate analysis procedures and remedial actions. Arguably, the most effective means of modeling pillars in the complex post-peak phases is by explicit dynamic finite element/discrete element (FDEM) simulations. These would provide a more realistic representation of the mine-pillar system, and better capture the loads transferred to the support system. Elmo and Stead (2010) used the FDEM approach to incorporate fracture networks to study pillar failure. They found that failure in slender pillars is predominantly controlled by naturally occurring discontinuities, with wider pillars failing through a combination of brittle and shearing processes.  Generally, when a system is no longer in static equilibrium, velocity and acceleration components are introduced into the equilibrium equation: 𝑭(𝒕) = 𝒌𝒖 + 𝒄?̇? + 𝒎?̈?   (15) Equation 14  The simplified SDOF spring idealization introduced in the previous sections oversimplifies the phenomena that occur in the dynamic phase for several reasons:  Different portions of the mass components consist of different accelerations.  The vibrating motions are complex and consist of a series of modes of motion rather than a single simple harmonic.   During the dynamic phase the stiffness of the pillar is not constant, as the fracturing continues to degrade the stiffness throughout the dynamic phase.   Various damping mechanisms that attenuate the motion such as plastic deformations and fracturing cannot be simulated or predicted using a simple SDOF analysis. 59  Nevertheless, it is important to recognize the design stages that need to be addressed for post-peak analysis. For this, the spring idealization can be used to provide a framework for better understanding the design procedure. The following sections discuss considerations that should be made for pillar post-peak design, and pillar FoS in the residual state.  4.5.1. Dynamic Phase In the case of brittle failure of pillars, kinetic energy may be transformed into fracturing and spalling. The spalling of the pillar poses an immediate threat to mine personal operating nearby the pillar. In addition, the spalling process further reduces the pillar strength and stiffness. It is argued that in order to eliminate the spalling hazard, a support system consisting of wire mesh and/or shotcrete layer is necessary.  Esterhuizen et al. (2011) surveyed pillar failures in order to identify different modes of instability. They found that in hard brittle rock failure of the pillars was observed to be related to spalling hard brittle rocks, as well as shearing along pre-existing angular discontinuities or progressive extrusion of soft infill materials on bedding planes. In order to assess the support required for the dynamic impact of spalling, knowledge of the volume and velocity of spalling is required. Zipf (2011) discussed different design approaches for coal mining at shallow depth in order to prevent sudden violent collapse of room and pillar mines.  Mitelman and Elmo (2014) carried out FDEM simulations with the code ELFEN (Rockfield, 2007). They found simulated results of spalling of rock induced by blasting to be in agreement with physical tests in terms of spalling volume and spall velocities. In the context of pillar failure, publications of measurements of those is generally lacking, thus spalling results using FDEM simulations cannot be validated. More work for spalling prediction is required to allow proper design for the pillar dynamic phase, and current design would have to depend on engineering judgment. For such initial assessment, a dynamic FEM simulation of the spalling impacts such as illustrated in Figure 27can be carried out in order to assess if a given thickness of shotcrete or wire mesh is sufficient. It is recommended that models with a range of different spall sizes and velocities are undertaken due to the uncertainty of these.  60   Figure 27 Illustration of spalling impact model  4.5.2. Residual State Following the dynamic phase, the gravitational loads are redistributed between the pillar and mine. In order to guarantee stability of the mine-pillar system at the residual state, it is required to assess both the residual deformational state, the demand on the pillar, and the residual strength of the pillar.  Considering the SDOF spring idealization discussed in the pre-peak conditions, initiation of the dynamic phase would be analogous to a sudden replacement of the spring representing the pillar by a weaker spring. In this case the system would accelerate downwards and then oscillate around the new point of equilibrium. Figure 28illustrates the transition of the system from the pre-peak to post-peak state where 𝐹𝑃−𝑅𝑒𝑠 and 𝐹𝑀−𝑅𝑒𝑠 denote the portions of the load imposed on the pillar and mine in the residual state, respectively.  The maximum amplitude of motion in the dynamic phase is: 𝑨𝒎𝒑𝒍𝒊𝒕𝒖𝒅𝒆 = ∆𝑷−𝑹𝒆𝒔 -∆𝑷  (16) Equation 15  where ∆P−Res is the final pillar displacement in the residual state. The amplitude of motion is proportional to the potential energy of the mine-pillar system. Thus, for further empirical or numerical studies, it is suggested that the amplitude be used as an index for development of 61  relationships for pillar dynamic failure. ∆P−Res can be assessed by inserting a residual rock mass Young’s modulus E𝑷−𝑹𝒆𝒔 into the pillar equations presented in Section  ‎4.2. Depending on ∆P−Res, the stress acting on the pillar in the residual state is: 𝝈𝑷−𝑹𝒆𝒔 =𝑭𝑬𝑸(∆𝑴−∆𝑷−𝑹𝒆𝒔)𝑨𝑷∆𝑴   (17) Equation 16  where 𝐴𝑃 is the pillar cross-section.  The degradation in stiffness of the pillar during the dynamic phase has two contrasting effects on the residual FoS; a large reduction in pillar stiffness causes a larger final deformation, which in turn causes the residual load on the pillar to be smaller on one hand (see Figure 28). On the other hand, a large reduction on stiffness implies a significant reduction in pillar strength.   Figure 28 SDOF Mine-pillar representation of post-peak dynamic phase. 𝑭𝑷−𝑹𝒆𝒔 and 𝑭𝑴−𝑹𝒆𝒔 denote the portions of the load imposed on the pillar and mine, respectively. ∆𝑷−𝑹𝒆𝒔 denotes the residual deformation of the pillar.  62  It is emphasized once more that for the post-peak analysis the SDOF spring idealization serves only the purpose of conceptualization of the problem, and to provide a framework for interpretation of results of more advanced numerical analyses.      63  5. CHAPTER 5  APPLICATION OF THE EB METHOD FOR VERTICAL                                                                         SHAFT PROBLEMS 5.1. Introduction Vertical shafts are constructed for various purposes, such as part of a mine layout, launching tunneling or pipe jacking operations, ventilation shafts, waste storage, and more. Shafts can be of different shapes; however, a circular shaft has the advantage of acting as a ring, i.e. under ideal conditions the stresses that develop are uniform compressive hoop stresses. Shaft construction methods can be divided into two categories: methods where the ground is first excavated, and the support is installed afterwards and methods where the support is installed prior to excavation. The latter category includes methods such as secant pile driving and slurry walls. The shaft construction method reflects on the probable displacements and stresses that will be imposed upon the liner. When excavation precedes support installation, some ground relaxation may occur and consequently the stresses that develop in the liner are reduced, and the final displacements are larger. In a ground-support interaction diagram, such as the one shown in Figure 21 for the mine-pillar system, ground relaxation would be reflected by movement of the support reaction line to the right. Different authors have derived equations for vertical and circular shafts that take into account elasto-plastic ground behavior, including Terzaghi (1943), and Jaeger and Cook (1969). Wong and Kaiser (1988) used a convergence-confinement approach based on different modes of soil yielding. In a recent paper by Ozturk and Guler (2016), the authors presented a program for the calculation of liner thickness for circular shafts in rock, based on regression analysis of results of numerical analyses. In the following section, the vertical shaft problem is simplified and idealized as a double ring system as shown in Figure 29, with one ring which represents the ground, and the other representing the liner. The derivation for the shaft problem is similar to the previous process for the mine-pillar problem. 64   Figure 29 Shaft idealized into a double ring under uniform compression analogy  The proposed simplification incorporates the following assumptions:  The ground is homogenous and the material behavior is elastic.  The construction method is such that little or no deformations occur prior to support installation.  The ground loading is due to gravitational forces, and the initial horizontal stresses are equal in all directions.  Frictional effects between the ground and liner are not accounted for.  Radial stresses acting through the liner are negligible. This assumption becomes less valid for very thick liners.  The overall bending resistance can be neglected. The following section presents simple equations for the calculation of the hoop stresses and deformations that develop in the shafts based on the above assumptions. The stresses in the longitudinal direction cannot be found using this method. The shaft-liner interaction analogy is limited to the surface area, where the maximum displacement occurs, so that the system can be simplified into a SDOF problem. 65  5.2. Derivation and Results In order to solve the shaft problem using a convergence-confinement approach, two of the three components: stiffness, force, and deformation, must be known for both the ground and the liner. Based on elastic theory, the stiffness of a ring subjected to a uniform radial pressure is: 𝑲𝒓𝒊𝒏𝒈 =𝑬(𝒓𝟐−(𝒓−𝒕)𝟐)𝒓(𝟏−𝒗)[(𝟏−𝟐𝒗)𝒓𝟐+(𝒓−𝒕)   (18) Equation 17  where 𝐸 is the Young's modulus of the ring material, r is the radius, t is the thickness, and v is the material Poisson’s ratio. Thus, by inserting the liner parameters into Equation 18 the shaft liner stiffness 𝐾𝑳𝒊𝒏𝒆𝒓 can readily be found. The process for determining the shaft response curve is non- trivial.  Numerical analyses were conducted in RS2 to first study the displacement profile of the unsupported shaft. These models assume axisymmetric conditions, i.e. the shaft centreline coincides with the rotated vertical axis (see Figure 30). The models were built according to the assumptions mentioned in the previous section (i.e. elastic and homogeneous material, support installed prior to excavation, loading is applied by gravitational body forces). Analysis results showed that the maximum displacement for an unsupported circular shaft would occur at or close to the surface level (i.e. in a distance of 1-2 shaft diameters from the surface level) and would be independent of the shaft depth).    66   Figure 30 Unsupported shaft axisymmetric model in RS2. Hs is the shaft height, and Ds is the shaft diameter   Figure 31shows a plot of the shaft strains against the normalized shaft depth for two different Young’s moduli for the ground. Note that the shaft strains are defined as the radial horizontal displacement divided by the shaft radius, and the normalized shaft depth is defined as the shaft diameter Ds over the shaft depth, or height, Hs. Thus, the ground displacement ∆𝑮𝒓𝒐𝒖𝒏𝒅 can be found based on the Young’s modulus of the ground 𝐸𝐺𝑟𝑜𝑢𝑛𝑑 and shaft radius r alone: ∆𝑮𝒓𝒐𝒖𝒏𝒅[𝒎] =𝟓∗𝟏𝟎𝟓𝑬𝑮𝒓𝒐𝒖𝒏𝒅[𝑴𝑷𝒂]𝒓  (19) Equation 18  where ∆𝑮𝒓𝒐𝒖𝒏𝒅 is in mand 𝐸𝐺𝑟𝑜𝑢𝑛𝑑 is in MPa.   67    Figure 31 Shaft strain vs. shaft normalized depth   Next, an expression for the ground stiffness must be found. The ground is assumed to behave similar to the liner: as a ring in compression. Therefore, a ring is chosen for the ground structural analogue, and the problem is reduced to a double ring problem. For the ring stiffness, the unknown which must be found is the fictitious ring thickness. A trial and error process is carried out, by assigning values and comparing to results obtained in FE models. The equivalent thickness of the ground ring 𝑡𝐺𝑟𝑜𝑢𝑛𝑑  is then found to be: 𝒕𝑮𝒓𝒐𝒖𝒏𝒅 = 𝟎. 𝟒𝟓 ∗ 𝒓   (20) Equation 19  where r is the shaft radius. The ground stiffness 𝐾𝐺𝑟𝑜𝑢𝑛𝑑 can then be calculated using Equation 17 by inserting 𝑡𝐺𝑟𝑜𝑢𝑛𝑑 along with the ground Young’s modulus and Poisson’s ratio. The fictitious force is given by the product of the displacement and the stiffness found in the previous two steps: 𝑭𝒇𝒊𝒄 = 𝑲𝑮𝒓𝒐𝒖𝒏𝒅 ∗ ∆𝑮𝒓𝒐𝒖𝒏𝒅 (21) Equation 16  68  The supported shaft displacement is then calculated by applying the fictitious force on the combined stiffness of the ground and support, similar to the step taken in Equation 6 for the mine-pillar spring analogue: ∆𝑺𝒉𝒂𝒇𝒕= 𝑭𝒇𝒊𝒄 (⁄ 𝑲𝑳𝒊𝒏𝒆𝒓 + 𝑲𝑮𝒓𝒐𝒖𝒏𝒅) (22) Equation 20  The liner axial hoop force 𝑭𝑳𝒊𝒏𝒆𝒓 is obtained by multiplying the liner stiffness by the supported shaft displacement.  𝑭𝑳𝒊𝒏𝒆𝒓 =∆𝑺𝒉𝒂𝒇𝒕𝑬𝑳𝒊𝒏𝒆𝒓𝑨𝑳𝒊𝒏𝒆𝒓𝒓                 (23) Equation 21  To complete the design process, the liner force must be compared to the ultimate strength of the liner material. For the case of a concrete liner, reinforcement may need to be added accordingly. FEM models were conducted in RS2 and compared to results using the proposed ring analogue equations. The liner thickness was varied from 100 mm to 400 mm in increments of 100 mm, while all other parameters were kept constant (see properties listed in Table 7). Results of axial hoop forces and shaft displacements are presented in Figure 32 and show good agreement between the proposed method and the FEM results. For the thicker liners, results begin to deviate, possibly due to additional effects that are not accounted for in the ring simplification. However, the deviation for very thick liner can still be considered as negligible, and therefore the proposed method can be used as a first approximation.   Table 7 Properties for shaft parametric study Ds [m] 6 Hs [m] 30 𝐸𝐺𝑟𝑜𝑢𝑛𝑑[MPa] 2,000 𝐸𝐿𝑖𝑛𝑒𝑟 [MPa] 20,000 𝑡𝐿𝑖𝑛𝑒𝑟[mm] 100, 200, 300, 400  69   Figure 32 Results showing comparison of (a) liner hoop axial force, and (b) displacements, vs. liner thickness using FE model and ring analogy     70   6. CHAPTER 6  APPLICATION OF THE EB METHOD FOR TUNNEL PROBLEMS 6.1. Introduction Tunneling projects are considered to be high cost civil engineering projects. Tunnel support is a major component of the costs associated with tunnel construction and can greatly vary depending on the geological conditions (Hoek et al., 2006). For the design of tunnel support there are a number of methods available, including analytical, empirical, and numerical methods. Amongst the empirical methods, the Q and RMR classification systems are commonly used (Barton et al., 1974; Bieniawski, 1989).  The Q system in particular relies on extensive databases, and provides a valuable aid for support selection. However, both systems do not account for the stresses and deformations that develop in the support system and are therefore a very limited tool in terms of actual tunnel support design. Furthermore, engineers have applied those classification systems in the context of a probabilistic/risk approach to design without considering the limitations imposed by the nature of the classification systems, which are based on scales of measurements (nominal and ordinal measurements) that are inappropriate for statistical analyses. The limitations imposed on nominal and ordinal measurements in terms of statistical analysis are generally ignored by geotechnical engineers, since the process of assigning numbers to geology produces the misconception that classification ratings are rock mass physical properties. Classification systems simply produce a qualitative ranking process for the rock mass (Elmo and Stead, 2010). The introduction of classification systems based on interval and ratio scales is therefore required if it is required to implement a statistical approach to address complex rock engineering problems. Other researchers in the literature have started questioning the use of classification systems like RQD, arguing that it should be replaced by fracture frequency, which is an interval scale measurement (for example, Pells et al., 2017). 71  Numerous analytical solutions have been developed for different conditions. Some notable solutions are mentioned in Chapter 2. Most of the analytical solutions for tunnels assume circular shaped tunnels, as the circle shape is an ideal shape mathematically. Furthermore, most tunnels are circular or near circular. The two main reasons for this are: 1) in circular tunnels the stress flow around the tunnels is ideal, unlike square excavations that attract higher stress concentrations, and 2) many tunnels are TBM driven tunnels. It is noted that TBM tunnels are currently circular, owing to the advantage of perfect rotational torque applied at the excavation face.  The term "deep tunnels" refers to tunnels where the ratio of the tunnel diameter D to the tunnel depth Z is relatively small. In this case, a constant stress field can be assumed, rendering development of closed-form solutions considerably simpler. As the ratio D/Z decreases, the effect of the increase in ground stresses with increasing depth cannot be neglected. Studies have shown that when the ratio D/Z is 0.5 or less, the tunnel can be considered as deep, for the purpose of support calculations (Ngan et al., 2017) Table 8 lists the notable analytical solutions for deep circular tunnels mentioned in Chapter 2, along with their limitations. Currently, there are no analytical methods for deep circular tunnels that provide a full solution, not limited to either one of the simplifying assumptions: non-hydrostatic stress fields, plastic ground, and full interaction with support.  Table 8 Notable analytical solutions for circular and deep tunnels and their limitations Full interaction with support Plastic ground Non-hydrostatic stress field Solution ✕ ✕ ✓ The Kirsch equations ✓ ✓ ✕ Convergence-confinement method Limited to tunnel liner ✕ ✓ Einstein and Schwartz (1979) ✕ ✓ ✓ Detournay's solution (1985)  As stated in the research objectives, the method developed in this thesis relies on available solutions for unsupported tunnels, in order to develop the double-beam analogue, and then compute the stresses and displacements of the tunnel with the addition of a support system. 72  Thus, the method provides an efficient means of tunnel analysis which includes non-hydrostatic stresses, plastic ground, and the full interaction with different types of supports and their combinations. A discussion regarding non-circular tunnels is given in Section ‎6.6. 6.2. Derivation for Elastic Ground  The following section provides a detailed outline of the derivation and computational process of the proposed method regarding the case of deep circular tunnels in an elastic medium. The MATLAB code is given in Appendix 1 so that it may be reproduced by others. The basis of the proposed method for tunnel support analysis is deriving a beam (hereinafter referred to as the Equivalent Boundary Beam, EBB) that behaves in a manner equivalent to the tunnel boundary.  For a circular tunnel under elastic conditions and a constant gravitational stress field, the displacements can be found using the Kirsch equations. The equations for the radial and tangential displacements of the excavation boundary are:   𝒖𝒓 = −𝒑𝒂𝟐𝟒𝑮𝒓[(𝟏 + 𝑲) − (𝟏 − 𝒌) {𝟒(𝟏 − 𝝂) −𝒂𝟐𝒓𝟐} 𝒄𝒐𝒔𝟐𝜽]  (24)  Equation 22         𝒖𝜽 = −𝒑𝒂𝟐𝟒𝑮𝒓[(𝟏 − 𝑲) {𝟐(𝟏 − 𝟐𝝂) +𝒂𝟐𝒓𝟐} 𝒔𝒊𝒏𝟐𝜽]                   (25)    Equation 23  where p is the initial vertical pressure, K is the coefficient of lateral pressure, G is the ground's shear modulus, a is the tunnel radius, r is the radial distance from the center, and 𝜃 is the angle measured counter clockwise from the horizontal axis.  It is noted that solutions derived for circular tunnels, such as the Kirsch solution, commonly use a 2D polar coordinate system (r- θ system) rather than a Cartesian system (x-y system). The choice between these two systems is mainly a matter of convenience. The transformation between the two systems can be achieved using available formulae. In terms of computational convenience, there are a number of advantages using the polar system for the given problem: 73   It obviates the need to transform the polar displacements obtained from the Kirsch solution into the x-y system.   For beam elements in the polar system the stiffness matrix can be reduced to a more compact 4X4 matrix (rather than the 6X6 matrix, as shown in Equation 28).  Modification of the method to include slip conditions (i.e. assuming no friction between the liner and ground, so that relative tangential movements are permitted), is more convenient in the polar system derivation, as the slip process is a tangential phenomenon.  Rock bolts in circular tunnels are usually installed in the radial direction, which coincide with the polar system.  Disadvantages of using the polar system include:  It necessitates solving the problem separately for radial and tangential displacements, and then superimposing both solutions. When dealing with plastic material, this matter is more than an issue of computational convenience, as the principle of superposition does not apply to plastic material analysis.   The derivation based on the polar system is limited to circular tunnels. Contrastingly, a formulation using the Cartesian system enables simple modification of the method for more generalized tunnel geometries (i.e. non-circular excavations),   For these reasons, the derivation described in this Chapter assumes a Cartesian system.  For an infinite cylinder shaped beam, subjected to a uniform radial pressure p, the elastic solution for the radial displacement is: 𝒖𝒓 = −𝒑𝒓𝟐𝑬𝒕   (26) Equation 24  where E is the Young's modulus of the beam material , and t is the thickness of the beam.  The objective is to find the equivalent thickness of the beam which yields the identical displacement as a circular tunnel. By assigning K=1 and a=r into Equation 23 (i.e. uniform pressure condition), and equating Equations 24 and 26, the thickness of the equivalent beam 𝒕𝒆𝒒 is found to be: 74   𝒕𝒆𝒒 = 𝒓 (𝟏 + 𝝂)⁄   (27)  Equation 25   Figure 33 Conversion of the plane strain problem to a double beam problem  Using 𝒕𝒆𝒒, the ground can be represented by an equivalent beam, and the plane strain problem can be converted into a simpler and more intuitive double cylinder problem, as illustrated in Figure 33. Under hydrostatic loading conditions, the initial in-situ stress will be distributed between the ground and liner according to the ratio of each component's axial stiffness (equal to the product of the Young’s modulus and cross section area), to the total stiffness.  In order to address the more complex and realistic condition of non-hydrostatic in-situ stresses, the fictitious forces that act upon the tunnel-support boundary must be found. Note that these forces are fictitious, as they are not identical to the initial forces acting upon the excavation in the plane strain problem, as these two structural entities respond differently to equal loading.  75  The two main computational steps (see Figure 34 and Figure 35) used to solve the problem are:  1. The unsupported tunnel displacements are initially used to obtain the vector of fictitious forces; and 2. The fictitious forces are than applied upon the combined stiffness matrix of the ground and support to find the displacements of the supported tunnel. Once the displacements of the supported tunnel are found, the inner stresses in the support components can be deduced using basic mechanical theory.   Figure 34 Summary of the two main EBB computation steps. Kg represents the global stiffness matrix of the equivalent beam, Dus is the unsupported tunnel displacement vector, Ffic is the vector of fictitious forces, Ks ad Kst are the stiffness matrices of the support and supported tunnel, and Ds is the supported tunnel displacement vector.  76  a)  b)  Figure 35 Boundary beam elements. 𝑭𝒏 and 𝑴𝒏are the fictitious nodal forces and moments acting on: a) the equivalent boundary beam at the first stage, and b) the lined tunnel at the second stage.  77  It is noted that the proposed method would apply to both concrete liners applies and steel sets, as the inputs would be the same. For composite supports such as steel sets embedded in shotcrete or concrete reinforced with rebar, an additional analysis step must be applied for finding the distribution of stresses between the concrete and steel components. This technique was proposed by (Carranza-Torres and Diederichs, 2009) and called the “equivalent section approach."  There are various methods for solving beam problems. For this thesis, the Direct Stiffness Method (DSM) is used. The DSM is a well-known numerical method that makes use of the members' stiffness relations for computing member forces and displacements in structures (Felippa, 2001). DSM is implemented in different commercial software for structural analysis. The procedure for deriving DSM will not be reviewed here, but the basic steps required to formulate the computational process discussed in this chapter will be briefly outlined. The elements used for assembly of the boundary beam stiffness matrix are beam elements. As explained in Section ‎3.2, it is both more accurate and convenient and to convert the ground into structural elements identical to the support system. The tunnel liner behaves as a beam (i.e. its response includes axial forces and bending moments), therefore the ground is converted into a beam as well. The integration of bolts will be discussed in Section ‎6.4.   Beam elements consist of a line with two nodes, where the degrees of freedom for each node are for translational, axial and rotational movement (see Figure 36). Thus, there is a total of 6 degrees of freedom for each beam element.   Figure 36 Degrees of freedom for a beam element 78  Accordingly, the beam stiffness matrix is a 6-by-6 matrix, and is: 𝑲𝑩𝑬 =[      𝑨𝑬/𝑳 𝟎 𝟎 −𝑨𝑬/𝑳 𝟎 𝟎𝟎 𝟏𝟐𝑬𝑰/𝒍𝟑 𝟔𝑬𝑰/𝒍𝟐 𝟎 −𝟏𝟐𝑬𝑰/𝒍𝟑 𝟔𝑬𝑰/𝒍𝟐𝟎 𝟔𝑬𝑰/𝒍𝟐 𝟒𝑬𝑰/𝒍 𝟎 −𝟔𝑬𝑰/𝒍𝟐 𝟐𝑬𝑰/𝒍−𝑨𝑬/𝑳 𝟎 𝟎 𝑨𝑬/𝑳 𝟎 𝟎𝟎 −𝟏𝟐𝑬𝑰/𝒍𝟑 −𝟔𝒍𝑬𝑰/𝒍𝟐 𝟎 −𝟏𝟐𝑬𝑰/𝒍𝟑 −𝟔𝑬𝑰/𝒍𝟐𝟎 𝟔𝑬𝑰/𝒍𝟐 𝟐𝑬𝑰/𝒍 𝟎 −𝟔𝑬𝑰/𝒍𝟐 𝟒𝑬𝑰/𝒍 ]       (28) Equation 26  where l is the beam element length,  E is the Young's modulus, I is the inertial moment, and A is the cross section of the beam element.  The derivation of the beam stiffness matrix can be found in many standard text books, and will not be repeated here. In general, expressions are derived so that a single component kij in the stiffness matrix corresponds to the force reaction in the ith degree-of-freedom, to a unit displacement of the jth degree-of-freedom. As the objective is to formulate an equivalent circular beam comprised of many small beams. It is clear that increasing the number of beam elements causes the series of beams to become closer to an ideal circle.  The single beam element stiffness matrix corresponds to a local axis. Since each beam element has a different orientation, it is necessary to transform each beam element matrix from the local coordinate system to a global coordinate system. This is achieved by using a local-global transformation matrix, L, equal to: 𝑳 =[      𝑳𝒙 𝑳𝒚 𝟎 𝟎 𝟎 𝟎−𝑳𝒚 𝑳𝒙 𝟎 𝟎 𝟎 𝟎𝟎 𝟎 𝟏 𝟎 𝟎 𝟎𝟎 𝟎 𝟎 𝑳𝒙 𝑳𝒚 𝟎𝟎 𝟎 𝟎 −𝑳𝒚 𝑳𝒙 𝟎𝟎 𝟎 𝟎 𝟎 𝟎 𝟏]        (29)  Equation 27  where Lx and Ly are the vector components of each beam element equal to cosθi and sinθi, where θi is the angle of the ith  beam element relative to the global horizontal axis. Each global beam element stiffness matrix is equal to: 𝑲𝒈𝒍𝒐𝒃𝒂𝒍 = 𝑳𝑻𝑲𝒍𝒐𝒄𝒂𝒍′𝑳  (30) Equation 28  79  Once each beam element stiffness matrix is computed according to its angle θi, the global structural stiffness matrix needs to be assembled.  The process of assembly requires that element matrices components kij  that share a global degree-of-freedom with the neighboring elements are summed. Equation 31 shows an example of the assembly of two neighboring beam elements. For each matrix component the superscript denotes the element number, and the two-digit subscripts denote the row and column of each component corresponding to the single element matrix given in Equation 28. The stiffness matrix beam of beam elements #1 and #2 are marked in blue and red respectively, and correspond to the example shown in Figure 37. It can be seen that the overlap in the global matrix corresponds to the shared degrees-of-freedom 4, 5 and 6, in this example. (31)  Equation 29  80   Figure 37 Example of two beam elements   Clearly, solving such matrices manually for more than 2-3 elements would become a highly cumbersome task. Nonetheless, this method is ideal for computer implementation.  For each structural problem, it is essential to assign proper boundary constraints. As the problem is symmetric in both axes, it is possible to treat only one quarter of the tunnel. Constraints are defined at the tunnel crown and wall accordingly, as shown in Figure 38: at the tunnel wall, or spring-line, owing to symmetry and continuity, translation is permitted in the horizontal direction, and vertical translation and rotation are prohibited. Similarly, at the tunnel crown, vertical translation is permitted, and horizontal translation and rotation are prohibited.  81   Figure 38 Application of symmetry and constraints  For the element and formulation process, first, the number of elements, n, is defined, and the angle θ is defined as 90o/ n. The ith coordinate (xi, yi) is then computed as: 𝒙𝒊 = 𝒓 ∗ 𝐜𝐨𝐬 ((𝒊 − 𝟏) ∗ 𝜽)   (32) Equation 30  𝒚𝒊 = 𝒓 ∗ 𝐬𝐢𝐧 ((𝒊 − 𝟏) ∗ 𝜽)   (33) Equation 31  where i is the coordinate number (i=1,…,n). Figure 39 shows an example for a formulation of two elements. 82   Figure 39 Example of coordinate and element formulation for the assembly of two equivalent beam elements  In numerical methods, results must converge from a certain number of elements. Analyses to obtain result convergence vs. the number of elements were conducted. Figure 40 shows a plot of result convergence vs. number of elements, and it can be seen that from a number of elements in the range of 50, results deviate from the final convergence result by only 2-3 percent. As a comparison, 2D plane strain FEM models of tunnels typically consist of thousands of elements and nodes. 83   Figure 40 Result convergence vs. number of boundary-beam elements  In addition to the stiffness matrix assembly procedure, the vector of nodal displacements must be constructed. As discussed, the Kirsch equations are used as for computing the displacements for elastic ground conditions. The displacement vector for n elements consists of a 3*(n+1) rows, so that corresponding to each coordinate a horizontal, vertical and rotational component (denoted respectively as: x, y, and θ) are assigned: 𝑫 =[      𝒖𝟏−𝒙𝒖𝟏−𝒚𝒖𝟏−𝜽⋮𝒖𝒏−𝒙𝒖𝒏−𝒚𝒖𝒏−𝜽]         (34) Equation 32   84  The nodal rotations 𝑢𝑖−𝜃 can be obtained by differentiation of the radial displacements given in Equation 23: 𝒖𝒊−𝜽 =𝒅𝒖𝒓𝒅𝜽= −𝒑𝒓𝟐𝑮[(𝟏 + 𝑲) − (𝟏 − 𝒌){𝟒(𝟏 − 𝝂) − 𝟏}𝒔𝒊𝒏𝟐𝜽]  (35) Equation 33  The procedure for the assembly of the global stiffness matrix of the tunnel liner is identical to that of the equivalent boundary beam, only that the inputs are the actual input parameters of the tunnel liner. These parameters include the Young's modulus, thickness, and Poisson ratio of the concrete, or shotcrete, used for the liner. Finally, the two computational steps shown in Figure 34 can be carried out. The first step is to find the fictitious forces by multiplying the stiffness matrix of the equivalent boundary beam by the unsupported displacements. These forces are fictitious, as they are not identical to the initial forces acting upon the excavation in the plane strain problem, as these two structural entities respond differently to equal loading.  Figure 41(a) shows an example of a plot of the displacements calculated using the methodology outlined in this chapter. As shown in the legend, along with the original tunnel contour, the unsupported and supported tunnel deformation contours are plotted. These deformations contours have been magnified by a factor of 100 so that they are noticeable. As mentioned, the unsupported displacements are given by the Kirsch equations. The supported tunnel deformation contour lies between the original tunnel contour and the unsupported displacement contour.  Another such example is shown in Figure 41(b), where the K ratio has been reduced relative to Figure 41(a). It can be seen there that the outward displacement of the tunnel at the horizontal axis is greater with support than it is without support. 85  a)  b)  Figure 41 Examples of plots of displacement contours generated through MATLAB and magnified by a factor of 100 using the proposed methodology, where a) the K ratio is smaller than 1, and b) the K ratio is larger than 1.  86   Once the supported displacements are found, the internal forces of the support system can be deduced using standard procedures based on beam theory. The change in length in each beam element dictates its inner axial force, according to the relationship based on Hooke's law. Hence, for a single beam element, the axial force is equal to: 𝑵 = ∆𝑫𝒄𝑳   (36) Equation 34  where N is the axial force, ∆ is the difference in element length, and 𝐷𝑐 is the compressibility coefficient. Under plain strain conditions 𝐷𝑐 is equal to:  𝑫𝒄 =𝑬𝑨𝟏−𝝂   (37) Equation 35  where E is the Young's modulus of the liner material, A is the cross section area of the liner element, and 𝜈 is the Poisson ratio of the liner material. In addition to axial forces which occur at the element level, when the angle between two elements is altered, bending moments are introduced. This change in angle (or slope) can be found by differentiating the expression of the beam deflection. The relationship between the bending moment M and the beam deflection is: 𝑴 = 𝑲𝒄𝒅𝟐𝒘𝒅𝒙𝟐   (38) Equation 36  where w is the beam transverse deflection, x is the incremental beam length, and Kc is the flexibility coefficient, which under plain strain conditions is equal to: 𝑲𝒄 =𝑬𝑰𝟏−𝝂   (39) Equation 37  where E is the Young's modulus of the liner material,  I is the inertial moment of the liner element, and 𝜈 is the Poisson ratio of the liner material. Based on Equation 35, the second derivative of the function of the tunnel displacements can be used to find the liner bending moments. It is noted that for the computations undertaken in 87  the thesis, the simpler Bernoulli-Euler beam theory was used. As opposed to the more generalized Timoshenko beam formulation, the Bernoulli-Euler formulation neglects deformations due to shear forces. This assumption is generally valid when the length of the beam is significantly larger than its thickness, as is the case for practical applications involving tunnel liners.  The liner maximum and minimum inner stresses can be calculated by superimposing the stresses due to axial forces and bending moments; accordingly: 𝝈𝒎𝒊𝒏𝒎𝒂𝒙 =𝑵𝑨±𝑴𝑰𝒕 𝟐⁄    (40) Equation 38  where N is the axial force,  A is the liner cross section area, M is the bending moment, I is the moment of inertia of the liner section, and t is the liner thickness.  As shown in Figure 35, the stiffness matrix of the support tunnel Kst is computed by summation of the stiffness matrices of the two beams: the equivalent beam Kg which represents the ground, and the beam that represents the support Ks  (Kg+ Ks= Kst). In order to perform this summation, it must be assumed that the neutral axis of both beams coincide with each other. For most cases, the neutral axis of a single beam lies at the geometric center of its cross section.  Figure 42(a) shows an example of two cross sections which share an identical neutral axis, and Figure 42(b) shows an example of two attached beams which have separate neutral axes. Obviously, in reality, the tunnel liner is attached below (or inwards) of the ground, as is the case in Figure 42(b). The radius of the equivalent beam is the external radius re, and the radius of the support is the internal radius, ri. This difference in radii must be accounted for when computing the final liner stresses.  Hence, initial results of axial forces must be corrected by multiplying by a correction factor: 𝒄. 𝒇. = 𝒓𝒆 (𝒓𝒊⁄ − 𝟎. 𝟓 ∗ 𝒉𝒍)   (41) Equation 39  where c.f. is the correction factor, and ℎ𝑙  is the liner thickness.  88  For computing the liner bending moments no additional adjustment is necessary, as both beams (i.e. the equivalent ground beam and tunnel liner) follow the same curvature.  Figure 42 a) Example of two beam cross sections which share the same neutral axis.        b) Example of two beams with different neutral axes.  Finally, the factor of safety FoS of the support system can be calculated according to the ratio of the ultimate strength of the liner material 𝜎𝑢𝑙𝑡 , to the maximum resultant stress, 𝜎𝑟𝑒𝑠: 𝑭𝒐𝑺 =𝝈𝒖𝒍𝒕𝝈𝒓𝒆𝒔   (42) Equation 40  In the context of finding the liner FoS, it is important to examine both ends (i.e. top and bottom fibers) of the tunnel liner. Under axial pressure, the whole cross section is subjected to a uniform stress. Under bending, the stress varies along the cross section from maximum compression at one end, through zero stresses at the neutral axis, and maximum tension at the other end. Figure 43 shows illustrations of the superposition of axial and bending stresses. It can be seen, that according to the values of the axial and bending stresses, the cross section of 89  the beam can either be entirely in compression, or simultaneously in compression and tension at the top and bottom fibers. In the latter case, Equation 40 must be checked for both the ultimate compressive and tensile ultimate stresses, and the FoS is determined according to the lowest of the two results.  As concrete is generally weak in tension, additional calculations can be used to determine the extent of steel reinforcement that can sustain the tensile force. Such calculations are part of standard practice in civil engineering, and appear in textbooks and local building codes.     Figure 43 Examples of the superposition of axial and bending stresses, when a) the entire cross section is in compression, and b) the bottom fiber is in tension.   To summarize what has been discussed thus far, based on the problem input parameters, the EB method can be implemented to obtain the supported tunnel displacements. These results can then be used to compute liner stresses based on known analytical relationships between stress, strain, and problem geometry.  90  In order to properly assess overall tunnel stability, it is essential to gain knowledge of the maximum stresses that develop in the ground, as well as the liner. By entering the parameters of the equivalent boundary beam into Equation 34 and 37, the axial forces and bending moments that develop in the equivalent boundary beam can be found. However, this beam is a fictitious entity, and therefore these results have no practical use, i.e. they are not equal to the actual stresses that develop in the ground. Apparently, there is no straightforward method for re-converting the resultant stresses of the equivalent boundary beam to the stresses that develop in the original plane-strain space. It is possible that a mathematical formulation can be derived for this specific purpose. Alternatively, a simple procedure for obtaining approximate results of the rock stress is hereby proposed. This procedure includes the following steps: 1. The maximum tangential stress of the unsupported tunnel is calculated using the equation given by the Kirsch solution: 𝝈𝜽𝜽 =𝒑𝟐[𝟐(𝟏 + 𝒌) + 𝟒(𝟏 − 𝒌)𝒄𝒐𝒔𝟐𝜽]   (43) Equation 41  2. The displacement of the unsupported tunnel at a given reference point (e.g. the tunnel crown or wall) is calculated using Equation 24 given by Kirsch. 3. The displacement of the supported tunnel at the same reference point is given by the proposed EB method, as outlined in this chapter. 4. The ratio of supported tunnel displacement (from Step 3) over the unsupported tunnel displacement (from Step 2) is computed and multiplied by the maximum stress (from Step 1). The summary of steps 1-4 in equation form gives: 𝝈𝟏−𝒈𝒓𝒐𝒖𝒏𝒅 ≅ 𝝈𝜽𝜽−𝒖𝒏𝒔𝒖𝒑𝒑𝒐𝒓𝒕𝒆𝒅 ∗𝒖𝒊−𝒔𝒖𝒑𝒑𝒐𝒓𝒕𝒆𝒅𝒖𝒊−𝒖𝒏𝒔𝒖𝒑𝒑𝒐𝒓𝒕𝒆𝒅 (44) Equation 42  where 𝜎1−𝑔𝑟𝑜𝑢𝑛𝑑 is the major stress that develops in the ground in the supported tunnel, 𝜎𝜃𝜃 is the tangential stress in the unsupported tunnel (from Equation 43),  𝑢𝑖−𝑠𝑢𝑝𝑝𝑜𝑟𝑡𝑒𝑑 is the displacement of the supported tunnel at a specific reference point, and 𝑢𝑖−𝑢𝑛𝑠𝑢𝑝𝑝𝑜𝑟𝑡𝑒𝑑 is the displacement of the unsupported tunnel at the same reference point. 91  In order to examine the degree of accuracy of given by Equation 44, an example was investigated. Table 9 lists results of the maximum major principal stresses that develop in the ground, at the ground-support interface. The results refer to a tunnel with a radius of 2.5 m, Young’s moduli of 10 GPa for the ground and 30 GPa for the liner, a vertical pressure of 10 MPa, a K ratio of 0.5, and  liner thickness values that vary from 100 to 600 mm. It can be seen that the results of the major principal stresses obtained using Equation 44 are in agreement with the results from FEM models.  Table 9 Comparison of results of maximum major principal stresses for lined tunnels 𝜎3−𝑔𝑟𝑜𝑢𝑛𝑑 [MPa] 𝜎1−𝑔𝑟𝑜𝑢𝑛𝑑 [MPa] Liner thickness [mm] Error [%] RS2 EB (Equation 44) 2.9 3.08 22 22.7 100 3.6 4.29 20.1 21 200 4.2 2.04 19.2 19.6 300 4.6 1.09 18.2 18.4 400 4.9 -2.86 18 17.5 500 5.9 -3.61 17.2 16.6 600  It is important to note that the stresses in the ground at the ground-support interface consist of radial stresses (or minor principal stresses) as well as tangential stresses (or major principal stresses). As opposed to the unsupported tunnel, where the ground stresses at the excavation boundary are purely tangential, the introduction of the tunnel lining introduces additional radial stresses, imposed by the liner upon the ground. In Table 9, the minor principal stresses are listed alongside the major principal stresses. Apparently, as the minor principal stresses increase, the deviation between results increases as well.  For typical conditions, the results obtained via Equation 44 will provide a conservative upper bound value. This is because the minor principal stresses add confinement and increase the ultimate strength of the ground, when assuming Mohr-Coulomb behavior. This increase in strength is a function of the ground friction angle 𝜃 and the unconfined compressive strength UCS of the geological material, as given in Equation 45: 𝝈𝟏 = 𝝈𝑼𝑪𝑺 +(𝟏+𝒔𝒊𝒏𝜽)(𝟏−𝒔𝒊𝒏𝜽)𝝈𝟑   (45) 92  Equation 43  Finally, after determining the maximum principal stresses in both the ground and liner, engineering judgment is required to assess whether the FoS can determined based solely on the results of the elastic analysis, or if additional analyses are required. Generally, if the maximum principal stress that develops in the ground is below the ground strength, the tunnel FoS can be calculated with Equation 42, and according to the stress and strength of the support system alone. The complex nature of geological materials and environments may require additional engineering judgment in order that the ground strength is properly assessed. For example, in tunnels in brittle rock under high in-situ stresses, rock failure may be triggered well below the UCS. For a discussion and practical example regarding the application of the proposed EB method to tunnels in brittle rock, see Chapter ‎7. For plastic and deformable rock, if preliminary elastic analysis shows that the rock exceeds its elastic limit, additional analysis that accounts for plastic behavior must be carried out. For application of the proposed EB method for elastoplastic analysis, see Section ‎6.5; there, a methodology is derived so that the tunnel-support interaction analysis accounts for the non-linear behavior of plastic rock.  Figure 44 shows a flowchart that briefly summarizes the discussion given here regarding the integration of EB elastic analysis into the process of tunnel design.  93   Figure 44 Summary of design process with EBB elastic analysis  Another fundamental aspect pertinent to tunnel-support interaction analysis is the effect of the timing of support installation. As discussed in the literature review (Chapter ‎2), the Longitudinal Displacement Profile (LDP) is a term used for describing the deformation curve from the tunnel face to the point of support installation. Knowledge of the LDP can assist in assessment of the magnitude of initial displacements induced prior to support installation. In order to simulate the delay in support installation in a 2D model it is possible to stage the project and split the loads according to the anticipated initial deformation. Given the elastic condition assumption, the procedure for load split implementation in the proposed EB method 94  is simple. The computation process is divided into two steps: 1) the pre-support installation stage, and 2) the support installation stage. The proper portion of the total stress is assessed according to the LDP, and assigned to each stage accordingly. For example, if the initial assessment leads to the conclusion that 40% of the total unsupported tunnel deformation will occur prior to support installation, then the computation will be as the following: 1. The deformations of the unsupported tunnel are calculated by multiplying the initial vertical stress by a factor of 0.4. 2. The final stresses acting on the support system are calculated using a factor of 0.6 for the initial vertical stress. These stresses are the final stresses. 3. The final displacements are the sum of the displacements given by the first and second steps.  6.3. Validation for Tunnel Liner Analysis  Comparison of results obtained from 2D numerical analyses conducted via RS2 and the proposed EBB method are presented in the following tables.  The models generated in RS2 consisted of 6-noded triangle elements, with 6,000-8,000 nodes. Parameters were varied, not including the vertical pressure and Poisson ratio that were kept constant in all models: a vertical pressure of 10 MPa, and Poisson ratios of 0.25 and 0.2 for the rock and concrete liner, respectively. All models consist of a single stage, with no application of initial relaxation prior to support installation. Table 10 lists shows the inputs and results for a lined tunnel. In this table, a comparison of results of the maximum and minimum displacements, which occur at the tunnel crown and wall, is presented. In addition to results from RS2 FE models, solutions obtained through Einstein and Schwartz's analytical solution (1979) were calculated and listed as well. As can be seen, results using the proposed EBB method are generally closer to the ideal analytical solution than results obtained through FEM modeling.  Minor deviations between different methods are inevitable, for a number of reasons:  95  1. FEM codes do not provide solutions in exact agreement with the analytical solution and deviations in the range of 5-10% are not uncommon (see RS2 verification).  2. Minor deviations have been found even between results using different analytical solutions (Hung et al., 2009). 3. The proposed EBB method itself includes some form of inaccuracy due to the division of the beams into discrete elements, and possibly also due to the reliance on a fictitious equivalent entity.  Nevertheless, given the observable heterogeneous nature of physical parameters and/or geological processes, these deviations could be considered as negligible.    96  Table 10 Result comparison for lined tunnels. E&S refers to results obtained via the analytical solution given by Einstein and Schwartz (1979). Crown displacement [mm] Wall displacement [mm] Input parameters  EBB E&S RS2 EBB E&S RS2 K ratio Tunnel radius [m] Liner thickness [cm] E Liner [GPa] E Rock [GPa] Model # 2.8 2.7 3 0.4 0.4 0.45 0,5 2 10 30 10 1 1.46 1.45 1 0.25 0.25 0.24 0,5 2 10 30 20 2 1 1 1 0.18 0.18 0.17 0,5 2 10 30 30 3 0.75 0.75 0,7 0.14 0.14 0.12 0,5 2 10 30 40 4 2.67 2.65 3 0.33 0.34 0.36 0,5 2 10 40 10 5 2.58 2.55 2.1 0.28 0.28 0.26 0,5 2 10 50 10 6 2.5 2.5 2.1 0.23 0.22 0.19 0,5 2 10 60 10 7 2.5 2.5 2.1 0.23 0.23 0.2 0,5 2 20 30 10 8 2.3 2.3 2.1 0.1 0.1 0.1 0,5 2 30 30 10 9 2.13 2.1 2.1 -0.01 -0.01 -0.01 0,5 2 40 30 10 10 5.9 5.8 6 1 1 1 0,5 4 10 30 10 11 9 8.9 9 1.6 1.6 1.7 0,5 6 10 30 10 12 12 12 10 2.2 2.3 2.2 0,5 8 10 30 10 13 2.1 2.1 2 2.1 2.1 2 1 2 10 30 10 14 0.8 0.8 0.7 5.5 5.5 5 2 2 10 30 10 15 -0.5 -0.5 -0.4 8.9 8.9 8,6 3 2 10 30 10 16  The results listed in Table 10 are maximum results, which for circular tunnels occur at the points of either the tunnel wall and/or crown (i.e. θ=0 and 90o according to the definition shown in Figure 45(a)). Figure 45 shows a comparison of the distribution of liner axial forces and bending moment along the tunnel liner, from wall to crown, for the problem parameters listed in Table 11. With respect to the bending moments, it can be seen that the results EBB are smoother than those obtained by the FE model.   97  Table 11 Parameters for model comparison Vertical pressure [MPa] 10 K ratio 3 Rock Young’s Modulus [GPa] 40 Concrete liner Young’s Modulus [GPa] 60 Ground Poisson ratio 0.25 Concrete Poisson ratio 0.2 Liner thickness [cm] 20   a)  98  b)  c)  Figure 45 a) RS2 numerical model, showing the angle measured from horizontal (wall) to vertical (crown) b) Result comparison between RS2 and EBB of axial forces (using the compression-negative sign convention)  c) Result comparison between RS2 and EBB of bending moments 99   6.4. Rock Bolt Analysis Rock bolts are a common type of support used in underground mines and civil tunnels. Rock bolts contribute to the resistance to inward tunnel deformations. This action of the bolts is studied using analytical and numerical solutions for tunnel support interaction analysis. There are additional roles and benefits to rock bolts, for example stabilization of loose rock wedges (Hoek, www.rocscience.com/education/hoeks_corner).  However, these additional roles are generally not accounted for within regular stress and strain analyses, and are treated separately, commonly using kinematic analysis for the detection and stabilization of loose wedges. Generally, support can be classified as being either passive or active (Brady and Brown, 2006). Active support imposes a predetermined load to the rock surface at the time of installation, and passive support gains pressure as the rock mass deforms. Bolts can act as either passive or active support. Pre-tensioned rock bolts apply to the first category, while un-tensioned grouted bolts or cables are considered passive support.  There are different types of passive bolts, with the most common used in practice being:   Grouted bolts rely on the adhesive bonding created by the injection of grout, or resin, into the void between the rock borehole and the bolt. The interaction between the rock and bolt occurs along the full length of the grouted bolt. Hence, it is necessary to derive continuous or discretized solution procedures to account for the variation in stresses and strains along the bolt (see Figure 46(a)).   Mechanically anchored bolts rely on an expansion shell which expands against the wall of the borehole. Thus, the interaction occurs only at two points: 1) the face of the bolt, and 2) the expansion shell. Hence, solution procedures for mechanically anchored bolts can rely on a single force or element to account for the rock-bolt interaction (see Figure 46(b)). 100  a)  b)  Figure 46 Illustration of a) grouted bolts, and b) mechanically anchored bolts  The incorporation of the action of passive grouted and ungrouted bolts, as well as active bolts, into the proposed EB method for tunnel analysis, will be addressed in the following sections. The following sections assume elastic behavior of the ground, as the formation of a plastic radius around the tunnels further complicates the interaction with the bolts. 6.4.1. Pre-tensioned Bolts The procedure for integrating pre-tensioning forces into the proposed method, can be achieved by simply adding the pre-tensioning forces into the vector of fictitious forces (found from Step 1 in Figure 34), and then applying this modified force vector upon the supported tunnel (Step 2 in Figure 34). The bolt pre-tensioning forces resist the gravitational loads and therefore they must be subtracted from the force vector components. The bolt force vector must be assembled so that the forces correspond to the degree-of-freedom of the bolts' 101  locations. Hence, for a bolt at the kth degree-of-freedom with a pre-tensioning force 𝑓𝑏𝑜𝑙𝑡, the fictitious force vector 𝐹𝑚 is adjusted accordingly: [𝑭𝒎] =[     𝑭𝟏𝑭𝟐…𝑭𝒌…𝑭𝒏]     +[     𝟎𝟎…𝒇𝒃𝒐𝒍𝒕…𝟎 ]       (46) Equation 44   6.4.2. Mechanically Anchored Bolts In the case of passive support in the form of ungrouted mechanically anchored bolts, the bolts are assumed to resist the ground deformation and to share the load imposed by gravity with the surrounding ground. This action of the bolts can be integrated into the methodology discussed previously in this chapter by adding spring elements to the equivalent beam representing the ground. Based on elastic theory, the stiffness of a bolt, without accounting for the interaction with the ground, is equal to: 𝑲𝒃𝒐𝒍𝒕 =𝑨𝒃𝒐𝒍𝒕𝑬𝒃𝒐𝒍𝒕𝑳𝒃𝒐𝒍𝒕  (47) Equation 45  where 𝐴𝑏𝑜𝑙𝑡 is the bolt cross section area, 𝐸𝑏𝑜𝑙𝑡 is the Young's modulus of the bolt , and 𝐿𝑏𝑜𝑙𝑡 is the length of the bolt.  However, the actual stiffness of the bolt is less than the expression given in Equation 47 as some movement occurs at the point of the bolt end. To account for this effect, the radial displacement given by the Kirsch solution (see Equation 24) can be used to derive a reduction factor. This reduction factor is calculated separately for each bolt according to the bolt location defined by the angle 𝜃 measured from the horizontal axis. The reduction factor rf  is equal to: 𝒓𝒇 = 𝟏 −𝒖𝒂𝒏𝒄𝒉𝒐𝒓𝒖𝒆𝒙𝒄𝒂𝒗𝒂𝒕𝒊𝒐𝒏 (48) Equation 46  102  where 𝑢𝑎𝑛𝑐ℎ𝑜𝑟 and 𝑢𝑒𝑥𝑐𝑎𝑣𝑎𝑡𝑖𝑜𝑛 are the displacements of the unsupported tunnel calculated at the points of the bolt anchor and excavation boundary, respectively. The values of 𝑢𝑎𝑛𝑐ℎ𝑜𝑟 and 𝑢𝑒𝑥𝑐𝑎𝑣𝑎𝑡𝑖𝑜𝑛 can be obtained by assigning  𝑟 = 𝑎 + 𝐿𝑏𝑜𝑙𝑡  and r=a into Equation 24, where a is the tunnel radius, and 𝐿𝑏𝑜𝑙𝑡 is the bolt length.  A comparison for varied parameters of a tunnel supported by both liner and bolts is presented in Table 12, along with the constant parameters listed in Table 13.   Table 12 Result comparison for bolted tunnels Max. Displacement [mm] Max. Bolt axial force [MN] Min. Bolt axial force [KN]  EBB RS2 EBB RS2 EBB RS2 Pre-tension force [MN] K ratio Liner thickness [mm] Bolt length [m] 4.3 4.2 64 67 21 23 0 0.5 200 2 4.3 4.2 51 52 12 14 0 0.5 200 4 4 4 47 49 8.9 11 0 0.5 300 4 3.8 3.8 45 46 6 8 0 0.5 400 4 8.5 8.4 101 103 25 27 0 2 200 4 13.5 14 156 163 9 11 0 3 200 4 4.1 4.1 NA NA NA NA 0.5 0.5 200 4  Table 13 Liner and bolt parameters used for combined liner-bolt support comparison  Bolt spacing [m/m] Bolt diameter [mm] E bolt [GPa] E liner [GPa] E rock [GPa] Tunnel radius [m] 1/1 25 200 30 10 3.2    6.4.3. Grouted Bolts In the case of passive support in the form of grouted bolts, the bolts are assumed to resist the ground deformation and share the load imposed by gravity with the surrounding ground. Constructing confinement curves that represent grouted bolts is associated with some 103  theoretical difficulties. The main difficulty is that the underlying assumption of the convergence-confinement method is to treat only the tunnel boundary, whereas the interaction of the bolts with the rock varies in a non-linear fashion in the radial direction. This is different and more complicated than the case of the mechanically anchored bolts, which can be treated as single elements, due to their simpler boundary conditions that dictate a constant force response. Attempts have been made to simplify bolt behavior by assigning equivalent support pressures from within the tunnel wall (Bischoff and Smart, 1975). Grasso et al. (1991) introduced the concept of a reinforced rock mass and modeled the bolt as a fictitious increase in rock cohesion. Oreste and Peila (1996) developed a numerical method that accounts for additional factors such as the distance of the bolting pattern to the tunnel face. Nie et al. (2014) developed rockbolt models within the frame work of the two-dimensional discontinuous deformation analysis (DDA) by treating the rockbolt behavior as elastic, linear, strain-hardening. The references mentioned all assume the simplifying assumption of a hydrostatic stress field. Under a non-hydrostatic stress field, the rock-bolt interaction is further complicated as each bolt is subjected to different stresses and displacements.  The procedure proposed in this thesis is to treat the passive bolts as a series of unequal springs. The spring stiffness is calculated for each bolt at a separate preliminary stage. In this preliminary stage, the calculation process is similar to the concepts outlined in the previous sections. Fictitious “rock bar” elements are assumed, in order to represent the rock surrounding the bolts (see Figure 47). In reality, the behavior of the rock surrounding the bolts is different and more complex than bars, which consist of pure axial behavior. However, bars best represent the action of the bolts; therefore, the surrounding ground is converted into equivalent bars.  The calculations in the preliminary stage include the following steps: 1. The rock surrounding each bolt is assumed to behave as a separate "rock bar". 104  2. The bolt parameters must be defined: bolt diameter, length, Young's modulus, and location. The location is defined as a function of the bolt’s angle from the horizontal axis. 3. Both the bolt and "rock bar" are discretized into small elements, each consisting of a length of dL (see Figure 47). Each "rock bar" element is assigned a bar element stiffness of 𝐸𝑅 ∗ 𝐴𝐹−𝑅𝐵 𝑑𝐿⁄ , where 𝐸𝑅 is the rock’s Young’s modulus, and 𝐴𝐹−𝑅𝐵 is an equivalent cross section of the “rock bars.”  4. Global stiffness matrices for the rock bars and bolts (Krock-bars and Kbolts) are assembled according to the parameters defined in the previous steps. The stiffness matrix of the combined “rock bar” and bolt elements is calculated by simple summation of the two matrices. 5. The “rock bar” displacements are computed for each bolt according to the radial displacements of the unsupported tunnel given by the Kirsch equations (Equation 24). Using the relationship between force, stiffness and displacement (F=K*D), these displacements, along with the "rock bar" stiffness matrix, are used for finding the vector of fictitious forces. 6. The fictitious forces found in Step 5 and the stiffness matrix of the combined rock and bolt system found in Step 4 are used to compute the final displacements along the bolts. This is achieved by using the relationship between force, stiffness and displacement, this time accounting for the stiffness of the combined system: 𝑫𝒇𝒊𝒏𝒂𝒍  = 𝑭𝒇𝒊𝒄𝒕𝒊𝒕𝒊𝒐𝒖𝒔/(𝑲𝒓𝒐𝒄𝒌−𝒃𝒂𝒓𝒔 + 𝑲𝒃𝒐𝒍𝒕𝒔)  (49)  Equation 47  7. The inner stresses acting along each bolt can be found by multiplying the displacements found in the previous step, by the bolt stiffness. The above steps refer to a tunnel supported solely by bolts. For tunnels supported by a combination of liner and rockbolts, the following additional steps must be implemented: 8. The force acting on each bolt is divided by the final bolt displacement at the rock-bolt interface so that the equivalent spring stiffness of each bolt is found.  9. The spring elements with the stiffness found in the Step 9 are assigned to the combined rock and liner model. 105  10. The EBB analysis is computed according to the procedure presented in Section ‎6.2.  Using this methodology within the EBB method, the final solution for a tunnel support system consisting of both liner and bolts captures the effects due to the interaction between the different components, rock, bolt and liner: the bolts cause a reduction in the maximum stresses and displacements acting on the liner; in addition, the bolts cause the liner to develop localized bending moments around the area of the bolts. The equivalent cross section 𝐴𝐹−𝑅𝐵 of the "rock bar" must first be found in a preliminary stage through a trial and error process. Accordingly, this cross section is found to be: 𝑨𝑭−𝑹𝑩 ≅ 𝟎. 𝟓 𝒎𝟐   (50) Equation 48    Figure 47 Schematic representation of the equivalent “rock-bar” and bolt system, and the displacements and fictitious forces along a single bar.  For the purpose of method validation, results where compared to results obtained through FEM models assuming varying properties. 106  Using the properties listed in Table 10 for Model #3, the axial force acting along the bolts using the proposed method are compared to results from RS2, as shown in Figure 48. It can be seen results obtained via the EBB method are smoother. The apparent reason for this is that in EBB the bolts are modeled as ideal springs, whereas in FEM the surrounding tunnel is explicitly modeled using plane-strain elements. Thus, inaccuracies due to mesh effects are more pronounced in FEM.      Figure 48 Comparison of axial forces along bolts computed in EBB and RS2. Table 14 shows the inputs and results for a bolted tunnel, with no liner.  Table 15 shows the constant inputs for a combined liner and bolt support system, and  Table 16 shows the varied inputs and results. It can be seen that results using the methodology proposed in the chapter are in good agreement with the results from the FEM models.  107  Table 14 Result comparison for bolted tunnels Max. Displacement [cm] Max. axial force [MN] Input parameters EBB RS2 EBB RS2 Bolt diameter [mm] Bolt length [m] E bolt [GPa] K ratio Tunnel radius [m] # 9.3 9.3 2.13 2.1 25 2 200 3 2 1 18.8 18.7 2.13 2.1 50 4 200 3 4 2 37.9 39.2 2.13 2.1 50 8 200 3 8 3 4.62 4.6 1 1 50 4 200 1 4 4 5.9 5.8 0.77 0.8 50 4 200 0.5 4 5 8.8 9.1 3.67 3.58 25 2 400 3 2 6 9 9.1 2.13 2.1 25 4 200 3 2 7 8.2 8.5 5.73 5.74 50 2 200 3 2 8  Table 15 Liner and bolt parameters used for combined liner-bolt support comparison Pretension force [MN] Bolt diameter [mm] E bolt [GPa] Liner thickness [cm] E liner [GPa] 0.3 25 200 40 30  Table 16 Result comparison of combined liner-bolt support system Max. Displacement [cm] Max. liner axial force [MN] Input parameters EBB RS2 EBB RS2 Lb K ratio Tunnel radius [m] 0.66 0.69 27.1 27 2 3 2 1.6 2 37.8 36.8 4 3 4 2.5 3 43.6 42.1 6 3 6 0.63 0.57 13.1 11.8 6 1 6   108   6.5. Derivation for Plastic Ground Throughout the last decades, great efforts have been invested in attempts to obtain closed-form solutions for the problem of a circular tunnel in plastic ground under non-uniform loading. Detournay (1986) derived an approximate semi-analytical solution for the elasto-plastic interface of a Mohr-Coulomb material. Detournay’s solution relies on a number of simplifying assumptions; mainly, that the tunnel is fully encircled by an elasto-plastic interface throughout the progression of plastic failure. The latter assumption renders the problem statically determinate. This topic of research is still of interest; in a recent paper, Lu et al. (2018) published a solution for the widely used Hoek-Brown failure criterion.  However, knowledge of the tunnel plastic displacements alone still provides an insufficient basis for calculating the actual stresses that develop in the support system. Detournay and St. John (1988) developed design charts based on Detournay’s solution (1986) in which they assumed simplified uniform pressure of the support in order to assess the ground-support interaction. In this section, the procedure outlined in Section ‎6.2 for elastic ground is used as a basis for treating a circular tunnel in plastic ground, so that the actual non-uniform ground-support interaction problem can be solved. The plastic displacements need to be known and can be derived from any of the available solutions, or found from a single finite-element analysis of an unsupported tunnel. In numerical analyses involving plasticity, where the stress-strain relationship is non-linear, an iterative solution process must be applied. It is clear that ground weakens as it deforms, as the zone of plastic failure increases. Hence, it is proposed to reduce the stiffness of the equivalent beam in each stage of analysis. In order to obtain the proper degree of reduction, the following steps are proposed: 1. The problem is divided into n stages, so that the non-linear problem is properly represented by n linear stages.  2. From the elastic analysis procedure discussed in Section ‎6.2, the fictitious forces are divided by the number of stages n. 109  3. The global stiffness matrix of the equivalent beam representing the ground is taken from the elastic procedure as well, where 𝑡𝑒𝑞 is initially equal to the same value given for the elastic assumption (see Equation 26). 4. The displacements for each stage are computed, initially equal to the elastic displacements; using an iterative process, the stiffness matrix representing the ground is gradually decreased, until the displacement results are approximately equal to the known plastic displacements.  It is important to note that the decrease in ground stiffness mentioned in Step 5 is non-trivial, as there are an infinite number of ways to decrease the stiffness of the equivalent beam. Preliminary analysis showed that decreasing 𝑡𝑒𝑞 by a single reduction factor failed to yield proper displacement profiles. Through a trial and error process, it was found that good results can be obtained by applying two different factors to 𝑡𝑒𝑞: one factor for the initial reduction of  𝑡𝑒𝑞; a second factor for a linear reduction, so that 𝑡𝑒𝑞 gradually decreases for each beam element, when the global stiffness matrix is constructed from wall to crown.  Once steps 1-5 are completed, the ground support interaction can be computed. This requires an additional computation process, as illustrated in Figure 49. It is emphasized that for non-hydrostatic loading conditions, the total ground-support interaction cannot be fully represented by convergence-confinement curves. The curves showed in Figure 49 represent a reference point corresponding to a specific location of the tunnel (e.g. the tunnel crown).  In this example, the problem was divided into 4 stages. In order to avoid unnecessary iterations, the load increment 𝑃𝑖 used for each stage is increased so that the support response at the reference point will intersect the points of displacement marked at each stage.  The computation is repeated until 𝛴𝑃𝑖 is equal to 𝑃𝑡𝑜𝑡𝑎𝑙, where 𝑃𝑡𝑜𝑡𝑎𝑙 corresponds to the initial load at the reference point; for the example shown in Figure 49, 𝑃𝑡𝑜𝑡𝑎𝑙 = 𝑃1 + 𝑃2 + 𝑃3.   110   Figure 49 Computational process for plastic ground-support interaction analysis  As an example, a circular tunnel in weak plastic rock obeying the Mohr-Coulomb failure criterion is computed. Table 17 lists the tunnel and rock parameters. In practice, tunnel support is regularly installed after some displacement has occurred, dependent primarily on the distance of the point of support installation to the unsupported tunnel face (Vlachopoulos and Diederichs, 2009). For this example, a relaxation factor of 25% is assumed. Accordingly, the problem was divided into 4 loading stages, with the support applied at stage 2.   Table 18 lists the unsupported tunnel crown and wall displacements for each stage, used as inputs for solving the problem. Calculations were repeated for 4 different liner thicknesses, of 200, 400, 600, and 800 mm. The Young's modulus and Poisson ratio of the liner were taken as 30 GPa and 0.2, respectively.   111   Table 17 Tunnel parameters   16 Tunnel diameter [m] 5 Field stress [MPa] 0.5 K ratio 4 Young's Modulus [GPa] 0.25 Poisson ratio  37 Friction angle [degrees] 0.65 Cohesion [MPa] 25 Relaxation [%]  Table 18 Tunnel crown and wall displacements per stage Wall displacement [mm] Crown displacement [mm] Stage 0.78 3.9 1 2.46 8.2 2 5.6 13.3 2 10 19.0 4  Table 19 Results for liner axial forces, bending moments, and displacements, obtained from numerical models in RS2 and the proposed EBB method Max. Displacement [mm] Max. Bending moments [kNm] Max. Axial force [MN] Liner thickness [mm] EBB RS2 EBB RS2 EBB RS2 16 15 3.3 3.5 7.7 7.5 200 14 13 27 28 11.4 10.8 400 13 12 89 90 13.8 13.3 600 12 11 206 220 16.1 15.2 800  Table 19 shows a comparison of results obtained from numerical models carried out in RS2 and the proposed EBB method. These results include maximum axial forces, bending moments, and displacements that develop in the liner. The maximum values in this case all occur at the tunnel crown. The results for EBB are in good agreement with results from the FEM models. The results obtained through EBB are consistently higher in this case, for the 112  reason that the tolerance criterion for the calibration process in which the equivalent beam stiffness was reduced, was set to the conservative side. Once the maximum axial force and bending moment that develop in the liner are known, the maximum inner stresses can then be calculated using the relationship from beam theory: 𝝈𝒎𝒂𝒙 =𝑵𝟏∗𝒕+𝑴𝑰𝒕 𝟐⁄   (51) Equation 49  where N is the axial force, M is the moment, I is the moment of inertia of the liner section, and t is the thickness of the liner. The maximum stress can then be used to assess the suitability of the support, by comparison to the allowable concrete stresses.   Based on the results listed in Table 19 and Equation 51, the liner maximum stresses are plotted and shown in Figure 50. Results obtained using the convergence-confinement method were computed and added to the plot. It can be seen that the stresses calculated using the convergence-confinement method, which assumes hydrostatic pressure, are significantly smaller. This is somewhat counter-intuitive, as the overall volume of loading is larger for the convergence-confinement method in this example (K ratio equal to 1 in the convergence-confinement method, vs. K ratio equal to 0.5 in the actual problem). For the case where the K ratio is greater than one, the difference in results will increase. This demonstrates the importance of considering the non-hydrostatic loading effects on tunnel support.   113   Figure 50 Comparison of liner maximum stresses computed using finite-element models in RS2, the proposed EBB method, and the Convergence-Confinement method        (denoted C-C in the legend)  In terms of solving and post-processing time, results for the proposed method are generated almost instantaneous; FEM results are calculated in less than 10 minutes. One may argue that the difference is not large enough to conclude there is a significant advantage in using the proposed analytical method. However, a large number of iterations would be required in order to assess the cost of different support systems under varying ground conditions.  For such extensive analyses the difference in the time required for both solving and analyzing results would increase considerably. A new approach which requires the application of extensive probabilistic analysis is presented in ‎7.   114   6.6. Horseshoe Shaped Tunnel  Much of the work related to tunnel-support interaction assumes a circular tunnel cross section. This simplifying assumption is not limited to hypothetical cases, as many tunnels are indeed circular, or near circular. Examples of such tunnels include circular TBM tunnels, and near circular mechanically excavated tunnels. However, many tunnels have non-circular cross sections. Tunnels constructed using conventional excavation methods by and large consist of horseshoe shaped cross sections. The term “horseshoe” here, does not refer to a certain definitive geometry. Horseshoe shaped tunnel geometries may vary and include different width to height ratios and non-constant arc curvatures (see Figure 51).  It is noted that although at the current state of practice TBMs produce circular cross sections, there have been some recent technological advances in developing non-circular TBMs. These machines are comprised of a combination of circular and/or quasi-circular faces (Li, 2017).     In general, the horseshoe shape has a number of advantages over the circular shape; mainly, the cross section area is better utilized, and the flat surface provides a more convenient working space. For determining the ideal tunnel cross section for a given project, different factors must be considered, such as functionality, permanent equipment installation (e.g. air supply ducts), constraints posed by the excavation method, etc. In addition to these constraints, the tunnel engineer may seek to modify the cross section due to stability considerations. For example, in some cases it is important to increase the tunnel curvature to avoid high stress concentrations at the edges of the tunnel.  115   Figure 51 Example of different two-lane road tunnel cross sections (modified from Thewes et al., 2013).   As outlined in Section ‎6.2, in the derivation process for the proposed EB method, the stiffness of the equivalent beam which represents the ground was found by assuming a circular cross section for the tunnel. The idealized circular shape allows for a simpler mathematical formulation. In order to address the problem of non-circular tunnels, it is possible that the stiffness expression for circular tunnels may not yield acceptable results.  The objective of this chapter is to assess the validity of the proposed methodology for circular tunnels for the analysis of horseshoe shaped tunnels, without further modification of the expression for the equivalent stiffness. Hence, the stiffness found for circular tunnels is used here as an input, with a single adjustment; this expressions is related to the tunnel diameter (see Equation 25); hence, the average between the total width and height of the horseshoe tunnel was used instead to compute the equivalent stiffness.  Additionally, the two following modifications are made with regards to the computational process: 1. The coordinates assigned for the analysis follow the horseshoe shape. 2. The unsupported tunnel displacements are taken from a preliminary FE model of the unsupported tunnel.  116  With respect to the latter modification, unlike the case of circular tunnels where closed-form solutions are used for computing the unsupported tunnel displacements, there are no such available solutions for horseshoe shaped tunnels. Therefore, it is required that an initial FE model is executed. This requirement may raise the question: if it is necessary to run a FE model, what is the practical benefit of applying the proposed method?  There are two possible answers to this question: 1. If in the future researchers will develop solutions for the displacements of horseshoe tunnels, these solutions could be readily integrated into the proposed method for the purpose of tunnel support analysis.  2. For cases where various ground conditions and/or support systems need to examined, the proposed method is greatly advantageous. In this case, a single FE model of the unsupported tunnel serves numerous instant computations of varying conditions with the EB method. This point will be demonstrated in the example given in Chapter ‎7.  In order to examine the suitability of the proposed methodology for horseshoe tunnels, a practical example is presented based on a large horseshoe shaped cavern currently being constructed in Menuhot Mountain, Jerusalem. Results from analysis using the proposed method are then compared with 2D numerical models for method validation. The authors have published a paper on the Menuhot project study case (Elmo et al., 2019). This paper focuses on a specific aspect of the project related to scale effects in rock mechanics. In this project, an array of small diameter holes is being drilled along some portions of the caverns' sidewalls (see Figure 52 and Figure 53). During construction, issues of instability were encountered through the process of the drilling of the small holes. In contrast, there were no problems associated with the excavation of the large caverns, and the good tunneling conditions allowed minimal support to be applied. The authors argued that the apparent explanation for this phenomenon stems from the influence of the different scale of the cavern and holes. A hybrid finite-discrete (FDEM) approach was used for modeling and demonstrating these effects. The analysis and conclusions from this paper will not be reviewed in this thesis. Instead, the relevant project information and data gathered from the site will be briefly reviewed here in order to provide sufficient background for the given example.  117  The “Menuhot” tunneling project is currently in its final stages of construction and includes a network of large horseshoe shaped caverns. These caverns are 14 m in width, and 16.5 m in height (see Figure 52). The caverns are situated in dolomitic limestone, and the topography provides for an average overburden of approximately 50 metres.  The excavation has been primarily executed by roadheaders with occasional drilling and blasting.  Table 20 lists the rock input parameters used for analysis. Figure 54 shows the typical geological profile.  Figure 52 “Menuhot” cavern cross section. 118   Figure 53 Drilling of small diameter holes in the cavern sidewalls  Table 20 Rock input parameters  Value Unit Uniaxial Strength 60 MPa Rock mass Young’s Modulus 30 GPa Hoek-Brown parameter mi 8 - Poisson’s ratio 0.23 - Unit weight 23 kN/m3  119   Figure 54 Typical geological profile cross-section.  As mentioned, the objective of the following study is to examine the suitability of the effect of the proposed EB method for horseshoe shaped tunnel-support interaction analysis. In order to avoid over complication, elastic behavior is assumed for the current study. It is noted that for the current project, this assumption is fairly reasonable, as the rock mass is considered to be massive, and the uniaxial strength is significantly greater than the in-situ stresses. A comparative study was carried out, where 4 different methods were used to compute the supported tunnel displacements:  2D numerical analysis via the program RS2.  The proposed EB method for the horseshoe shaped tunnel.  The proposed EB method for a circular shaped tunnel.  The convergence-confinement method. For the RS2 model, two materials were used: 1. Table 20 Dolomitic limestone, with the elastic properties listed in Table 20. This material was used from the level of the tunnel floor and upwards.  2. Stiff rock, with an unrealistic and exaggerated Young's modulus of 300 GPa. This material was used from the level of the tunnel floor and downwards (see Figure 55).  120  The reason for the unrealistic stiffening material from the floor level downwards relates to a constraint posed by the proposed EB method. As explained in Section ‎6.2, proper boundary constraints are necessary for obtaining meaningful results. From a mathematical perspective, boundary constraints reduce the number of unknown displacements, and render the problem solvable. For the circular tunnel, symmetry allowed for the application of fixed rollers at both the vertical and horizontal axes. In the case of the horseshoe tunnel, symmetry cannot be applied upon the horizontal axis; therefore, an alternative constraint is necessary. By using an exaggerated stiffness for the stiff rock material, near zero displacements will occur at the tunnel surface. Hence, a boundary constraint of zero vertical and horizontal displacements can be assigned to the corresponding node. It is argued that this is a reasonable assumption; some practitioners artificially stiffen the material of the tunnel surface as part of their standard practice, as in reality significant heaving at the tunnel floor seldom occurs.   121   Figure 55 RS2 cavern model  At the tunnel crown, the boundary constraints are identical to the circular tunnel, i.e. no horizontal displacement and no rotation can occur. Figure 56 shows these boundary constraints. In addition, Figure 56 shows the portion of the tunnel contour computed with EB, and the radii of curvature along this portion. Both the numerical models and the EB analysis for the horseshoe tunnel shape follow the exact shape of the actual tunnel cross section. For the analyses using the EB method for circular tunnels and the convergence-confinement method, a diameter of 15.25 m was assigned, which is the average of the tunnel width and height (14 m and 16.5 m, respectively). 122    Figure 56 Portion of tunnel contour computed with the EB method, and queried in RS2 for comparison. The boundary constraints for EB are shown as well.  For the comparative study, all parameters are kept constant, and only the liner thickness varies from 150 mm, to 250 mm and 350 mm. Note that these thicknesses coincide with the actual support classes used in the project, which have been designed to address good, fair and poor tunneling conditions, respectively. The RS2 model results are assumed to be most accurate, and the results of the other methods are assessed according to their deviation from these results.  6.6.1. Results and Discussion  The results of most interest are the maximum vertical displacement that occurs at the tunnel crown, and the maximum horizontal displacement that occurs at the tunnel wall. Table 21 and Table 22 shows the results of crown and wall displacements obtained via the different methods.    123  Results of crown displacements show that analyses using the EB method for both circular tunnels and horseshoe tunnels are in good agreement with the numerical models, with results of the circular tunnel analysis being slightly more accurate. Conversely, the displacements calculated by the convergence-confinement method are significantly smaller. Results of wall displacements show that the analysis using the EB method for the horseshoe shaped tunnel are in very good agreement with the numerical models. The results from the EB method for circular tunnels are significantly smaller, and the results obtained from the convergence-confinement method are considerably larger.  Figure 57 shows the results of total displacements compared along the full length of the query line (as defined in Figure 56), obtained by RS2 models and the EB method for the horseshoe shaped tunnel. It can be seen that results are in agreement, with the results obtained by the EB method being slightly smaller. These results indicate that by applying the proposed EB method upon horseshoe shaped tunnels, acceptable results can be attained with no need for further calibration of the equivalent beam stiffness. It is noted that as elastic material behavior is assumed for the rock mass, the selected material properties are roughly proportionate with the elastic modulus and in-situ stresses. Hence, conclusions based on the current example are not limited to the specific parameters of the selected case.  However, the conclusions are indeed limited to geometries that resemble the current case study. In this context, the horseshoe shape that was investigated consisted of a relatively small difference between its height and width. While this is the more common case, in some cases the tunnel height and width may be very different. Under such circumstances, it is possible that using the equivalent stiffness derived for circular tunnels would yield unacceptable results. 124  Table 21 Comparison of results using different methods, for maximum displacements at the tunnel crown  Max. crown displacement [mm] Liner thickness [mm] C-C EB- circular EB- horseshoe RS2 387 487 486 486 150 382 482 476 480 250 375 476 469 475 300  Table 22 Comparison of results using different methods, for maximum displacements at the tunnel wall    Max. wall displacement [mm] Liner thickness [mm] C-C EB circular EB horseshoe RS2 387 93 205 204 150 382 91 204 203 250 375 86 203 202 300  a)  125  b)  c)   Figure 57 Comparison of displacement results along the tunnel boundary (as defined in Figure 56) using RS2 and EBB, for different liner thickness of: a) 150 mm, b) 250 mm, and c) 350 mm     126   7. CHAPTER 7  METHODOLOGY FOR TUNNEL SUPPORT PROBABILISTIC ANALYSIS 7.1. Introduction From a practical perspective, the computational resources for FE analysis of a single tunnel are negligible, as models are set up and solved in a matter of minutes. However, for probabilistic analyses that require thousands of iterations, resources may become a considerable issue. In these cases, simplified solutions are highly advantageous.  This chapter presents a new methodology for the purpose of tunnel support cost estimation. As this topic is brought in the context of demonstrating the advantage of using the EB method subject of this thesis, it is also constitutes a stand-alone idea, i.e. this methodology can be used in conjunction with any available analytical or numerical method used for tunnel support analysis.    Perhaps the greatest challenge in underground engineering is dealing with the inherent uncertainty that stems from the heterogeneous nature of geological formations. In tunneling and underground caverns, the strength characteristics of the rock and/or soil dictate the support required for stabilization.  Naturally, the degree of uncertainty is greatest at the prefeasibility stage. Decision makers must determine whether or not to carry out a project, or to choose between different construction alternatives, based on very limited data. Considering the high cost and risks of underground projects, it is essential that reliable analysis tools are made available to them.  Currently, studies have shown that cost over-runs in tunnel projects are a global phenomenon. For example, Efron and Read (2012) examined 158 tunnel projects from 35 countries, and concluded that in every case the final cost was higher than the initial estimation. While there are a number of different risks that cause cost over-runs in tunneling projects, the costs of support have a highly significant impact on the overall project budget (Paraskevopoulou and Benardos, 2012).  127  Statistical tools and probabilistic analyses are widely applied in geotechnical practice for the purpose of coping with uncertainty. Probabilistic design tools in geotechnical engineering have been discussed and developed by many researchers, amongst them: (Hatzor and Goodman, 1993; Carter, 1992; Hoek, 2006).  Perhaps the most rigorous method used for managing uncertainty is the Monte Carlo method. The Monte-Carlo method is an iterative analysis process which uses input parameter values created by random number generators based on statistical distributions. This is in contrast to the more frequent and simple deterministic analysis, where single values are used as inputs.  The great advantage of the Monte-Carlo method is its ability to work with any type of distribution and variables, whether independent or not. The drawback of the Monte-Carlo method is that it requires a great number of iterations, compared to other probabilistic methods, such as the the Generalised Point Estimate Method, developed by Rosenbleuth (1981), and the First Order Reliability Method (Mirzaeian et al., 2015). In the context of geotechnical design, this analysis allows obtaining results of the Probability of Failure (PoF). Figure 58 shows an example of a PoF histogram plot, obtained via a Monte-Carlo analysis. The reliability of the system, or percentage of success, can be found by summation of the bars where the Factor of Safety (FoS) is greater than one (coloured blue in the figure). The PoF is the sum of the red bars on the failed side (FoS ≤ 1). The Factor of Safety is defined as the structural capacity over the demand. As an example, in a 1 km length tunnel, a PoF of 5% indicates that failure of a total length of 50 m would be expected. A clear and concise introduction to the theory of probability in the context of geological engineering is given by Hoek (www.rocscience.com/learning/hoek-s-corner).   128   Figure 58 Example of probability of failure histogram  It is important to be aware of human bias and false intuitions when dealing with such concepts as PoF and FoS. By the simple mentioning of words such as “safety” vs. “failure”, different messages may be conveyed. The term FoS creates an impression of safety redundancy, whereas the term PoF conveys the opposite; the PoF term implies that hazardous or even fatal outcomes are acceptable. Thus, organizations may be reluctant to use this definition. Nevertheless, Hoek argues that the PoF gives a more rational assessment of risks than the classical FoS. In this chapter, a methodology is proposed for utilizing the Monte Carlo analysis to better estimate the magnitude of support required for a tunnel. A distinction is made between tunnels that will be supported by a constant support, and tunnels where the support system varies, and is selected according to the ground conditions. 129  In order to demonstrate the implementation of the proposed methodology, a practical example is presented. The example is based on an analysis conducted for the comparison and cost estimation of different construction methods: TBM vs. SEM tunneling. In this example the tunnel is to be excavated through hard brittle rock susceptible to spalling. By the way of presenting the methodology and example, a brief review of the related topics is given, including the influence of tunneling methods on analysis, and brittle rock failure.  7.2. Tunneling Method Influence on Probabilistic Analysis Tunnels can be constructed using either conventional excavation methods or with a fully mechanized operation using a Tunnel Boring Machine (TBM). In conventional tunneling it is most common that the excavation and support are sequenced, and referred to by some as SEM (Sequential excavation method).    The choice of the optimal method, or combination of methods, depends on many factors, including tunnel length, cross section size, schedule, availability of skilled workers, the preparatory measures required in the site, and more. Different guidelines and rules of thumb have been published in order to help select the optimal construction method (Gütter et al., 2011; Hemphill, 2012). In many cases, a trivial decision cannot be made and a rigorous analysis is required in order to choose the favorable option. This decision is made either by the owner at the pre-feasibility stage, or by the contractor, when the bidding conditions enable both construction methods.  It is noted that in some cases, a constant tunneling support pattern is adopted for conventionally excavated tunnels. For these tunnels, the methodology proposed for TBM tunnels applies.  One important factor that is sometimes overlooked in the decision process is the magnitude of support required for both methods. In a TBM tunnel the support remains constant throughout the tunneling project, whereas in SEM a different support class is fitted for each step of advance. Hence, it is argued that different probabilistic analysis methods should be adopted for estimating the magnitude of support for each method.  130  The magnitude of support for a tunneling project will hereinafter be referred to as the Volume of Support (VoS). The VoS includes the overall volume of concrete (or shotcrete) for the tunnel lining, the overall volume of rock bolts, and of any other support components, such as steel sets and wire mesh. Obviously, knowledge of the VoS allows making a good estimate of the costs that will be invested in supporting the rock. In addition to the direct costs of support which include material purchase and the labor required for installation, indirect costs that are influenced by the VoS can also be assessed. These may include the cost effects of factors such as advance rate and logistics. The proposed methodology for probabilistic analysis allows better cost estimation of the VoS based on the tunneling method. The constant support in TBM tunnels renders the VoS simpler to calculate, and is: 𝑽𝒐𝑺𝑻𝑩𝑴 = 𝑻𝑳 ∗ 𝑺𝒖𝒑𝒑𝒐𝒓𝒕  𝒄𝒓𝒐𝒔𝒔 𝒔𝒆𝒄𝒕𝒊𝒐𝒏𝒂𝒍 𝒂𝒓𝒆𝒂  (52) 50 where TL is the tunnel length, and n is the number of iterations used in the Monte-Carlo analysis. The VoS for SEM tunnels is: 𝑽𝒐𝑺𝑺𝑬𝑴 = ∑ 𝑺𝒖𝒑𝒑𝒐𝒓𝒕 𝒄𝒓𝒐𝒔𝒔 𝒔𝒆𝒄𝒕𝒊𝒐𝒏𝒂𝒍 𝒂𝒓𝒆𝒂𝒏𝟏𝑻𝑳𝒏  (53) Equation 51  In each iteration of the SEM analysis, the support class is recalculated according to the anticipated ground properties that the tunnel will cross. The anticipated ground conditions are a function of the probabilistic analysis, i.e. they are computed via the Monte Carlo method. For finding the proper support system in each method a failure criterion must be defined. It is argued that the failure criterion most suitable for the support system in a TBM tunnel is a PoF with a near zero value. This is due to the fact that a single support pattern must have the sufficient capacity to tolerate even the most severe conditions. This requires careful consideration of the ranges of the statistical distributions selected for the ground parameters. If a near zero PoF is found to yield a very heavy and inefficient support system, than the cost of a limited extent of failure could be weighed and considered. For the SEM tunnels, it is suggested that a moderate FoS be considered. As tunneling progresses, the inner rock surfaces are revealed, mapped and classified. The support is then 131  selected according to this classification. In this process, the degree of uncertainty is reduced, but still does not completely diminish. There is always some degree of subjectivity and uncertainty involved in the rock classification process, and therefore the geological engineers should add a moderate margin of safety to their judgment.  The computational process for finding the VoS for each of the tunneling methods is as the following:   For TBM driven tunnels, the failure criteria and acceptable PoF are defined, as well as the ground parameters and their statistical distribution; a Monte-Carlo analysis is then executed in order to calculate the PoF for an assumed support class; the process is repeated until the assumed support class satisfies the acceptable PoF; finally, the VoS is found by multiplying the support class by the tunnel length, as given in Equation 52.    For SEM tunnels, the failure criteria and acceptable FoS are first defined, as well as the ground parameters and their statistical distribution; a Monte-Carlo simulation is executed, and in each iteration the support class that satisfies the FoS is found; the process is repeated, and the distribution of support is found; The VoS is than calculated according to the distribution of support, as given in Equation 53. Figure 59 shows flowchart diagrams for these two processes.     132  a)  133  b)  Figure 59 Flowchart diagram for: a) Tunnels with constant support (TBM tunnels)  b) Tunnels with varying support (SEM tunnels)  Figure 60 shows an example of a convergence-confinement plot for a single iteration within the SEM computation scheme. In this process, the tunnel convergence curve is first constructed based on the random ground parameters. Then, the support (or confinement) curve of support class #1 is constructed, consisting of an elastic component and a flat plastic component.  The linear portion of the support curve class #1 in this example does not intersect the convergence curve, and this indicates that the support will fail (FoS<1). Consequently, an 134  additional curve according to support class #2 is generated. This time the elastic component intersects the convergence curve; the FoS is calculated, and is equal to the ratio of the final support pressure to the support yield pressure. If this FoS satisfies the failure criteria, the iteration is completed, and subsequently an additional tunnel convergence curve is generated. It is emphasized that Figure 60 section is shown only for illustration, as the computation here will be undertaken using the proposed EB method, and not by the convergence-confinement method. As stated in Chapter ‎1, for non-hydrostatic loading, convergence-confinement plots cannot be generated. Nonetheless, the EB computation scheme follows the same general concept.  Figure 60 Example of convergence-confinement plot for SEM computation scheme  Once the processes shown in Figure 59 are completed, the VoS for each method can be obtained. It is noted that after completing this computation process and obtaining results, support classes can be reassumed in order to investigate other options that may be more cost effective. For example, if the distribution of support results show that results are skewed to the minimal support class, this may indicate that it would be more cost effective to create 135  lighter support classes. The total cost should then be compared to previous results in order to examine if the modified support classes are indeed beneficial.  Currently, this optimization process can be done manually, after viewing results, applying judgment, and modifying the inputs for the selected support classes. Additional work can be carried out to further improve the procedure so that the optimization process is an inherent part of the computational scheme. The probabilistic methodology proposed in this chapter requires numerous iterations. Nested loops are needed for finding the VoS of both tunneling construction methods, SEM and TBM (see Figure 59). Commercial FEM programs could be modified to carry out such analysis; however, this would require an extremely long duration of time for computing and interpreting the results.  The proposed EB method for tunnels is hereby used for demonstrating the implementation of the proposed methodology by a practical example. As shown in Section ‎6.2, the EB method requires as little as 50-100 elements for accurate computation of the tunnel-support interaction. Therefore, this method is significantly more efficient compared to FEM based numerical models, which require thousands of elements for a single analysis.   In the example presented in Section ‎7.4, a comparison of the solution time of the Monte Carlo simulation for EB and FEM is given. It is noted that the bottleneck for systematic utilization of the Monte Carlo approach is not necessarily the simulation time, but rather computer graphics and flow simulators in FEM codes that at the present are much too slow. This problem is avoided when simplified methods are implemented.   7.3. Brittle Failure Criterion When tunneling through massive hard rock under high stresses, brittle failure may occur. This is the assumption for the example given in the following Section ‎7.4. Brittle failure is manifested by sudden cracking and ejection of rock fragments. This phenomenon is often referred to as spalling, and in its more severe form as strainbursting, or rockbursting.  136  It is noted that the issue of rockbursting is not an essential part of both proposed methodologies presented here, i.e. the EB method for ground-support analysis and probabilistic analysis for support cost estimation. The issue of brittle failure is incorporated into the example in order to render the problem and results more interesting and less trivial.  Note that in the paper published as part of the work on this thesis, a simpler example was discussed. There, a tunnel in weak rock subjected to hydrostatic loading was assumed, and the standard convergence-confinement method was used for analysis (Mitelman and Elmo, 2019).  The following provides a brief discussion of brittle failure criterion. The spalling phenomenon has been widely studied, from a number of different aspects:  Mechanistic study- the spalling mechanism has been studied at the laboratory scale, and explained by fracture mechanics theory. Eberhardt et al. (1998) introduced the technique of acoustic emission measurements in order to find the crack initiation and propagation thresholds for rockburst occurrence.  Empirical prediction- based on documentation of several study cases, different methods are available for the estimation of spalling and strain burst potential, and spalling depth of failure, or notch size. One method for assessment of spalling susceptibility is based on the ratio of compressive to tensile strength, where the potential for spalling increases for strong rocks with a high compressive to tensile strength ratio (Diederichs, 2007).   Numerical modeling- given that spalling occurs in massive brittle rock, elastic modeling is suitable for the prediction of spalling failure (Kaiser 1996). A value of 0.6-0.8 of the UCS is recommended for numerical modeling. Attempts have been made to use different failure models in order to simulate the extent of failure (Diederichs, 2003).  For near-circular excavations, Martin et al. (1999) proposed that the Depth of Failure (DoF) near the wall can be obtained by a semi-empirical equation: 𝑫𝒐𝑭 = 𝒓 ∗ (𝟏. 𝟐𝟓𝝈𝒎𝒂𝒙𝝈𝒄− 𝟎. 𝟓𝟏 ∓ 𝟎. 𝟏)   (54) Equation 52  137  where r is the tunnel radius,  𝜎𝑚𝑎𝑥 is the maximum stress in the rock, and σc is the unconfined compressive strength of the intact rock. The DoF geometry is shown in Figure 61.   Figure 61 Depth of Failure (DoF) for a tunnel with a radius r   The uncertainty of both the occurrence and severity of spalling can be incorporated into any probabilistic analysis. For example, Lee et al. (2013) used fuzzy probability theory to systematically assess the probability of spalling in the early design stage. The FoS is generally defined as the capacity over the demand. For a tunnel support system this can be defined as the support material strength divided by the maximum stress in the support. In order to determine the FoS for the brittle failure criterion, two different design approaches can be used: 1. The stress induced in the rock is limited so that spalling failure is prevented, i.e.  𝜎1UCS< 0.4. This criterion is then used for determining the FoS of the rock.  2. The potential for spalling is permitted, and in turn, the support must have sufficient capacity to resist this effect. For this approach, The FoS is established based on the FoS of the support alone, without defining a limit for the rock. 138  A final FoS can be selected after considering both approaches, based on the lower result between the two. It is well known that some relaxation is inevitable prior to the installation of support, and the stresses in the stress concentration zones increase. Subsequently, additional movement will occur during the ground-support interaction process, and the rock stresses are again increased. Considering these two effects, the first design approach may not be cost effective, and in some cases even unfeasible.  Essentially, spalling and strainbursting are a dynamic process. Static analyses are significantly more computationally efficient than dynamic analyses, and therefore it is desirable to develop a pseudo-static design approach for the analysis of brittle failure, rather than conducting highly complicated dynamic analyses which simulate the actual spalling process. In order to achieve this, it is necessary to develop an expression for the Equivalent Static Force (ESF). Equivalent static forces are widely used in engineering. One such notable and widely used example is equivalent forces for seismic design of buildings.  The ESF as defined here, is the static load that would inflict a similar response upon the support system, in terms of displacements and stresses, as would the ejection of the failure notch at a given ejection velocity. The ESF is defined as: 𝑬𝑺𝑭 = 𝒏 ∗ 𝒎𝒏𝒐𝒕𝒄𝒉 ∗ 𝒈  (55) Equation 53  where n is an impact factor relating the dynamic load to the static load, 𝑚𝑛𝑜𝑡𝑐ℎ is the mass of the failure notch that is formed within the brittle failure process, and g is the acceleration of gravity. The mass of the failure notch is dependent on the DoF (see Figure 61).  For the general case of an impact problem, the impact factor n is equal to: 𝒏 = √𝝁𝒗𝟐𝒈∗𝜹𝒔𝒕𝒂𝒕𝒊𝒄    (56)  Equation 54  where v is the velocity of the impacting body, δstatic  is the static displacement due to a static force operating at the point of impact,   µ is a coefficient of efficiency which is influenced by different factors, such as inelastic energy loss, the ratio of the striking mass to the impacted 139  mass, and more. In general, for the impact of two hard elastic bodies, the value of µ is close to one; the theoretical lower bound impact factor is 2, and can be as high as 10, and even over 100 (Pytel and Kiusalaas, 2003).  Currently, there is a lack of data available to assist in the selection of the impact factor n for brittle failure spalling. The peak velocities measured in rock bursting are in the range of 3 m/sec, and can exceed 10 m/sec (Kaiser and Cai, 2012). However, these values of velocity would not be realistic when using Equation 54. A previous paper by the authors discussed the subject of spalling in rock induced by blasting (Mitelman and Elmo, 2015). There, it was argued that a much lower impact factor must be assumed for practical applications, due to a number of reasons:   Plastic deformations and fracturing processes are energy absorbing phenomena.  Additional energy losses occur within the spalling process due to inner collisions and rotation of fragments.  The spalled fragments of rock do not strike at the same time, as the deeper spalled layers are ejected with lower velocity.  More work based on field observations and data analysis is required for establishing guidelines for the selection of the proper impact factor and equivalent static force and their integration into the support design process.  In order to demonstrate the effect of the impact factor on results, a numerical investigation was carried out. The EB method was used for this investigation. Generally, one of the fundamental concepts in all methods for tunnel design is that the surrounding rock carries some portion of the load. The assumption for this analysis is that the equivalent static force is imposed upon the liner structure alone, as great portions of the rock during the spalling process unravel and lose their load bearing capacity.  A preliminary sensitivity analysis was conducted in order to examine the effect of the impact factor on the liner stresses. For this purpose, a circular tunnel with the mean values listed in Table 23 in Section ‎7.4, supported by a 15 cm thick concrete lining was computed. Based on the results of this analysis, Figure 62 shows a plot of the increase in maximum liner stress vs. the impact factor. As concrete strength is usually in the range of 25-50 MPa, these results 140  show that selection of the impact factor has a significant influence on the support FoS. This highlights the argument that additional work is essential for establishing proper guidelines for impact factor selection. Lacking this, conservative assumptions must be adopted.   Figure 62 Increase in liner stress vs. impact factor   The increase in stresses due to the spalling effect cannot be computed in a single and straightforward stress and strain analysis. The ESF depends on the dimensions of the failure notch; the dimensions of the notch depend on the maximum stresses that develop in the rock-support system; finally, the ESF affects the stress and strain analysis by inducing additional stresses upon the tunnel. Thus, the following computational process is proposed: 1. An initial stress and strain numerical analysis is carried out.  2. Based on the resultant maximum stresses calculated in Step 1, the tunnel radius, and the UCS, the Depth of Failure (DoF) is calculated according to Martin's empirical formula. 141  3. According to the DoF, the weight of the failure notch is computed. The ESF is then found by multiplying the weight of the notch by the impact factor, according to Equation 55.  4. The support stresses due to the application of the ESF are calculated. 5. The maximum support stress is determined by summation of the maximum stresses computed in Steps 1 and 4.   6. The FoS is determined, and is equal to the stress capacity of the support divided by the maximum stress calculated in the Step 5. This computational process is repeated for each iteration within the probabilistic analysis described in the following section.  7.4. Practical Tunnel Analysis Example Two twin tunnels are to be excavated through brittle dolomite under an overburden of 250-700 meters. The dolomite formation is assumed to be generally massive, with few widely spaced discontinuities. It is assessed that under deeply covered areas, spalling failure may initiate. The project can be executed as an SEM operation by drilling and blasting, by TBM tunneling, or by some combination of the two methods. With respect to the type of TBM, it is assumed that a shielded TBM is necessary, due to the regions where the rock is more fractured (Hemphill, 2012). The required dimensions of the tunnel are 7.60X7.80 m (for SEM), or a diameter of 8.0 m (for TBM), as shown in Figure 63. For the SEM alternative, initial numerical models carried out showed that the differences between results of the elliptical tunnel shape and a circular tunnel with an average diameter of 7.70 m are negligible. The distance between the tunnels is assumed to be great enough so that the interaction between the two tunnels can also be neglected. 142   Figure 63 Tunnel required inner dimensions (in mm units) for SEM (solid line). The dashed line represents the circular TBM cross section of an 8 m diameter.  Based on a site investigation and laboratory tests, mean values and standard deviations given in Table 23 are assumed for probabilistic analysis. All parameters are assumed to be normally distributed. In order to avoid numerical instability, the normal distribution was truncated to two standard deviations from the mean values.      143  Table 23 Parameters and standard deviations for analysis Parameter Mean value Standard deviation GSI 70 10 UCS [MPa] 80 10 Overburden [m] 450 125 Rock weight [MN/m³] 0.027 - Poisson’s ratio 0.2 0.02 K ratio 2 0.1 Relaxation factor SEM: 0.45 TBM: 0.5 0.02  The overburden too is assigned a distribution, due to two reasons:  1. When tunneling through mountainous environment, the overburden pressure is influenced by variations in topography and geology. 2. At the preliminary stage of a tunneling project, it is common that the exact route of the tunnel has yet to be determined. The relaxation factors given in Table 23 refer to the percentage of relaxation that occurs in the rock prior to support installation. There are different methods for assessing the relaxation factor, but since none are exact it would be sensible to apply a statistical distribution to this factor as well. As this factor depends on the timing of installation, it is clear that the different tunneling method dictate a different relaxation factor. For the TBM, this factor relies on machine operation, whereas in SEM this factor can be controlled. Hence, a lower factor is selected for the SEM analysis, because under brittle failure conditions support should be installed as immediate as possible. In the Monte Carlo analysis, the deformation modulus is computed based on the UCS and GSI randomly generated in each iteration. Different empirical equations have been developed for the rock mass deformation modulus. For this example, the empirical equation proposed by Hoek et al. (2002) is used, where the rock mass deformation modulus 𝐸𝑟𝑚 is: 𝑬𝒓𝒎(𝑮𝑷𝒂) = (𝟏 −𝑫𝟐)√𝝈𝒄𝒊𝟏𝟎𝟎𝟏𝟎(𝑮𝑺𝑰−𝟏𝟎𝟒𝟎)   (57) 144  Equation 55  where D is a disturbance factor related to the quality of blasting, and assumed zero for this work.  The spalling failure criterion outlined in Section ‎7.3  is used. An impact factor of 10 is used in all analyses for this example. It is possible to incorporate the impact factor into the probabilistic process; however, as stated, the knowledge required for the proper selection of the impact factor is generally lacking; Therefore, it is argued that doing so does not contribute to attaining better results. The acceptable PoF for the TBM analysis is assumed to equal 3%, and the acceptable FoS for the SEM analysis is assumed to be 1.15. In order to carry out the analysis, support classes must first be assumed. For the SEM, the support classes are given in Table 24. For the sake of simplicity, a one passage liner is assumed for SEM, i.e. no additional layer is considered as load carrying beyond the primary shotcrete liner. All bolts are assumed to be passive (no pre-tensioning), fully grouted, with a 25 mm diameter and a length of 3 m. The increase in liner compressive strength can be achieved by either using different strength concrete, and/or by the addition of reinforcement, such as fibers and rebar.  The final column in Table 24 shows the cost estimation for each support class, measured in dollars/m of tunnel length. These estimates are based upon unit prices for each support component (shotcrete liner, bolts and steel sets). These unit prices include the total estimated costs, i.e. cost of both material and labor.   145  Table 24 Support classes for SEM modeling Support class Liner thickness [mm] Liner compressive strength [MPa] Bolt pattern [mxm] Comments Cost estimation [$/m] 1 50 30 2.5x2.5 - 1,010 2 100 35 2x2 - 1,920 3 150 40 1.5x1.5 - 3,285 4 200 45 1x1 - 5,850 5 250 60 1.5x1.5 Addition of steel sets: 310 mm I-beam every 0.8 m 14,340  For the final design of segmental lining for TBM operation several factors must be considered (Herrenknecht and Bäppler, 2003). For this preliminary analysis, it was assumed that the minimal properties for the liner are a thickness of 30 cm and strength of 40 MPa, which are required in order to withstand handling loads, such as lifting and stacking. As part of the computational process, the liner strength was increased by increments of 2.5 cm in thickness and 5 MPa in strength, until the PoF criterion is satisfied.  In shielded TBM driven tunnels the lining is made of precast concrete segments. In numerical modeling, the connections between the lining segments can be implemented by adding hinges. The TBM liner was assumed to consist of 6 segments, and hinge elements were added to calculations accordingly. Hinges dictate a zero moment transfer from one element to the other. The hinge beam elements should be added into the global stiffness matrix of the liner according to the locations that correspond to their actual position in the tunnel. The element matrix 𝐾𝐻𝑖𝑛𝑔𝑒 for a hinged beam element is:  𝑲𝑯𝒊𝒏𝒈𝒆 = 𝑲𝒄 𝒍𝟑⁄ [𝟑       𝟑𝒍      − 𝟑       𝟎 𝟑       𝟑𝒍𝟐    − 𝟑𝒍      𝟎−𝟑 − 𝟑𝒍      − 𝟑     𝟎𝟎        𝟎             𝟎      𝟎]  (58) Equation 56  where l is the element length and 𝐾𝑐 is the flexibility coefficient (see Equation 39).  146  A constant bolt pattern of 2x2 m was used for all calculations. For the bolts, the price units for the SEM analysis were multiplied by a factor of 0.75, as the installation process in TBM is fully mechanized and therefore cheaper. For the TBM liner, the price units were multiplied by a factor of 1.25, due to manufacturing and handling expenses associated with prefabricated elements. The resultant distribution into support classes for the SEM alternative is given in Figure 64. Based on this distribution, the cost of support is calculated (see Table 25). It is noted that these results can be further used in order to estimate the project schedule. By assuming the advance rate for the different support classes, the results of support distribution according to their length can be divided by the advance rate to find the estimated project schedule.   Figure 64 Support class distribution per 1 km of tunnel for the SEM alternative  For the TBM alternative, a 40 cm thick liner with concrete strength of 50 MPa together with a 2X2 m bolt pattern was found to yield a PoF equal to 2.6%, and was therefore selected as the tunnel support system. Figure 65 shows the FoS distribution for this support system.    147   Figure 65 Probability of Failure of support system for the TBM alternative  Table 25 shows the cost estimation per kilometer of support for both options. The cost of support for SEM is lower, as the TBM support is fitted for the lower bound parameters, and in a SEM project better ground conditions allow the selection of lighter support classes. In addition to the difference in direct costs, the volume of liner required for both options shows that the TBM alternative requires 7,424 additional cubic meters of rock to be excavated for each kilometer of tunnel. The excavated rock generates additional costs associated with mucking and haulage, equipment wear, and other logistics. It is emphasized again, that the proposed methodology refers solely to costs associated with excavation and support, and that there are additional crucial factors that must be weighed in order to compare the two tunneling methods.    148  Table 25 Cost estimation result comparison Tunneling method Cost of support   [M$/km] Volume of liner [m³] TBM 8.89 10,556 SEM 3.48 3,132 Difference 5.41 7,424   Tunnel support costs may vary depending on a number of factors and assumptions. Costs will increase as the degree of heterogeneity of the rock increases. In our example, if the values of standard deviation given in Table 23 are reduced by a factor of 0.5, the support costs per kilometre would decrease by 25% and 18% for the TBM and SEM alternatives, respectively. Support costs also depend on the selected failure criterion (FoS and PoF), the partitioning of support into classes for the SEM, and in some cases the minimal thickness for the TBM liner segments. For the Monte Carlo analysis 10,000 trials were used. Due to nested iterations (as shown in Figure 59), the number of times a stress and strain computation was carried out was 40,000 for the TBM analysis, and over 50,000 for the SEM analysis. For comparison, a Monte Carlo analysis of 50,000 trials was executed via the program RS2 (RocScience, 2015). The FEM model for a circular tunnel in elastic ground, consisted of 8,828 elements. Both computations were carried out on a personal computer with a 3.2 GHz processor and 16 GB RAM.  A comparison of the running time and file size is brought in Table 26, where it can be seen that there is a considerable difference for both parameters. The implication of file size is slow computer response when viewing results in FEM. This considerable difference demonstrates the advantage of adopting a simplified computation method for Monte Carlo analyses.  149  Table 26 Comparison of the required computational resources for Monte Carlo analysis of EB vs. FEM (RS2)  EB RS2 Running time 3 minutes 60.5 hours File size 9 KB 353.5 GB  Additional and more complex factors and calculations may be incorporated into such analyses, such as additional failure criteria, more sophisticated pricing techniques, and more. Also, additional conclusions may be drawn from the results, such as schedule estimation and distribution to support classes in SEM.    150  8. CHAPTER 8  CONCLUSIONS AND RECOMMENDATIONS This thesis has presented a novel approach termed ‘Equivalent Boundary’ (EB) for the analysis of ground-support interaction problems. The basic idea of the proposed method is to simplify these problems by representing the ground with analogous structural entities. The approach allows for increased efficiency by focusing on the ground-support boundary, rather than simulating a great portion of the surrounding ground, as is required in finite-element models. 8.1. Summary of key contributions The EB method developed by the author was applied to study different ground-support problems, including room and pillar mines, vertical circular shafts, and tunnels:  The mine-pillar system is idealized to a spring system. According to this idealization, a formula for finding the maximum force acting in the pillar is derived. This formula can be used for evaluating the effect of support on pillar factor of safety. This idealization is verified using FEM analysis and also compared to the Middleton mine study case. It is shown that the tributary method traditionally used for initial pillar assessment is highly over-conservative, since it does not account for the room and pillar interaction. It is shown that in shallow mines bolting of the roof can add a significant contribution to the pillar FoS. A general and conceptual methodology is proposed and outlined for pillar post-peak analysis based on a distinction between two stages: the dynamic phase and residual state. In the dynamic phase the pillar support must withstand spalling. In the residual state, the pillar must continue to support the room after its strength has been reduced.   In circular vertical shafts, the liner behaves as a ring subjected to axisymmetric pressure. An equivalent boundary ring is found through a calibration process, so that the problem can be solved as a double ring system. For both the pillar and shaft problem, elastic and homogenous ground conditions were assumed.   For tunnel support analysis, an equivalent beam is used to represent the ground surrounding the tunnel boundary. This allows converting the plane-strain problem to a simpler double beam problem. This approach represents an enhancement of the 151  convergence-confinement method to include the common and more complex condition of non-hydrostatic loading conditions. Similar to the process derived for pillar and shaft analysis, the solution for circular tunnels relies on the deformations of the unsupported tunnel as an input. For elastic analysis, the Kirsch equations were used to find the unsupported tunnel deformations. For plastic analysis, the solution derived by Detournay (1986) can be used. The procedure for solving the plastic problem is based on a staged reduction of the initial stiffness of the equivalent boundary beam through an iterative trial and error process. Different support systems can be incorporated into the proposed method, including steel sets, grouted bolts, mechanically anchored bolts, and pre-tensioned bolts. Some types of support can be modeled using a straightforward approach (e.g. steel sets), while other types require a more complex derivation process (e.g. grouted bolts).  By comparing different methods, this thesis has demonstrated that the EB method yields results that are considerably more accurate compared to the standard convergence-confinement method. Moreover, results for elastic analysis are found to be in better agreement with precise analytical solutions compared to FE models.   The methodology for circular tunnels was modified to examine its applicability for the analysis of horseshoe shaped tunnels. As opposed to circular tunnels, there are no available solutions for the displacements of horseshoe tunnels. Therefore, an initial FE model is required so that the displacements of the unsupported tunnel are used as an input. A horseshoe cross section of a large cavern from an actual project was analyzed. Based on the results, it was found that the expression for the equivalent beam derived for circular tunnels yields acceptable results for the horseshoe tunnel.  In general, closed-form solutions for unsupported excavations do not account for the effect of ground-support interaction, and therefore cannot be used for the practical purpose of support design. The proposed EB method provides a framework for any available and future solutions so that the knowledge of the unsupported displacement is utilized for ground support design.     Finally, mmethodologies are proposed for the analysis of constant support projects (TBM tunnels) vs. adaptive support projects (SEM), for the purpose of tunnel support cost 152  estimation. Utilizing the Monte Carlo analysis method for the assessment of support for SEM tunnels, the anticipated support distribution is found, and then used to compute the estimated costs. It is shown that the difference in support costs between methods may be significant, and is influenced by the degree of heterogeneity of the ground.  While this proposed methodology can be implemented using any available tunnel support analysis method, the proposed EB method developed for this thesis is applied to calculations, and its advantage in terms of computational resources is clearly demonstrated.     8.2. Recommendation for Further Studies This thesis focuses mainly on ideal geometries. Nevertheless, there are still many practical instances where results can differ significantly from idealized geometries.  In order to extend the proposed method for solving such problems, a number of difficulties must be addressed. First, for these problems, there may be no available solutions for the unsupported displacements. Additionally, estimation of the stiffness of the equivalent elements of irregular excavation shapes is presumably more difficult to obtain than for idealized geometries. The feasibility of such modification is questioned.  Nevertheless, it is possible that rather than attempting to derive analytical or semi-analytical solutions, these problems could be solved by developing efficient algorithms that calibrate the equivalent boundary to match FEM results. The development of such algorithms may be a computing challenge in and of itself.    Figure 66 shows a flowchart diagram proposed for such cases. This process could be developed as a subroutine for FEM programs dedicated for the purpose of ground-support problems in geotechnical analysis.  In addition, it is possible that this approach may be suitable for additional engineering problems that include the interaction between two distinct bodies, where the larger body is converted to an equivalent and simplified body. For example, for the problem of the impact of a small body with a half-infinite surface, the surface could be converted to an equivalent beam. It is noted that such a problem includes a further degree of complexity, as it is a dynamic problem. 153    Figure 66 Flowchart for generalized application of equivalent boundaries   In addition to the general approach presented in Figure 66, further studies to investigate the application of the methodology presented in this thesis may include specific ground-support problems such as: 154   Multiple room and pillar mines- in this thesis, the problem of the interaction of a single pillar in a single room was addressed. A more accurate solution would account for the scenario of multiple pillars. This would require developing an equivalent slab to represent the rock surface interacting with the pillars. Such an analogue would reduce the 3D problem into a simpler "slab-on-springs" problem, which requires considerably less computational resources.  Retaining walls- the problem of retaining walls is typically addressed using simple analytical solutions, such as the Coulomb or Rankine solutions, Winkler spring solutions, or FEM analysis. Similar to the tunnel problem, deriving an equivalent beam to represent the ground at the ground-wall interface may lead to more accurate solutions than solutions obtained by deriving equivalent springs, and without significant increase in computational resources. This is due to the continuous nature of the beam, as opposed to independent springs.    Foundation slabs- in practice, Winkler springs are most commonly used for foundation slab analysis. Deriving an equivalent slab may contribute to obtaining more accurate results.  Three-dimensional tunnel problems- while for most tunnel design applications the 2D plane-strain assumption is valid, there are different specific cases where a 3D solution is necessary. For example, in very weak and plastic rock, it is essential to account for the three-dimensional effects of the interaction between the tunnel face and walls. At the current state of practice, 3D models are often time consuming in terms of both construction of the geometry, and computer solving. In order to simplify this problem using the methodology proposed in this thesis, an equivalent 3D surface can be constructed, so that the problem remains a three-dimensional problem, but has significantly less elements.  To summarize, for any applications that include analysis of the interaction between different bodies, it may be efficient to convert the boundary of the larger body into an equivalent interface.  Additional work can be carried out to address secondary issues. This work could include: 155   Further developing the initial framework discussed in the context of post-peak pillar failure analysis. For the dynamic phase, field data of pillar spalling velocities and spalling volume are required in order to design support with sufficient capacity to absorb the spalling impact; for the residual state, improved guidelines for the selection of post-peak strength based on fracturing field observations and residual displacements are required.  Enhancing the proposed methodology for tunnel support cost estimation to include a process of support optimization.  Developing a method for finding the equivalent static force for loads induced by brittle failure in tunnels in hard rock.     156  References American Association of State Highway and Transportation Officials (AASHTO). 2010. Technical Manual for Design and Construction of Road Tunnels. AASHTO, Sep. 2010 Edition. Anandarajah, A., 2010.Computational methods in elasticity and plasticity: Solids and porous media. First ed. New York: Springer.  Barczak, T. M., Esterhuizen, G. S. and Dolinar, D. R. 2005. Evaluation of the impact of standing support on ground behavior in longwall tailgates. Proceedings of the 24th International Conference on Ground Control in Mining, 23–32. Barton, N.R., Lien, R. and Lunde, J. 1974. Engineering classification of rock masses for the design of tunnel support. Rock Mech. 6(4), 189-239. Bieniawski, Z. T. 1989. Engineering rock mass classifications: a complete manual for engineers and geologists in mining, civil, and petroleum engineering . New York : Wiley .  Bischoff J. A. and Smart J. D. 1975. A method of computing a rock reinforcement system which is structurally equivalent to an internal support system. Proc. of 16th Syrup. on Rock Mechanics, Univ. of Minnesota, 179-184. Brady BHG., Brown ET. 2006. Rock mechanics for underground mining. Springer, 3rd Edition, Netherlands. Brown, E.T., Bray J.W., Ladanyi B., Hoek E. 1983. Ground response curves for rock tunnels. J. Geotechnical Eng. ASCE, 109: 15-39.  Burland, J. B. 1995. Assessment of risk of damage to buildings due to tunneling and excavation. Int. Earthquake Geotechnical Engineering, Tokyo, 3,1189–1201. Carranza-Torres C., Diederichs M., 2009. Mechanical analysis of circular liners with particular reference to composite supports. Tunnelling and Underground Space Technology, 24, 506–532 Carranza-Torres, C. 2004.Elasto-plastic solution of tunnel problems using the generalized form of the Hoek-Brown failure criterion. International Journal of Rock Mechanics and Mining Sciences, 41(1), 1-11 Carranza-Torres, C., Fairhurst, C. 1999. General formulation of the elasto-plastic response of openings in rock using the Hoek-Brown failure criterion. International Journal of Rock Mechanics and Mining Sciences. 36 (6), 777-809. 157  Carranza-Torres, C., Fairhurst, C., 2000.Application of the convergence–confinement method of tunnel design to rock masses that satisfy the Hoek–Brown failure criterion. Tunnelling and Underground Space Technology 15(2), 187–213 Carranza-Torres, C., Rysdahl B., and Kasim M. 2013. On the elastic analysis of a circular lined tunnel considering the delayed installation of the support. International Journal of Rock Mechanics and Mining Sciences, 61, 57-85. Carter, T.G. 1992. Prediction and uncertainties in geological engineering and rock mass characterization assessments. Proc. 4th. int. rock mechanics and rock engineering conf., Torino.  Chern, J.C., Shiao, F.Y., Yu, C.W. 1998. An empirical safety criterion for tunnel construction. Proceedings of the Regional Symposium on Sedimentary Rock Engineering, Taipei, Taiwan. 222–227. Crowder J. J., Bawden W.F. 2004. Review of post-peak parameters and behaviour of rock masses: current trends and research. Rocnews, Fall, 13. Cui, L., Zheng, J.-J., Zhang, R.-J., & Lai, H.-J. 2015. A numerical procedure for the fictitious support pressure in the application of the convergence–confinement method for circular tunnel design . International Journal of Rock Mechanics and Mining Sciences, 75, 336-349. Detournay E., St. John. 1988. Design charts for a deep circular tunnel under non-uniform loading. Rock Mechanics & Rock Engineering, 21, 119-137. Detournay, E. 1986. An approximate statical solution of the elastoplastic interface for the problem of Galin with a cohesive-frictional material. International Journal of Solids and Structures, 22(12), 1435–1454.  Diederichs, M. S. 2003. Manuel Rocha Medal Recipient Rock Fracture and Collapse Under Low Confinement Conditions. Rock Mechanics and Rock Engineering, 36(5), 339–381.  Diederichs, M. S. 2007. The 2003 Canadian Geotechnical Colloquium: Mechanistic interpretation and practical application of damage and spalling prediction criteria for deep tunnelling. Canadian Geotechnical Journal, 44(9), 1082–1116.  Dobry R, Gazetas G, Tassoulas JL. 1985. Vertical Response of Arbitrarily Shaped Embedded Foundations. Journal of Geotechnical Engineering. 111(6), 750-771. Duncan Fama, M. E. 1993. Numerical Modeling of Yield Zones in Weak Rock. Analysis and Design Methods. Oxford: Pergamon, 49-75 158  Eberhardt, E., Stead, D., Stimpson, B., & Read, R. S. 1998. Identifying crack initiation and propagation thresholds in brittle rock. Canadian Geotechnical Journal, 35(2), 222–233.  Efron, N., & Read, M. 2012. Analysing International Tunnel Costs. An Interactive Qualifying Project. Report, Worcester Polytechnic Institute, USA. Einstein H.H., Schwartz C.W. 1979. Simplified analysis for tunnel supports. Journal of the Geotechnical Engineering Division, American Society of Civil Engineers 105, 499–518 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. PhD thesis, University of Exeter, Cornwall Elmo D., Stead D. 2010. An integrated numerical modeling– discrete fracture network approach applied to the characterization of rock mass strength of naturally fractured pillars. Rock Mechanics and Rock Engineering, 43, 3–19  Elmo, D., Mitelman, A. Rozen, A., & Tzuker, D. 2019.Design and Construction of a Dense Array of Small Diameter Holes in a Large Cavern Sidewall . 53rd U.S. Rock Mechanics/Geomechanics Symposium. June, New York  City, New York. Esterhuizen, G.S., Mark C., Murphy M. 2010. The ground response curve, pillar loading and pillar failure in coal mines. Proceedings of the 29th International Conference on Ground Control in Mining, Morgantown, WV, 19-27. Esterhuizen G. S., Dolinar D. R., Ellenberger, J. L. 2011. Pillar strength in underground stone mines in the United States. International Journal of Rock Mechanics and Mining Sciences, 48, 1, 42-50. Exadaktylos G.E., Stavropoulou M.C., 2002. A closed-form elastic solution for stresses and displacements around tunnels. International Journal of Rock Mechanics & Mining Sciences, 39, 905–916 Fang, Q., Zhang, D., & Wong, L. N. Y.2012. Shallow tunnelling method (STM) for subway station construction in soft ground. Tunnelling and Underground Space Technology, 29, 10–30.  Felippa, Carlos A. 2001.A historical outline of matrix structural analysis: a play in three acts. Computers & Structures, 79 (14), 1313–1324.  Ferreira A.J.M. 2008. MATLAB codes for finite element analysis: solids and structures. Springer Science & Business Media, Vol. 157. 159  Filonenko-Borodich M. 1946, On a certain system of functions and its application in the theory of elasticity. Prikl. Mat. Mekh. 10, 193-208.  Gazetas G. 1991. Shear modulus of rockfill from lateral pile load test. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, 28(6). Gazetas, G; Dobry, R; Tassoulas, J. 1986. Vertical response of arbitrarily shaped embedded foundations. Journal of Geotechnical Engineering, 111(6), 750–771.  Gerolymos N, Gazetas G. 2006. Development of Winkler model for static and dynamic response of caisson foundations with soil and interface nonlinearities. Soil Dynamics and Earthquake Engineering. 26(5), 363-376. Grasso P, Mahtab A, Ferrero AM, Pelizza S. 1991. The role of cable bolting in ground reinforcement. Soil and Rock Improvement in Underground Works, 1,127-38. Gütter, W., Jäger, M., Rudigier, G., Weber, W. (2011). TBM versus NATM from the contractor’s point of view / TBM versus NÖT aus Sicht des Unternehmers. Geomechanics and Tunnelling, 4(4), 327–336.  Hatzor. Y., Goodman. R.E. 1993. Determination of the ‘design block’ for tunnel supports in highly jointed rock. Analysis and design methods, Oxford: Pergamon, 2, 263-292. Hedley D.G.F, Grant F. 1972. Stope-and-pillar design for the Elliot Lake uranium mines. Bull Canadian Institution of Mining and Metallurgy, 65, 37-44. Hemphill, G. B. 2012. Practical Tunnel Construction. Wiley.  Herrenknecht, M., & Bäppler, K. 2003. Segmental concrete lining design and installation. Soft Ground and Hard Rock Mechanical Tunneling Technology Seminar. Hoek, E., Marinos, P., Kazilis, N., Angistalis, G., Rahaniotis, N., & Marinos, V.. 2006. Greece's Egnatia Highway Tunnels. Tunnels & Tunneling International, 32-35. Hoek, E. 1989. A limit equilibrium analysis of surface crown pillar stability. Surface crown pillar evaluation for active and abandoned metal mines, 3-13.  Hoek, E. and Brown, E.T. 1997. Practical estimates or rock mass strength. International Journal of Rock Mechanics and Mining Sciences, 34(8), 1165-1186. Hoek, E., Marinos, P. 2000. Predicting Tunnel Squeezing. Tunnels and Tunnelling International. Tunnels and tunnelling international, 32(11), 45-51. Hoek, E., Carranza-Torres, C. and Corkum, B. (2002). Hoek-Brown failure criterion – 2002 Edition. 5th North American Rock Mechanics Symposium and 17th Tunneling Association of Canada Conference, 1(1), 267-271. 160  Hoek, E., Carranza-Torres, C., Diederichs, M.S., Corkum, B. 2008. Kersten Lecture: Integration of geotechnical and structural design in tunnelling. Proceedings University of Minnesota 56th Annual Geotechnical Engineering Conference, Minneapolis, 1-53 Hoek, E., Carter, T. G. & Diederichs, M. S. 2013. Quantification of the Geological Strength Index chart. 47th U.S. Rock Mechanics/Geomechanics Symposium. Hoek, E. Practical rock engineering. Chapter 7: A slope stability problem in Hong Kong. Hoek’s corner, Internet edition. www.rocscience.com/learning/hoek-s-corner Hoek, E. 1998. Practical rock engineering. Chapter 8: Factor of safety and probability of failure. Hoek’s corner, Internet edition. www.rocscience.com/learning/hoek-s-corner Hung, C. J., Monsees, J. E., Munfah, N., Wisniewski, J. 2009. Technical Manual for Design and Construction of Road Tunnels - Civil Elements, AASHTO, 1–702. Itasca Consulting Group, Inc. 2016. FLAC — Fast Lagrangian Analysis of Continua, Ver. 8.0. Minneapolis: Itasca. Jaeger, J. C., N. G. W. Cook. 1969. Fundamentals of rock mechanics. Wiley, Methuen, London, 395-399;  Jibson., Randall W. 2011. Methods for assessing the stability of slopes during earthquakes- a retrospective. Engineering Geology, 122(1).  Kaiser P. K, Cai., M., 2012. Design of rock support system under rockburst condition, Journal of Rock Mechanics and Geotechnical Engineering, 4 (3), 215–227 Kaiser, P. K., MacCreath D. R., Tannant D.D. 1996. Canadian rockburst support handbook : prepared for Sponsors of the Canadian Rockburst Research Program 1990-1995. Geomechanics Research Centre. Sudbury, Ontario.  Kaiser, P. K., Cai M. 2012. Design of rock support system under rockburst condition. Journal of Rock Mechanics and Geotechnical Engineering, 4(3), 215–227.  Kerr A.B. 1964. Elastic and viscoelastic foundation models. Journal of Applied Mechanics, 491–498. Kirsch C. 1898. Die Theorie der Elastizität und die Bedürfnisse der Festigkeitslehre. Zeitschrift des VereinesdeutscherIngenieure, 42, 797–807. Krauland N., Soder P.E. 1987. Determining pillar strength from pillar failure observations. E&MJ-Engineering and Mining Journal, 188(8), 34-40. 161  Lee, K. H., Bang, J. H., Lee, I. M., Shin, Y. J. 2013. Use of fuzzy probability theory to assess spalling occurrence in underground openings. International Journal of Rock Mechanics and Mining Sciences, 64, 60–67.  Leonhardt, F. 1973. Vorlesungen uber Massivbau. Berlin: Springer, Verlag. Lewin, K. 1945. The research center for group dynamics at Massachusetts Institute of Technology. Sociometry, 8(2), 126-136. Li, J. 2017. Key Technologies and Applications of the Design and Manufacturing of Non-Circular TBMs. Engineering, 3(6), 905-914. Lu, A., Wang, S., Zhang X., Zhang N. 2018. Solution of the elasto-plastic interface of circular tunnels in Hoek–Brown media subjected to non-hydrostatic stress. International Journal of Rock Mechanics and Mining Sciences, 106, 124–132.  Lunder PJ, Pakalnis R. 1997. Determination of the strength of hard-rock mine pillars. Bull. Canadian Institution of Mining and Metallurgy Bulletin. 90, 51-55. Martin C.D., Maybee W.G. 2000. The strength of hard-rock pillars. International Journal of Rock Mechanics & Mining Sciences, 37(8), 1239-1246. Martin, C. D., Kaiser, P. K., McCreath, D. R. 1999. Hoek-Brown parameters for predicting the depth of brittle failure around tunnels. Canadian Geotechnical Journal, 36(1), 136-151. Mathworks 2012. MATLAB. User's guide: http://www.mathworks.com/products/matlab 2012 Mirzaeian Y., Shahriar K., Sharifzadeh, M. 2015. Tunnel Probabilistic Structural Analysis Using the FORM. Journal of Geological Research. Mitelman A., Elmo D. 2014. Modelling of blast-induced damage in tunnels using a hybrid finite-discrete numerical approach. Journal of Rock Mechanics and Geotechnical Engineering, 6(6), 565-573.  Mitelman A., Elmo D. 2016. Analysis of tunnel support design to withstand spalling induced by blasting. Tunnelling and Underground Space Technology, 51, 354-361. Mitelman, A. & Elmo, D. 2018. A Proposed Probabilistic Analysis Methodology for Tunnel Support Cost Estimation Depending on the Construction Method. 52nd U.S. Rock Mechanics/Geomechanics Symposium 8, Seattle, Washington. 162  Mitelman, A. Elmo, D. 2019. Analysis of tunnel-support interaction using an equivalent boundary beam . Tunnelling and Underground Space Technology incorporating Trenchless Technology Research  84, 218–226. Mitelman A., Elmo D., Stead, D. 2018. Development of a spring analogue approach for the study of pillars and shafts. International journal of Mining Science and Technology, 28(2), 267-274. Ngan V., Broere, W., Bosch, J. W. 2017. Structural analysis for shallow tunnels in soft soils. International Journal of Geomechanics, 17(8), 04017038.  Nie, W., Zhao, Z. Y., Ning, Y. J., & Guo, W. 2014. Numerical studies on rockbolts mechanism using 2D discontinuous deformation analysis. Tunnelling and underground space technology, 41, 223-233. Oravecz, K.I. 1997. Analogue modelling of stresses and displacements in bord and pillar workings of coal mines. International Journal of Rock Mechanics and Mining Sciences & Geomechanics, 14(1), 7-23. Oreste P.P., Peila D. 1996. Radial passive rockbolting in tunnelling design with a new convergence-confinement model. International journal of rock mechanics and mining sciences & geomechanics abstracts, Pergamon, 33(5), 443-454.  Oreste, P. P. 2003. Analysis of structural interaction in tunnels using the convergence–confinement approach. Tunnelling and Underground Space Technology, 18(4), 347-363.  Oreste, P. P. 2007. A numerical approach to the hyperstatic reaction method for the dimensioning of tunnel supports. Tunnelling and Underground Space Technology, 22(2), 185-205.  Ozturk H., Guler E. 2016. A methodology for lining design of circular mine shafts in different rock masses. International Journal of Mining Science and Technology, 26, 76-8.  Panet, M. 1993. Understanding deformations in tunnels. Comprehensive Rock Engineering, Pergamon, 1,663–690. Paraskevopoulou, C., Benardos, A. 2013. Assessing the construction cost of Greek transportation tunnel projects. Tunnelling and Underground Space Technology, 38, 497–505.  Pasternak P.L. 1954. New method of calculation for flexible substructures on two-parameter elastic foundation, Gosudarstvennoe Izdatelstoo, Literaturi po Stroitelstvu Arkhitekture, Moskau (in Russian), 1–56. 163  Peck, R.B., Hanson W.E., Thorburn T.H. 1974. Foundation engineering, 2nd edn. New York: Wiley Pells P.J., Bieniawski Z.T., Hencher S.R., Pells S.E. 2017. Rock quality designation (RQD): time to rest in peace. Canadian Geotechnical Journal, 54(6), 825-834  Pytel, A., Kiusalaas, J. 2003. Mechanics of Materials. Brooks/Cole Thomson Learning.  Ritter, W. 1879. Die Statik der Tunnelgewölbe. Berlin: Springer. Roberts D.P., Merwe J.N., Canbulat I., Sellers E.J., Coetzer S. 2002. Development of a method to estimate coal pillar loading. Safety in Mines Research Advisory Committee. Report No. 2001-0651. Rockfield (2007) Rockfield Software Ltd. Technium, Kings Road, Prince of Wales Dock, Swansea, SA1 8PH, UK. http://www.rockfield.co.uk Rocscience Inc. 2015, RocSupport 4.0 -Rock support interaction and deformation analysis for tunnels in weak rock. www.rocscience.com, Toronto, Ontario, Canada. Rocscience Inc. 2015, RS2 (Phase2) Version 9.0 - Finite Element Analysis for Excavations and Slopes.www.rocscience.com, Toronto, Ontario, Canada. Rosenblueth, E. 1981. Two-point estimates in probabilities. Applied Mathematical Modelling, 5(5), 329-335.  Schurch R., Anagnostou G., 2012. Applicability of the ground response curve to tunneling problems that violate rotational symmetry. Rock Mechanics & Rock Engineering, 45,1–10 Simler, K., Hanson R. 2018. The Elephant in the Brain: Hidden Motives in Everyday Life Kindle Edition. New York: Oxford University Press. Sjoberg J. 1992. Failure modes and pillar behaviour in the Zinkgruvan mine. Proceedings of 33rd US Rock Mechanics Symposium, Santa Fe, Balkema, Rotterdam, 491-500. Sober, E. 2015. Ockham's razors: A user’s manual. Cambridge University Press. Svoboda, T., & Mašín, D. 2010. Convergence-confinement method for simulating NATM tunnels evaluated by comparison with full 3D simulations. Proceedings of International Conference Underground Construction, Prague, 795–801. Terzaghi, K. 1946. Rock defects and loads in tunnel supports. Rock tunneling with steel supports. R.V. Proctor and T.L. White, Harvard University. Terzaghi, K. 1943. Theoretical soil mechanics, Wiley, New York, 202-206. 164  Thewes, M., Maidl, B., Maidl, U., Sturge, D. S. 2013. Handbook of Tunnel Engineering I: Structures and Methods . Berlin, Wilhelm Ernst & Sohn, Wiley.   Itasca Consulting Group. 2005. UDEC Manual Version 4.0.  Vlachopoulos N., Diederichs M.S. 2009a. Improved longitudinal displacement profiles for convergence-confinement analysis of deep tunnels. Rock Mechanics and Rock Engineering, 42(2), 131-146. Vrakas A., Anagnostou G. 2014. A finite strain closed-form solution for the elastoplastic ground response curve in tunnelling. International Journal for Numerical and Analytical Methods in Geomechanics, 38(11), 1131–1148. WALLAP Version 5. 2007. User Manual for Propped, Anchored and Cantilever Retaining Wall Analysis Program. Geosolve, UK. Baumgart, F. 2000. Stiffness-an unknown world of mechanical science. Injury-International Journal for the Care of the Injured, 31(2), 14-23 Wong, R.C.K., Kaiser, P.K. 1988. Design and performance evaluation of vertical shaft: rational shaft design method and verification of design method. Canadian Geotechnical Journal. 25, 320–37 Zipf, R. K. 2001. Pillar Design to Prevent Collapse of Room-and-Pillar Mines. Underground Mining Methods: Engineering Fundamentals and International Case Studies. Society for Mining Metallurgy and Exploration (SME), 493–511.    165  Appendix 1  A.1 MATLAB Formulation for Circular Tunnels in Elastic Ground In MATLAB coding, a basic m-file is a plain text file with a list of commands to be executed by the computer. The percent sign % represents comments, which are not executed. The following m-file computes the final displacements of a circular tunnel supported by a liner. The user may alter the ground and liner variables as he or she pleases.  Following the m-file is a function script. Functions in MATLAB are coded to accept inputs and return outputs. Functions are generally used for convenience, so that the basic m-file is not overloaded with commands, and that specific tasks are executed via a separate independent file.  Here, the m-file calls for the function "formStiffnessTunnel" to assemble the stiffness matrices of the equivalent boundary beam and tunnel liner.    % EB code for computation of circular lined tunnel in elastic medium clear variables  % Ground properties E=10e9; % [Pa] Ground Young's modulus r=2.5; % [m] Tunnel radius k=0.5; % k ratio poisson=0.25; % Ground Poisson ratio h=1.0*r; % Equivalent thickness   I=h^3/12; % Equivalent inertial moment (for bending). A=r/(1+poisson); %  Equivalent cross section area. EI=E*I; EA=E*A; G=E/(2*(1+poisson)); % Shear modulus p=10e6; % Vertical in-situ pressure   %Liner properties E_l=30e9; % [Pa] Liner Young's modulus poisson_l=0.2;% Liner Poisson ratio h_l=0.3; % Liner thickness I_l=h_l^3/12; % Liner inertial moment   A_l=h_l; EI_l=E_l*I_l; EA_l=E_l*A_l; 166  EI_c=EI+EI_l; EA_c=EA+EA_l;         numberElements=100; % number of elements.  % Variables definition for coordinates and unsupported displacements alpha=0.5*pi/(numberElements); % angle increment  teta_n=0:alpha:(0.5*pi); % angle variables ll=r*alpha; % element length  disp=zeros(1,numberElements+1); % unsupported radial displacements u_teta=zeros(1,numberElements+1); % unsupported tangential displacements disp_y=zeros(1,numberElements+1); % unsupported vertical displacements disp_x=zeros(1,numberElements+1); % unsupported horizontal displacements disp_xy=zeros(1,numberElements+1); % unsupported total displacements nodeCoordinates=zeros(numberElements+1,2); % node coordinates   % Loop for generating coordinates and displacements according to Kirsch  for i=1:numberElements+1      nodeCoordinates(i,1:2)=[r*cos(teta_n(i)) r*sin(teta_n(i));];   disp(i)=-(p*r/(4*G)*((1+k)-(1-k)*(4*(1-poisson)-1)*...   cos(2*teta_n(i)))); % radial displacements  u_teta(i)=-p*r/(4*G)*((1-k)*(2*(1-2*poisson)+1)*sin(2*teta_n(i)));% tangential displacements disp_y(i)=disp(i)*sin(teta_n(i))+u_teta(i)*cos(teta_n(i)); disp_x(i)=disp(i)*cos(teta_n(i))-u_teta(i)*sin(teta_n(i)); disp_xy(i)=sqrt(disp_y(i)^2+disp_x(i)^2);  end   slope=(numberElements/(0.5*pi*r))*gradient(disp); % rotational displacements slope(1)=0; slope(end)=0;   xx=nodeCoordinates; for i=1:numberElements     elementNodes(i,1)=i;     elementNodes(i,2)=i+1; end      numberNodes=size(nodeCoordinates,1);      xx=nodeCoordinates(:,1);      yy=nodeCoordinates(:,2);      GDof=3*numberNodes; % General degrees of freedom              % Formulation of vector of unsupported displacements         disp_r=zeros(1,GDof);       disp_c=zeros(1,GDof);         disp_r(1:3:end)=disp_x(1:end);       disp_r(2:3:end)=disp_y(1:end);       disp_r(3:3:end)=slope; % stiffness function, including transformation from local to global % coordinates       [stiffness_r,stiffness_c,lls]=... 167           formStiffnessTunnel(GDof,numberElements,elementNodes,numberNodes,xx,yy,EI,EA,EI_c,EA_c,EA_l,EI_l); % calls for function which computed stiffness matrices of the ground and liner       % boundary constraints      prescribedDof=[2;3;GDof-2;GDof]; activeDof=setdiff([1:GDof]',[prescribedDof]);   % Computation of (fictitious) force vector  force=stiffness_r*disp_r'; %    % Computation of supported tunnel displacements disp_c(activeDof)=stiffness_c(activeDof,activeDof)\force(activeDof);  disp_cx=disp_c(1:3:end); % Supported horizontal displacements disp_cy=disp_c(2:3:end); % Supported vertical displacements disp_cxy=sqrt(disp_cx.^2+disp_cy.^2); % Supported total displacements xxx=1:length(disp); % Variable for plots     % Computation of coordinates of deformed shape  deformCoordinates=zeros(numberElements+1,2);  delta_ax=zeros(numberElements,1); lll=ll*(r-h_l)/r ;% liner element length, considering liner thickness for i=1:numberElements+1      deformCoordinates(i,1:2)=[(r)*cos(teta_n(i))+disp_cx(i); (r)*sin(teta_n(i))+disp_cy(i);];      if i>1          delta_ax(i-1)=sqrt((deformCoordinates(i,1)-deformCoordinates(i-1,1))^2+...              (deformCoordinates(i,2)-deformCoordinates(i-1,2))^2)-ll;      end       end xxd= deformCoordinates(:,1); yyd= deformCoordinates(:,2);    function [ stiffness_r,stiffness_c,kl ] =    formStiffnessTunnel(GDof,numberElements,elementNodes,numberNodes,xx,yy,EI,EA,EI_c,EA_c,EA_l,EI_l); % Calls for stiffness matrix assembly function     %Formulation of stiffness matrices of unsupported and supported tunnel, % including conversion of local to global coordinate system.  stiffness_r=zeros(GDof); % Unsupported tunnel stiffness matrix stiffness_c=zeros(GDof); % Supported tunnel stiffness matrix   for e=1:numberElements       indice=elementNodes(e,:);     elementDof=[3*e-2 3*e-1 3*e 3*e+1 3*e+2 3*e+3]; nn=length(indice); 168  xa=xx(indice(2))-xx(indice(1)); ya=yy(indice(2))-yy(indice(1)); length_element=sqrt(xa^2+ya^2); cosa=xa/length_element; sena=ya/length_element; ll=length_element;  L=[cosa sena 0 0 0 0;    -sena cosa 0 0 0 0;     0 0 1 0 0 0;     0 0 0 cosa sena 0;     0 0 0 -sena cosa 0;     0 0 0 0 0 1;];    kr=[EA/ll 0 0 -EA/ll 0 0;     0 12*EI/ll^3 6*EI/ll^2 0 -12*EI/ll^3 6*EI/ll^2;    0 6*EI/ll^2 4*EI/ll 0 -6*EI/ll^2 2*EI/ll;     -EA/ll 0 0  EA/ll  0 0;     0  -12*EI/ll^3 -6*EI/ll^2 0 12*EI/ll^3 -6*EI/ll^2;     0 6*EI/ll^2 2*EI/ll 0  -6*EI/ll^2 4*EI/ll ];    kc=[EA_c/ll 0 0 -EA_c/ll 0 0;     0 12*EI_c/ll^3 6*EI_c/ll^2 0 -12*EI_c/ll^3 6*EI_c/ll^2;    0 6*EI_c/ll^2 4*EI_c/ll 0 -6*EI_c/ll^2 2*EI_c/ll;     -EA_c/ll 0 0  EA_c/ll  0 0;     0  -12*EI_c/ll^3 -6*EI_c/ll^2 0 12*EI_c/ll^3 -6*EI_c/ll^2;     0 6*EI_c/ll^2 2*EI_c/ll 0  -6*EI_c/ll^2 4*EI_c/ll ];    stiffness_r(elementDof,elementDof)=stiffness_r(elementDof,elementDof)+L'*kr*L; stiffness_c(elementDof,elementDof)=stiffness_c(elementDof,elementDof)+L'*kc*L; end   

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}]}"
                            data-media="{[{embed.selectedMedia}]}"
                            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:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0390677/manifest

Comment

Related Items