3D FEM and DEM Analyses of Underground Openings in Competent Rock Masses

The paper is aimed at comparing the results of numerical analyses of underground openings in competent rock masses like the Carrara Marble (Italy) by considering a real and well documented case study. More specifically, 3D FEM and DEM analyses were carried out on a rock-mass model interested by two faults and three sets of discontinuities. The geometrical model is representative of deep underground openings where spalling-cracks and rock bursts can occur. PLAXIS 3D and 3DEC were used for the analyses. Intact rock and rock mass characterization of Carrara Marble was inferred from available technical literature. The analysis results were compared in terms of principal stresses and displacements in a number of monitoring points around the opening. The main practical interest is to find out a reliable approach for evaluating the stability of very large openings in a competent rock mass like Carrara marble. For such a purpose, a number of available in-situ stress measurements were used.


Introduction
: Schematic diagram suggesting the range of application of dis-continuum modeling in relation to Q value ([9] modified by [10]). (for the study case Q is about 7; RMR = 9*LNQ+44) Nowadays the continuum approach is still the most popular and few Authors have tried to compare the results of these two types of analyses in real cases [8][9][10][11]. On the other hand, in the case of layered structures (i.e. shale formations), where the rockmass behaviour is mainly controlled by the set of discontinuities, it could be convenient to use the continuum approach. Therefore, parametric studies have shown in which way the anisotropic behaviour of these types of rock -masses can be well simulated by using the continuum approach [12,13]. The case study does not belong to this category because it concerns a competent rock mass with few discontinuity sets. This paper deals with a real and well documented case study. According to Scavia [14] we believe that equivalent continuum and dis-continuum approaches may lead to very different results. Therefore, these approaches should be used for different purposes and objectives. Our aim is to compare, in terms of principal stresses and displacement vectors, the results of 3D continuous and discontinuous analyses that were carried out by using PLAXIS 3D [15]) and 3DEC [16]. The study model consisted of an "idealized" rock mass block (700x400x595m) interested by two faults and three sets of discontinuities. These structural features correspond to those observed around an existing deep quarry, as better explained later. The opening is 49x50x30m and the floor is 95m above the block-bottom ( Figure 2). The geo-mechanical characterization of the rock mass, intact rock and discontinuities was inferred from the technical literature. Carrara marble is extracted from the mining district of Carrara in the North-West Tuscany (Italy). It is a generic term indicating different commercial products of a carbonate metamorphic rock.
More specifically the so-called "Bianco di Carrara" is considered in the present study. The mining district is actually very wide and the mine exploitation is carried out in different ways. The paper firstly summarizes the mechanical characterization and classification indexes of Carrara marble from the available technical

Geo-mechanical characterization of the Carrara marble
The geo-mechanical characterization of Bianco di Carrara is based on a number of research papers [17][18][19][20]. Table 1 & 2 summarize the geo-mechanical characterization and the classification indexes of "Bianco di Carrara" as inferred from literature [17][18][19][20]. The strength and stiffness parameters specifically used for the numerical analyses are summarized later on. In these tables, values referred to x, y, z directions were obtained from laboratory testing on groups of specimens drawn from three perpendicular directions, of which z is directed orthogonal to the apparent marble layering [20]. Shear strength of intact rock (τ) was evaluated by means of direct shear tests on intact specimens that were subjected to a normal stress σ a =3.5 MPa.

ACET.000588. 4(3).2020
Tables 3-5 show the geo-mechanical characterizations of discontinuities as inferred from existing technical literature. Table   3 shows the characteristics of discontinuities surveyed in the Ravaccione & Fantiscritti pit [18]. The mechanical characterization of these discontinuities was obtained by laboratory tests, carried out according to ISRM suggested methods [18]. Table 4 concerns discontinuities surveyed in another quarry [19]. Table 5 refers to artificial discontinuities that have been prepared and tested in the lab [20]. It is worth mentioning that in Table 3 the reported value of cohesion was inferred from a specific analysis trying to account for the existence of rock bridges. The cohesion value of Table 5 was instead experimentally determined in the lab on artificial discontinuities.

Numerical methods
This paper aims to compare FEM and DEM methods in order to outline limitations and capabilities of these different approaches, in the specific field of underground excavations. The comparison concerns a 3D model of a blocky rock mass containing a cavity.
The case study represents an idealization/simplification of the Ravaccione Fantiscritti Quarry. Indeed, the sets of discontinuities correspond to those effectively observed [18], while the geomechanical characterization was inferred from those reported in Table 1-5.
As for the numerical analyses the adopted geo-mechanical characteristics are reported in Table 6  (K=bulk modulus, G=shear modulus).   (1)

ACET.000588. 4(3).2020
where τ and σ' n are the tangential and normal stresses on the failure plane at failure, respectively, c is the cohesion of the intact rock and Φ is the angle of shear resistance.
Mechanical parameters (Table 6) adopted for the rock matrix were inferred from those reported on Table 1 According to the Bienawski rating system [23], the rock mass is characterized by a geological strength index GSI=66 [18].
It is worth mentioning that, the set of parameters for DEM and FEM analyses were indepentely inferred fron available literature data. Consequently a different σ c was considered. In particular the uniaxial compression strength was considered equal to 80MPa for the DEM analyses while was equal to 99MPa in FEM analyses. In other words we would approach the analyses such as a practizing engineer using different computational tools that require different sets of input parameters. Anyway, the adopted MC (intact rock) and HB (GSI=100) criteria can be considered equivalent as better shown later on.

Model geometry and monitoring points
The model used for numerical analysis represents a rock mass block (700x400x595m) containing an opening/cavity (49x50x30m). The opening's sides are parallel to the block's one. Figure 2 shows coordinates (x, y, z) respectively of the whole block (at the left) and of the opening (at the right). The block is crossed by two faults and by three sets of discontinuities (K1, K2, K3). Table 8 shows orientation and mechanical characteristics of the faults. The  In order to perform finite element calculation, the mesh was established after a convergence study. Plaxis uses an automatic mesh generation process which considers all model's singularities.
Moreover, local refinement has been performed close to the cavity. A similar approach was used for DEM analyses, with the discretization of each block to the finite differences.

Results
The outputs of the numerical analyses were compared in terms of principal stresses and displacements. Figure 4 shows The linear Mohr-Coulomb envelope was cut at σ3 = -8MPa, that is the tensile strength of Carrara marble ( Table 1). As for the HB criterion, both the case for GSI = 100 (intact rock) and that for GSI = 66 are shown in Figure 4. The MC criterion (intact rock) and that by HB (GSI=100) mainly coincide for the relevant stress interval. In the following, stresses and cohesion are in MPa even if not explicitly indicated.    Table 10 shows the maximum displacements obtained by PLAXIS in the three directions x, y, z. Table 11 shows the values of the maximum principal stresses obtained by the two programs on the whole model with the relative point's coordinates. Graphic outputs of 3DEC show a max displacement of about 40cm, in the section parallel to the xz plane and passing through the center of the chamber, at the intersection between the chamber itself and the K3 discontinuity set. Still in the vicinity of the room, but outside the intersection with the K3 set, displacements vary between 1 and 2cm and are therefore comparable with those obtained from PLAXIS.     In this case the value of σ c was kept equal to 99MPa. The old and the new set of parameters are shown in Table 12. Numerical analyses were repeated by considering the new sets of parameters and the results are shown in Figure 8 that also shows the four strength envelopes. It is possible to see that DEM, in any case, gives higher stresses than that predicted by FEM analyses. Anyway, we performed an analysis by considering the above strength criterion. The results are summarized in Figure 9. The same figure also shows the stress state measurements (safe condition) carried out by Gulli [25].