Abstract
Tunnelling processes lead to stress changes surrounding an underground opening resulting in the disturbance and potential damage of the surrounding ground. Especially, when it comes to hard rocks at great depths, the rockmass is more likely to respond in a brittle manner during the excavation. Continuum numerical modelling and discontinuum techniques have been employed in order to capture the complex nature of fracture initiation and propagation at low-confinement conditions surrounding an underground opening. In the present study, the hybrid finite-discrete element method (FDEM) is used and compared to techniques using the finite element method (FEM), in order to investigate the efficiency of these methods in simulating brittle fracturing. The numerical models are calibrated based on data and observations from the Underground Research Laboratory (URL) Test Tunnel, located in Manitoba, Canada. Following the comparison of these models, additional analyses are performed by integrating discrete fracture network (DFN) geometries in order to examine the effect of the explicit simulation of joints in brittle rockmasses. The results show that in both cases, the FDEM method is more capable of capturing the highly damaged zone (HDZ) and the excavation damaged zone (EDZ) compared to results of continuum numerical techniques in such excavations.
1. Introduction
Significant changes in the stress regime and material properties of a rockmass are the result of the construction of underground openings [1]. Induced stresses because of an excavation and subsequent stress redistributions result in damage, fracturing, desaturation, and so on, hence leading to the creation of a zone around the underground opening that the rockmass is disturbed with properties (strength, deformability, permeability, etc.) that differ from the original material properties [2]. Understanding the geomaterial response under such conditions is critical to the design and construction of an underground project, especially at great depths. Therefore, assessment of the damaged zone around an excavation is of great significance in evaluating the effects on the support requirements and the underground opening stability.
The damage zone surrounding an underground opening is comprised of different subzones depending on the intensity of the induced damage, but these zones are usually referred to collectively as the excavation damage zone (EDZ). The intensity of the damage usually depends on the distance of an examined point of the rockmass from the excavation, with the increase in distance resulting in the decrease in the influence of the excavation to the surrounding rock [3]. By disregarding rockmass damage associated with the construction method employed, as it can be eliminated by taking necessary precautions [4, 5], damage associated with stress changes, the excavation geometrical features, preexistent joints within the rockmass, and so on can be evaluated based on fracture coalescence [6] and divided into the highly damaged zone (HDZ), the excavation damage zone (EDZ), and the excavation influence zone (EIZ) [3].
Given the specific requirements of a project, numerical modelling can be utilized in order to provide a better insight of the rockmass response during the excavation and assist in the geotechnical and geological design. Especially, regarding the numerical simulation of rock fracturing, this has been the focus of interest for various researchers [7–14]. However, the conducted numerical modelling is only as good as the applied input parameters and assumptions that are integrated into the model, which is required to be able to simulate the expected rockmass behaviour under specific conditions. Within this paper, the short-term mechanical response of an underground excavation within a hard, highly interlocked rockmass is examined by employing two different numerical techniques, the finite element method (FEM) and the hybrid finite-discrete element method (FDEM). For the developed FEM models, the constitutive assumptions adopted include the use of the Hoek–Brown criterion [15] and the damage initiation-spalling limit (DISL) model as proposed by Diederichs [16], while for the developed FDEM model, the finite-discrete element method as proposed by Munjiza [17] is adopted. Furthermore, two different scenarios including a “fracture-free” and a fractured model are utilized by integrating a discrete fracture network (DFN) within the numerical models. For the properties of the intact rock, the well-established case of the Atomic Energy of Canada Limited’s (AECL) Underground Research Laboratory (URL) [18] is used in order to examine the mechanical response of the geomaterial during the excavation. Based on the obtained numerical results, it is shown that continuum techniques can capture brittle failure in hard rockmasses at some extent, by making appropriate modifications in the constitutive assumptions. However, certain aspects such as the clear distinction between the highly damaged zone (HDZ) (collapsed material) and the excavation damage zone (EDZ) (fractured material that maintains its structural integrity and does not collapse) are not captured by continuum approaches adequately. On the contrary, the FDEM method appears to be a better fit for simulating brittle fracturing, as it is capable of capturing the complex phenomena associated with brittle failure in hard rock excavations in low-confinement environments.
2. Short-Term Mechanical Response of Hard Rockmasses and Constitutive Assumptions for Continuum Codes
Instability of underground openings occurs as a result of gravity-driven fallouts which are controlled by the rockmass structure or are stress driven as the rockmass strength is exceeded. In both cases, the rockmass behaviour is controlled by two major factors including the in situ stress regime and the rockmass degree of fracturing [19]. Common assumptions within the engineering design include full persistence of joints. However, at great depths, nonpersistent jointing environments away from zones where tectonic processes take place (faulting, shear zones, folding, etc.) are more likely to be encountered. Therefore, the presence of rock bridges and the existent joints are the two contributing factors controlling the behaviour of massive and/or rockmasses with nonpersistent joints.
For cases in which stress driven rockmass failure is expected to occur, constitutive models that are based on the shear strength of the examined materials are commonly applied. However, their ability to capture the behaviour of massive or moderately fractured rockmasses around underground excavations has been proven to be limited. As documented in the literature, damage within hard rocks is the result of extensile fracturing parallel to the direction of the maximum principal stress , as a result of exceeding the tensile strength of the rock [16, 19–21]. This results in a lower overall rockmass strength observed in situ, which cannot be predicted by shear failure-based criteria which consistently overestimate the rockmass strength in the numerical models; hence, shear failure based-criteria are not appropriate to use in such cases and other techniques need to be employed.
The Hoek–Brown failure criterion [22] (1) is best suited for rockmasses which can be described as ductile (elastic-perfectly plastic) or rockmasses exhibiting strain weakening (postyield strength decreases) [3].where and are the major and minor principal stresses, respectively, is the unconfined compressive strength of the intact rock, and is the reduced value of the material constant for the intact rock according to the following:where is the Geological Strength Index [23] and is a factor which depends on the degree of the ground disturbance to which the rockmass has been subjected to blast damage and stress relaxation.
and are constants for the rockmass given by the following equations:
For underground excavations, such conditions arise for shallow openings, and therefore low confinement conditions (as long as structurally driven failures such as wedge failures, unravelling, etc. are not expected), or deep openings for which the rockmass strength is relatively lower to the in situ stress regime (e.g., squeezing ground conditions). However, as previously mentioned, for brittle rockmasses with nonpersistent joints and of high strength under high stresses, the Hoek–Brown criterion is not appropriate to use.
Brittle damage as a result of induced stresses during an excavation is commonly simulated by employing a cohesion weakening and friction mobilization approach [19, 24–27], with this approach being used for massive and moderately jointed rocks. Diederichs [16] developed the DISL method in order to capture the response of brittle rockmasses by using the generalized Hoek–Brown criterion [22] as a base by using (1). This modified approach utilizes the peak and residual strength envelopes of the Hoek–Brown criterion, based on the parameters listed in Table 1, in order to simulate brittle fracturing within conventional and commonly available numerical packages. The method involves the damage initiation as representation of an elevated cohesion and low friction, transitioning to the spalling limit which is represented by cohesion loss and friction mobilization, in order to capture the brittle behaviour of the rockmass at low confinement environments around an underground opening as the tunnel advances. In Figure 1, the Hoek–Brown and DISL strength envelopes are demonstrated.
In this study, the Hoek–Brown criterion and the DISL method are used within the FEM numerical software RS2 [28], in order to be compared to the developed FDEM model as described in the following sections.
3. The Finite-Discrete Element Method
The finite-discrete element method (FDEM) [17, 29–31] is a numerical technique that combines the finite element method with a smeared crack model in order to capture the behaviour of systems that involve discontinuum mechanics under complex deformation, rotation, interaction, and fracturing conditions [32]. The method is based on the discretization of the medium domain into 3-noded, triangular elements that form the mesh. The constant strain triangular elements that comprise the mesh are connected by 4-node cohesive elements which are employed in order to simulate the nonlinear material behaviour ahead of the crack tip due to interlocking and microcracking [33]. The assignment of these cohesive elements between the edges of the adjacent triangular elastic elements allows for the initiation and propagation of fractures once their tensile strength (opening—Mode I), their shear strength (sliding—Mode II), or both (mixed mode—Mode I-II) are exceeded, depending on the stress and deformation state of the medium. The main advantage of this type of simulation is that the potential fracturing paths do not need to be determined a priori. On the contrary, the trajectories of the fractures are controlled by the paths of the induced stresses imposed by the loading conditions. The fracture trajectories though depend on the mesh topology, a limitation that can be overcome as long as the selected element size is small enough so that the fracture patterns can be considered independent of the mesh configuration [34–36].
Prior to the occurrence of any fracturing, and as long as the strength of the cohesive elements is not exceeded, in order to maintain the elastic response of the material, an artificial stiffness for the cohesive elements is introduced in the simulated system by using normal, tangential, and fracture penalty coefficients [13, 17]. Regarding the strength of the cohesive elements, this is defined by their tensile and shear strengths as previously mentioned. The tensile strength of the cohesive elements is controlled by the peak tensile strength and the fracture energy in tension . In a similar fashion, the shear strength is controlled by the peak shear strength and the fracture energy in shear . The peak shear strength is expressed as the friction coefficient , the cohesion , and the normal stress , as show in (4) [37]:
Once the cohesive elements reach their peak strength in Mode I or Mode II, they yield and energy is dissipated through the assigned fracture energy values depending on the failure mode of the cohesive element. Once the fracture energy is depleted, the cohesive element breaks and is removed from the simulation. The forces that act on the nodes of the triangular elements govern the motion of the finite elements according to (5) [17]:where is the nodal mass matrix, is the nodal displacement vector, are the internal nodal forces as a result of the deformation of the constant strain elements, and are the external nodal forces that are generated by the imposed external loads, cohesive bonding forces because of the deformation of the unbroken cohesive elements, and contact forces due to the interaction of the broken cohesive elements. The motion equations of the system are solved by an explicit time integration scheme. Once a cohesive element breaks, the two previously connected triangular elements, referred as the contactor and target elements, respectively, interact with one another along the fracture based on the penalty function method [31]:where is the contact force, is the outward unit normal to the penetration boundary , and and are the potential functions for the contactor and target elements, respectively.
As described in the previous section, the primary form of damage of hard rock materials is that of extensile fracturing because of defects and flaws within the geomaterial [21, 38, 39]. The progressive failure of hard rock materials under the low confinement conditions that are encountered around an excavation boundary can be captured from an FDEM model with crack initiation and propagation explicitly according to the nonlinear-elastic fracture mechanics principles [40, 41].
For the purposes of this study, in order to capture the brittle behaviour of hard rock excavations at great depths, the FDEM model was developed in the 2D FDEM numerical code Irazu [42], as described in the following sections.
4. Geological Conditions
In order to study the brittle behaviour of a hard rockmass during the excavation of an underground opening, the material properties of the Lac du Bonnet (LdB) granite and information collected from the case study of the URL Test Tunnel, located in Pinawa, Manitoba, Canada were used [18]. The URL Test Tunnel was excavated at a depth of 420 m comprising of a circular cross section of a 3.5 m diameter within the LdB granite in a virtually fracture-free environment [43]. Overall, the LdB granite at that depth and within the vicinity of the tunnel can be assumed as a hard, massive, brittle rockmass that is homogeneous and isotropic. Based on data obtained from laboratory testing on intact specimens [44], the mechanical properties of the intact LdB Granite are listed in Table 2. Additionally, the stresses measured in situ at the URL Test Tunnel [43], as shown in Table 3, were used for the FEM and FDEM models in this study.
5. Numerical Model Setup
5.1. Geometry, Mesh Configuration, and Excavation Sequence
The geometry of both the FEM and FDEM models that were developed follow the geometrical characteristics of the URL Test Tunnel. In Figure 2, it can be seen that both models are comprised of a 60 m × 60 m master domain, with the FEM model having a smaller 15 m × 15 m domain in which the employed mesh is finer. The size of the outer domain was selected as such so that potential boundary effects would be avoided [45]. In a similar fashion, the FDEM model is also comprised of smaller subdomains with elements varying in size per domain in order to optimize the computational cost without compromising the accuracy of the obtained results. The subdomain close to the vicinity of the excavation consists of approximately 0.03 m elements in order to secure that the potential fracture patterns would be independent of the mesh topology, as previously mentioned. The specifics of the models are listed in Table 4. The equations of motion for the discretized system were integrated with a time step of 6.1 × 10−8 s for the intact model and 4.2 × 10−8 s for the fractured model in order to ensure numerical stability for the explicit solver of the code.
Regarding the simulation of the tunnel advancement, the induced three-dimensional (3D) effects on the rockmass are simulated by applying the excavation-induced stresses as a distributed load on the excavation boundary for the FEM models, and the face replacement method [45] for the FDEM models. For the FEM models, after the initialization of the geostatic stresses, the excavation material is removed and substituted with a uniform load applied on the tunnel boundary. This load is gradually reduced to zero in order to cause the destressing effect of the tunnel excavation. For the FDEM models, the advancement of the tunnel involves the replacement of the material within the excavation boundary with unstressed, elastic material during each step, hence leading to the “softening” of the tunnel core and the tunnel circumference to converge. The stability of the numerical analyses conducted both with the FEM and the FDEM models was secured by employing a relatively large number of steps. For the FEM models, the simulation process was comprised of twenty stages in order to simulate the advancement of the excavation face, and at each stage, the induced stresses were decreased by 5% until no load was applied on the excavation boundary. For the FDEM models, a large number of steps was employed in order to secure the equilibrium at the establishment of the geostatic stress stage (600,000 steps), the elimination of dynamic effects during the advancing excavation face (2,200,000 steps), and the subsequent complete removal of the material within the excavation boundary (700,000 steps).
5.2. Field Stresses and Boundary Conditions
The field stresses assigned to the tunnel models replicate the in situ stress conditions encountered during the excavation of the URL Test Tunnel [43]. For the purposes of this study, only the in-plane stresses were used for the FDEM models since the adopted cohesive elements do not account for the influence of the out-of-plane stress on fracture nucleation and growth. On the contrary, the FEM models were assigned an out-of-plane stress of the same magnitude as the recorded intermediate principal stress in the URL Test Tunnel, which was 45 MPa. The effect of gravitational forces and gravity-induced stress gradients were not taken into consideration in the numerical model, and a uniform constant stress field was assigned instead (deep tunnel assumption). Regarding the far-field boundary conditions, displacements are fixed in the horizontal and vertical directions for both the FDEM and FEM models. Additionally, for the case of the FDEM model, an absorbing boundary condition was employed in order to minimize dynamic oscillations.
5.3. Mechanical Input Parameters, Model Calibration, and Discrete Fracture Networks
5.3.1. Deformability Parameters
In order to determine the deformability parameters that are required for the numerical analysis, the properties of the LdB granite and the in situ rockmass conditions at the URL were assessed. Since the rockmass within the tunnel was virtually free of fractures, it can be assumed that deformability properties (Young’s modulus and Poisson’s ratio ) of the rockmass may be approximately the same as these of the intact rock. For the conducted FEM analyses, the modulus and Poisson’s ratio values of the intact LdB granite, as listed in Table 1, were used directly into RS2.
However, for the FDEM model these values cannot be used directly, as the overall stiffness of the system is affected by the employed penalty coefficients [36]. Therefore, the selected modulus and Poisson’s ratio values need to be adjusted properly. In order to overcome this, a numerical model of an unconfined compressive strength (UCS) test is used in order to calibrate the required deformability parameters and penalty coefficients so that the modulus and Poisson’s ratio of the LdB granite can be obtained. In Figures 3 and 4, the UCS model configuration and the obtained stress-strain curves are illustrated, respectively. A viscous damping factor is introduced into the model in order to minimize numerical oscillations in the linear part of the stress-strain curves. As observed, before any fracturing occurs, the specimen has a linear-elastic response with the model behaviour governed by the 3-node elastic triangular elements, as long as the strength of the cohesive elements is not exceeded. The established parameters are listed in Table 5.
(a)
(b)
5.3.2. Strength Parameters
For the FEM models, there are two major assumptions in order to simulate the strength of the rockmass. The first assumption is that the rockmass is going to behave according to the Hoek–Brown criterion for a geological strength index (GSI) [23] value of 80 (massive rockmass with no or sparse fractures and moderate to good quality discontinuity surfaces) and a UCS value equal to 200 MPa, which falls within the range of the UCS of the LdB granite, as listed in Table 1. The second assumption is that the rockmass is expected to behave according to the DISL model, as it was described in Section 2, with the UCS of the intact rock being in this case as well 200 MPa, and the rest of the parameters of the modified Hoek–Brown adopted from [16]. The complete set of parameters used is shown in Table 5.
For the FDEM models, once the deformability parameters have been established, the strength properties of the cohesive elements need to be determined in order to replicate the field conditions and failure mechanisms observed in field. It has to be noted that for the FDEM method, no constitutive model is required, and the fracturing occurring depends solely on the strength parameters of the cohesive elements. In order to achieve an agreement between field and model observations, a trial-and-error calibration process was employed. In Figure 5, the formation of the “v-shaped” notch observed in situ and the spalling processes simulated in the FDEM model are illustrated. The established strength parameters of the cohesive elements are listed in Table 5. In Figure 6, a complete flowchart of the calibration methodology employed is demonstrated based on [46].
(a)
(b)
5.3.3. Discrete Fracture Networks
For the purposes of this study, two types of models were developed within RS2 and Irazu. The first type of model was created in order to simulate the fracture conditions encountered at the URL Test Tunnel; hence, no fractures are present [43]. The second type of model was assigned a discrete fracture network comprised of one subvertical (mean dip 80°) and one subhorizontal (mean dip 10°) joint sets. The number of the generated fractures is determined by determining two additional parameters: (a) the areal fracture intensity , defined as the sum of fracture trace lengths divided by the mapping area [47] and (b) the mean trace length of the fracture traces. The input parameters for the generation of the employed DFN geometries used for both the FEM and FDEM models are listed in Table 6. The generated DFN geometries were installed within the high-resolution areas of all the models, as described in Section 5.1. Additionally, the fractures are assumed to have a purely frictional behaviour in order to simplify the conducted analyses and the interpretation of the obtained results. It has to be noted that these values have been selected arbitrarily in order to investigate the effect of the rockmass structure within the FEM and FDEM models, on the development of the fracturing mechanisms, and the extent of the damage based on each numerical method.
6. Results and Discussion
6.1. Intact Models
In order to assess the impact of the selection of the numerical method on the rockmass response and whether or not it captures the in situ rockmass behaviour, based on the field observations from the URL Test Tunnel, initially the extent of the potential damage on the rockmass due to the excavation is examined for the models that do not have a DFN assigned to them. In the FEM models, the extent of the damage is initially assessed by identifying the yielded elements around the excavation circumference. Regarding the FDEM model, the numerical results are assessed based on the material that collapses as a result of the fracturing processes taking place. The yielded elements for the FEM models are illustrated in Figures 7 and 8. The collapsed material within the FDEM model is highlighted in Figure 9.
(a)
(b)
(a)
(b)
(a)
(b)
In Figures (7a), 8(a), and 9(a), results obtained for the same amount of deconfinement (40% decrease of the initial stress state) show that all numerical models still have a linear-elastic behaviour, as no fracturing has occurred yet. Based on the major principal stress contours, the magnitude of the monitored stresses is approximately the same, as expected, as the models are still in a prepeak state. However, as further destressing is taking place with the advancement of the face, yielding and fracturing occur. By recording the extent of damage based on the criteria above, in Figure 10, the damage profiles are compared for the fracture-free models between one another and the actual damage profile from the URL Test Tunnel. As it can be observed, the damage profile resulting from the FEM model using the DISL approach resembles the damage profile obtained from the FDEM model. Furthermore, both of them show that they are able to provide a reasonable estimate of the highly damaged area. However, the model using the conventional Hoek–Brown criterion is not able to capture the failure mechanism in a realistic manner.
As discussed in the previous sections, the excavation response of a brittle rockmass cannot be captured by conventional shear failure criteria since the failure mechanism is not controlled by the shear strength of the rockmass. More specifically, in Figure 10, it is shown that at the locations in which field observations suggest material collapse, the numerical model using the Hoek–Brown criterion shows the lowest extent of damage based on the recorded yielded elements. Furthermore, the shape of the damage extent is completely different from the estimations of the other two numerical models and the collapsed material monitored in situ. More specifically, it resembles more a “butterfly” shape than the characteristic “v-shaped notch.” This is attributed to the application of the Hoek–Brown criterion which essentially predicts a higher rockmass strength than that of the actual rockmass strength in situ. This higher strength allows the rockmass to withstand the high compressive stresses that manifest at the crown and floor of the excavation due to the tunnel advancement. This results in zero stress release at these locations, and it promotes failure at the sidewalls of the excavation which are under a tensile regime.
On the contrary, the DISL model efficiently captures the shape of the highly damaged zone (HDZ) and is able to simulate the brittle response of the rockmass based on the modified parameters of the Hoek–Brown criterion. Furthermore, the model predicts tensile failure at the sidewalls of the excavation. While in situ, no tensile cracks were observed, acoustic emission events [49] suggest tensile microcracking around the area that the elements are yielding (Figures 8 and 10). A similar observation can be done within the FDEM model as shown in Figures 9 and 10.
The DISL approach is able to capture the brittle response of the rockmass based on the instantaneous cohesion-weakening frictional strengthening model with the Hoek–Brown parameters. This method is relatively simpler than other approaches such as the one proposed by Hajiabdolmajid et al. [27] in which the cohesive strength reduction and frictional strength mobilization reach their post-peak values depending on the magnitude of the accumulated plastic strain. However, the DISL approach does not come without some limitations. As observed in Figure 10, while the FEM model using the DISL approach captures the highly damaged zone (i.e., the collapsed material as a result of the brittle fracturing), the rockmass beyond that area behaves elastically as no more elements are yielding past that extent. Microseismic events (MS) [49], however, suggest that stress-induced damage because of the excavation is extending beyond the highly damaged zone. In FEM models, the elements are assigned a constitutive model in order to simulate the material failure. As a result of this, the HDZ is captured quite well with the DISL approach which provides a good constitutive assumption for simulating brittle fracturing based on the yielding of the elements. However, beyond that yielded material area, the yielded elements cannot be used as a precursor of the damage, and other measured quantities (volumetric strain and principal stress concentrations) need to be assessed in order to evaluate the damage extent [3] and the captured fracturing processes. As observed in Figure 11, a distinction has been made for the FEM-DISL model between the HDZ and EDZ based on the major principal stress. By comparing the results with the recorded MS events, it is shown that the potential EDZ extent is underestimated as the proper mechanics of stress redistribution due to fracturing are not adequately captured. The FDEM approach, on the other hand, allows for the direct distinction between HDZ and EDZ, and it is more capable of capturing fracturing microprocesses as it allows for fracture propagation and interaction while the overall material maintains its structural integrity. As shown in Figures 9 and 10, this excavation-damaged zone (EDZ) can be captured by the calibrated FDEM model, hence overcoming the limitations of the FEM-DISL model. Furthermore, as observed in Figure 9, the collapsed material fails as a result of extensile fracturing (red lines) due to exceeding the tensile strength of the rock close to the excavation boundary where the confinement is low. Therefore, the developed FDEM model and the FDEM method show that are mechanisms capable of simulating the actual failure in order to replicate the in situ conditions.
The accurate simulation of the physical phenomena taking place and resulting in brittle fracturing under high stresses is of great importance, as they affect the design process and the selection of the appropriate temporary support measures. In Figure 12, a monitoring line starting at the crown of the excavation and going into the rockmass for a length of 2.25 m is shown. In Figures 11 and 13, stress measurements taken along this line highlight the difference between each model and field stress measurements conducted at the URL Test Tunnel (SM-5 stress cell) [50] and their subsequent implications in the design process. It can be observed that despite the fact that the in situ measured stresses are not in a very good agreement with the numerical results in terms of the absolute values and location, the FDEM model is capturing the possible stress conditions at the URL Test Tunnel. More specifically, in Figures 11 and 13, the principal stresses and which were recorded from the numerical models are higher than the stresses recorded at the SM-5 stress cell of the URL Test Tunnel. However, that stress measurement point was located in the disturbed zone around the excavation. The numerical results from the FDEM model show that the highly damaged zone along the monitored line extends up to 25 cm from the excavation boundary (blue dashed line). Within that range, the stresses monitored are zero, since the material is collapsing. Beyond that distance (blue box), however, a great fluctuation in the recorded stresses can be observed as the monitoring line now lies within the excavation damaged zone. The occurred fracturing results in a nonsmooth stress distribution along the line, with some points close to the SM-5 stress cell recording similar stress magnitude values. Of course, the exact same location is not the one monitored, a factor that also contributes in not obtaining the exact same value. On the contrary, from the FEM models, a distinction between the HDZ and the EDZ cannot be made instantaneously based on the extent of the yielded elements. In this case, as previously mentioned, the principal stress components were examined in order to determine the extent of the HDZ and EDZ, as suggested in [3]. The model using the Hoek–Brown criterion predicts a significantly smaller area affected by the excavation, along the monitoring line. Furthermore, the stresses beyond that distance correspond to these obtained from the material with a linear-elastic behaviour; from acoustic emission events [49], it is clear that the rockmass is damaged even beyond the collapsed material. For the FEM model in which the DISL approach was applied, according to the major principal stress measurements (change in the slope of the curve in Figure 11), it is shown that the estimated HDZ extent (red dashed line) is approximately the same as the one measured in the FDEM model (blue dashed line). However, the estimated extent of the EDZ is less than the one predicted from the FDEM model, which is more in agreement with the field observations. Furthermore, it fails to replicate the stress conditions which are monitored by the SM-5 stress cell. Therefore, the conducted stress measurements clearly show the merits of the creation of a more advanced numerical model in order to simulate such rockmass and stress conditions and promote a more efficient support design.
6.2. Fractured Models
In the previous section, the FEM and FDEM models were used in order to replicate the conditions of the URL Test Tunnel for the virtually “fracture-free” LdB granite, based on the use of different constitutive assumptions for the FEM models and the employment of the FDEM method. By assuming that the material parameters used in the previous section were representative of the intact rockmass, the same models were modified in order to include structure explicitly simulated within them. In order to achieve this, DFN geometry was assigned to each of the numerical models. Because the internal DFN generators incorporated within the Irazu and RS2 codes were used, while the initial DFN parameters were the same for both (Table 6), the exact same geometry was not able to be created for both codes. However, after generating a number of different geometries, two were selected, one in Irazu and one in RS2, with similar joint patterns in order to obtain comparable results.
In Figures 14–16 the high-resolution area of the numerical models containing the DFN geometries are illustrated. In the FEM model cases (Figures 14 and 15), the yielding elements, as in the previous section, are used in order to assess the stress-induced damage due to the tunnel advancement. Within the FDEM model (Figure 16), the excavation damage zone is divided into two separate groups including (a) the material collapsing because of the spalling processes (HDZ) and (b) the damaged rockmass that still maintains its structural integrity due to the confining stresses (EDZ). In all three different cases, it becomes evident that the explicit simulation of joints within the numerical models affects the stress-induced damage. These discontinuity elements act both as stress barriers and stress concentrators, and the length and orientation of the rockmass joints control both the shape and extent of the stress-induced damage.
More specifically, in Figure 17, the different damage profiles are compared to one another, and they are also compared to the excavation damage profiles obtained from the intact models, as discussed in the previous section. As it can be seen, the fractured models are strongly influenced by the explicit integration of the DFN geometries within the numerical models in all three cases. In the case of the FEM models (Figures 17(a) and 17(b)) which have the same fracture pattern, it can be observed that similar damage profiles are obtained, both in terms of magnitude and shape. The joints highlighted in Figures 14 and 15 control the inflicted damage to the rockmass (yielded elements) around the excavation. However, the use of a different constitutive approach in these two models results in a larger damage extent at the crown and floor of the excavation when the DISL approach is employed. On the contrary, when the conventional Hoek–Brown criterion is employed, the crown and floor damage appear to be minimal. In Figures 14 and 15, it can be seen that at the crown and floor of the tunnel, two relatively big areas are fracture free. Therefore, the rockmass behaviour in these two particular regions is controlled by the presence of these rock bridges (intact parts of rock), and the selection of the constitutive model for the rockmass simulation affects how the material is going to fail. More particularly, the DISL approach is promoting the brittle rockmass failure, as in the case of the intact model. On the contrary, the Hoek–Brown criterion overestimates the rockmass strength resulting in minimal damage at the crown and floor of the excavation.
(a)
(b)
(c)
As discussed in the previous section, the FDEM method can capture the behaviour of brittle rockmasses by simulating the extensile fracturing that is taking place, once the tensile strength of the intact rockmass is exceeded, and hence forming fractures that propagate along the direction of the major principal stresses . As a result of this, the damage profiles between the FDEM and FEM-DISL models have similarities due to their capability of simulating the brittle response of the excavation as the face advances. As observed in Figures 17(b) and 17(c), the damage inflicted on the rockmass at the crown and the floor of the excavation due to brittle failure, where intact parts of rock are present, is captured in both models. However, in the FDEM model, a clearer distinction between the HDZ and the EDZ can be made as result of the complete simulation of fracture initiation and growth.
Following the comparison between the fractured models with one another, each intact model is also compared to its corresponding fractured model in order to investigate both the impact of the constitutive model and the employed fracture pattern. In sparsely to moderately jointed rockmasses, the material behaviour is governed by the existence of intact rock parts between the joints present. In this study, joint structure was explicitly simulated, and therefore the behaviour of the medium and the damage extent due to the excavation are controlled by both the constitutive assumptions of the intact rock (rock bridges) and the overall geometry of the joint system. In Figure 17(a), it can be observed that the applied DFN geometry influences the damage profile at some extent, especially in areas that the fracture intensity is higher. However, in areas where the influence of the structure is not that significant (excavation crown and floor), the damage extent is similar in both models due to the higher strength predicted by the Hoek–Brown criterion. The overall shape of the damage profiles though is not significantly different between the two analyses as the overall rockmass behaviour is dictated by the employed constitutive model that promotes yielding towards specific directions because of the rockmass shear strength, hence decreasing the impact of the preexistent joints.
In the FEM-DISL model, however, that is not the case. As observed in Figure 17(b), by taking into account the brittle response of the intact material, through the application of the DISL approach, results in significantly different damage profile for the intact and fractured model. Since the constitutive model of the intact rock promotes brittle failure, this leads to the accumulation of damage at the intact rock parts located at the roof and floor of the tunnel in both cases. However, the presence of structure within the rockmass in the fractured model governs the propagation of the yielded elements. The subhorizontal and subvertical joints close to the excavation act as stress barriers and stress elevators at specific locations, hence shaping the formed damaged zone. This becomes more evident at the sidewalls of the excavation. In the intact model, the high horizontal stress leads to a high bending moment at the sidewalls, hence resulting in tensile failure due to bending. However, in the fractured model, that is not happening. The rockmass structure leads to a redistribution of the induced stresses in a way that this bending effect is not influencing significantly the excavation. On the contrary, it results in accumulation of damage at the upper and lower half of the tunnel.
A similar observation can be made for the FDEM model (Figure 17(c)). The effect of the preexistent joints in this case becomes evident on both the collapsed material (HDZ) and the damaged material (EDZ). As seen, the collapsed material in the case of the intact model (blue dotted line) is dictated by the applied anisotropic stress regime, with the formation of the notch following the general direction of the geostatic minor principal stress . On the contrary, for the fractured model, the notch formation (black continuous line) is influenced by the occurring extensile fracturing but at the crown is also controlled by a strong subhorizontal joint. This joint redirects the induced stresses and makes the notch to form to the right. That is not the case though for the floor of the excavation (black continuous line). In both models, a large area of intact rock can be seen. Therefore, the stress-induced fracturing is mainly the result of spalling, hence leading to similar HDZ profiles for the excavation floor. Regarding the EDZ, the influence of the DFN geometry is stronger. In the intact model, the EDZ follows the general direction of the HDZ. However, the damaged material is contained at the upper and lower parts around the excavation (blue dash-dot line), with minor tensile fracturing occurring at the sidewalls (blue dotted line) because of bending. In the fractured model, the preexistent joints result in a wider area around the excavation to be damaged (black dashed line). In Figure 16, the rockmass structure controlling the EDZ (white lines) creates specific blocks of intact material that restrain the stress-induced cracks. More specifically, strong subhorizontal joints at the crown guide fracture propagation to the right of the excavation, and the subvertical joints at this location act as constraints for the propagating cracks. For the excavation floor, fracturing is intensified due to the presence of some joints, but the damage of the material is mainly controlled by spalling due to the presence of rock bridges. Finally, the sidewalls are not imposed to significant damage due to the stress redistribution. However, it is evident that preexistent subvertical joints interact with one another resulting in the coalescence of adjacent discontinuities.
7. Conclusions
The numerical modelling of hard rockmasses has been attempted by employing different numerical techniques, including the use of continuum- and discontinuum-based approaches. In the present study, the numerical codes RS2 and Irazu were used in order to perform the numerical simulations using the FEM and the FDEM methods, respectively. The FEM models simulating the fracture-free rockmass were based on data obtained from the URL Test Tunnel for the LdB granite and relevant research that has been done on that specific site. Furthermore, the intact FDEM model was calibrated in order to replicate the failure mechanism observed at the URL Test Tunnel.
Once the intact numerical models were established, DFN geometries were integrated into them in order to investigate the impact of joints on the rockmass behaviour during an excavation when different constitutive assumptions are made and the joints are explicitly simulated. For the FEM models using the conventional Hoek–Brown criterion, for a GSI = 80, it was shown that the estimated damage profile does not correspond to the field observations for the intact model. Furthermore, the addition of joints does change the shape and the extent of the damage profile; however, it does not result in significant changes to the overall rockmass behaviour. The rockmass constitutive model in this case appears to be the main contributing factor in the rockmass response. However, this does not lead to a realistic damage profile. The FEM models using the DISL approach have the capability of capturing the brittle response of the rockmass when a continuum technique is employed. For the intact model, the URL field observations are replicated to an extent, and the simulation of the physical processes taking place is more accurate than using the conventional Hoek–Brown criterion. However, HDZ and EDZ cannot be distinguished directly and additional processing is required. On the other hand, the intact FDEM model is capable of simulating the extensile fracturing occurring during spalling, and a clear distinction between the HDZ and the EDZ can be made. While both the FEM-DISL model and the FDEM model provide similar predictions of the HDZ, the EDZ predicted by the FDEM model is in a better agreement with field observations, as the numerical results are more consistent with AE and MS data from the URL Test Tunnel. Additionally, in situ stress measurements correspond better with stress measurements from the FDEM model.
Once joints are added in the FEM-DISL and the FDEM models, it can be observed that similar damage profiles are obtained, since both the effect of the intact rock and the rockmass structure can be adequately captured. However, again the distinction between the HDZ and the EDZ is clearer in the FDEM model where the effect of preexistent joints on the rockmass response is more profound. This highlights the merits of such an advanced numerical technique which has the potential of assisting greatly in the design of an efficient support system for such rockmasses, as the stress and damage state of the material surrounding an excavation can be more realistically captured.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The authors would like to thank the Nuclear Waste Management Organization of Canada (NWMO), the Natural Sciences and Engineering Research Council of Canada (NSERC), the Ministry of National Defense, and the RMC Green Team for providing the founding and the resources for this work. Additionally, the authors would like to thank Dr. Omid K. Mahabadi and Dr. Andrea Lisjak (Geomechanica Inc.) for their meaningful discussions, assistance, and support in the preparation of this paper.