Stability analysis and support system design of Angouran underground mine batching room

Purpose. The aim of this paper is to develop a numerical model to investigate the stability of Angouran underground mine batching room, in which appropriate support systems are selected, designed, and compared to choose the most optimal case. Methods. Stability analysis is carried out by critical shear strain criterion using finite element method (FEM). The effect of support system on the stability of batching room has been determined for two types of composite support systems: one, incorporating rock bolt and shotcrete, and the other, comprising steel frame with one-meter distancing and shotcrete, based on comparing the displacements around stopes with the amount of critical and allowable strains resulted from Sakurai criterion. Findings. The results of the analysis show that the excavation of the batching room without a support system causes its instability and failure, so it is imperative to design a suitable support system according to the geomechanical and stress conditions of the area. Comparison between the two selected support systems also shows that there is not much difference between these two support systems and in both cases, batching room stability can be provided. Originality. A wide range of parameters involved in the calculation were explored. Results suggest that the failure strain seems to be a promising tool when used correctly. Considerable time saving can be achieved without losing accuracy. Practical implications. Numerical methods can be utilized to analyze tunnels with non-circular cross-sections, considering heterometric ground conditions around the tunnel, often heterogeneous in situ stresses, rock mass behavior and tunnel face effects. By selecting suitable rock failure criteria, it is possible to identify weaknesses around the excavations via numerical modeling and design an appropriate support system based on it.


Introduction
Issues such as the stability of large underground spaces, post-excavation deformation and stress analysis and excavation methods optimization are of particular importance in the underground mining and excavations [1]. The stability of underground excavations depends on their dimensions and geometry, excavation method, number and sequence of excavation stages, in situ stress conditions, support system and its installation time [2], [3]. With the development of computers, modeling methods and engineering knowledge, including the continuous and discontinuous environments mechanics, analytical and simulation methods have flourished rapidly, so that today, with the availability of high-speed computers and professional software, it can be claimed that any underground structure with any characteristics can be modeled and analyzed for pre-implementation stability or support system applicability [4], [5].
In numerical methods, the most essential factor in obtaining accurate results is the entry of correct information. Determining the behavior of the rock mass (stress-strain relationship), its geomechanical characteristics, and modeling external factors that are probably to cause instability is very important and at the same time, difficult [6]. For this reason, numerical analysis is a complex process whose results vary with changes in boundary conditions, model structure, rock mass properties, ambient temperature, groundwater conditions, and other factors that affect stability. In this regard, to model the engineering behavior of rock mass, a sub-model should be created to select the best behavioral model for rock mass by changes in the model and problem input. The input information of the problem should also be determined via geological surveying, rock mechanics laboratory and field experiments, and back analysis [7]. Numerical methods can be utilized to analyze tunnels by the non-circular crosssections, heterometric ground conditions around the tunnel, often heterogeneous in situ stresses, rock mass behavior and tunnel face effects, and three-dimensional tunnel behavior [8], [9]. By selecting suitable rock failure criteria, it is possible to identify weaknesses around the excavations via numerical modeling and design an appropriate support system based on it [10]. Ghafouri et al. (2010), utilized the experimental methods of RMR, Q, and GSI as well as the numerical method distinct elements, to design the temporary support system of the Darungar dam water transfer tunnel located in Razavi Khorasan province, northeastern of Iran, incorporating rock bolt, shotcrete and iron frames in four different tunnel segments. In the study, by changing the experimental to the numerical method, the accuracy and reliability of the design were improved. The support system of the water transfer tunnels' intersection with shock absorber tanks of Gotvand dam located in Khuzestan province, southwestern Iran, was designed utilizing experimental methods and finite difference method and showed that after the installation of the support system, the maximum displacement is less than the Sakurai's critical displacement [11]. Flex et al. (2008) modeled the tunnel excavation machine (TBM) with a compressed air shield based on the finite element method (FEM). In this study, loose soil around the tunnel is considered as a threephase material that consists of soil material, water in the pores, and air. The tunneling process is Three-dimensionally modeled by applying the effect of compressed air as a support factor of the tunneling face along with the installation of the segment.
The modeling results can be used in estimating the suitable tunnel face pressure to ensure the stability of the tunnel face and landslide control [12]. Bakhshandeh et al. (2014) utilized the convergence-confinement method and Hoek-Brown failure criteria to analyze the stability of the long tunnel support system of Sardasht dam via numerical software based on finite difference method (FLAC 3D) and thus showed the applicability of Hoek-Brown failure criteria [13]. Torkmanju et al. (2018) analyzed the stability of the tunnel face in mechanized excavation and heterogeneous conditions utilizing the resistance reduction method and the concept of safety factor via earth pressure balance method (EPB) with the aid of numerical modeling incorporating Plaxis 3D Tunnel and Phase2 software [14]. In another study, the stability of tunnels at the intersections of Y-shaped cross-sections considering their impact angle has been studied. The results of this numerical modeling show that by decreasing the impact angle, the additional support length required in the middle wedge and the crown of the tunnels increases [15]. Sabzevari et al. (2017) have performed stability analysis and design of water tunnel support system on Mashoureh dam utilizing experimental, analytical, and numerical methods. In the experimental method, two methods, Q and RMR are utilized and in numerical analysis 2-phase software, which is based on finite elements is used, and finally, the results of aforementioned methods presented a similar support system [16]. Basirat et al. (2018) have studied the numerical analysis of the segmental support system used in most tunnels bored with TBM and the effect of the earthquake on the resulted vertical and shear forces at the contact surface between the various parts. In this research, utilizing the information of Karaj water transfer tunnel to Tehran (Amirkabir) and the 2D UDEC software, the segmental coverage under two conditions of complete slip and non-slip, has been investigated and the shear and normal forces of the joints between the two segments have been analyzed [17]. Deformation and subsidence of the ground surface due to tunneling operations as well as displacement of the crown, walls, and floor of the tunnel are amongst the unstable factors of the tunnel. Ghiyasi and Kooshaki (2019) using numerical modeling have studied the subsidence of the earth in circular tunneling due to simultaneous changes in the geometric characteristics of the tunnel and the mechanical properties of the soil in the saturated environment. In this study, utilizing the finite element analysis method in ABAQUS software, considering the different values of tunnel depth and diameter in different conditions for soil adhesion, internal friction angles, permeability coefficients and their simultaneous changes, the effect of these variables on land subsidence has been investigated [18]. Nejad Shah Mohammad (2019) has evaluated and optimized the segmentation of these tunnels via utilization of the separate element method and simulating the body of large crosssection tunnels to minimize displacements at the ground level and reduce the expansion of shear strains. The variable utilized in this analysis is the ratio of the height of the progressive section to the bench height in the two-stage excavation method. Two twin D-shaped and circular tunnels with an average diameter of 13.7 meters, which are located at a distance of 33 meters apart, have been selected for case studies [19]. In another study by Obidor and Giant, the stability analysis of single and twin horseshoe tunnels in a rocky environment that is affected by additional pressure from the ground has been performed. By assuming that tunnels are long, planar strain laws have been utilized in the modeling. The Hoek-Brown fracture criteria are used to model rock mass fracture. To study the behavior of rock mass under different surface loading conditions, finite element numerical modeling has been utilized [20]. Trong and Jian investigated the interaction of twin tunnels in a joint rock mass using numerical modeling and discontinuity analysis (DDA). Land subsidence and redistribution of rock stress around tunnels are obtained from DDA simulations. The results show that the amount of tunnel interaction is significantly affected by the distance between the two tunnels of the joint rock mass. In addition, the slope angle of the rock layers significantly affects the patterns of subsidence and redistribution of stress around the excavation site [21]. In another study, Mollai et al. evaluated the effect of changing the material and slope of soil layers around the tunnel cross-section. The results of numerical modeling show that horizontal layering causes changes in displacement rate and if the soil in the upper floors of the tunnel has more resistance than the soil in the lower floors, the layered soil Pannet curve will be lower than the analogous soil curve. According to the results obtained for tunnel design in sloping layered conditions, the use of a modified Pannet curve is recommended [22].
Angouran lead and zinc mine is located about 125 km southwest of Zanjan city and according to the exploitation plan, at least 120000 tons of ore will be extracted annually from the underground section of the Angouran mine. To obtain the maximum amount of mineral (maximum grade and recovery) and taking into account the geometric properties of the deposit, geomechanical properties of the mineral, the surrounding rock, etc., the suggested suitable extraction method for this mine is stope and strip pillar with downward filling. By examining and selecting the type of filler to make a durable and integrated artificial hanging wall and the conditions for utilizing other filling techniques, the filling method of choice for this mine is concrete filling. To fill the stopes in the proposed design and ensure that no subsidence occurs during mineral extraction, a cement filler mass with a strength of at least 15 MPa (150 kg/m 3 ) is considered which is imperative to obtain the best mixture via utilizing and testing the available materials, various types of cement, stone materials and water in the workshop [23]. By examining various aspects of implementation, topographic conditions, mining geometry, and available equipment, it was suggested that to transport the filler and concreting the stopes at the end of the tunnel N. 2, two batching rooms should be designed and implemented. The dimensions of the batching rooms are designed to have a width of about 7 meters and a height of about 8 meters, the position of which to the preparatory excavations and stopes is shown in Figure 1.

Figure 1. The position of the batching room to the preparatory excavations and extraction stopes
The purpose of this study is to analyze the stability and design of the support system of the batching rooms based on the finite element method (FEM), in which appropriate support systems are selected, designed, and compared to choose the most optimal case (technically, economically and executive). Soil and rock behave non-linearly under load. This nonlinear stress-strain behavior can be simulated with several behavioral models. As the model becomes more complex, the number of input parameters of the analysis increases therefore, laboratory and field experiments are required. To model the elastoplastic behavior of soil and rock materials, the Mohr-Coulomb behavior model has been employed. The Mohr-Coulomb rupture criteria is a nonlinear, robust, and simple criterion and this soil/rock behavior model is defined by the five main parameters of elasticity modulus, Poisson's ratio, adhesion, friction angle, and expansion angle.

Geomechanical characteristics of Angouran lead and zinc rock masses
Angouran lead and zinc deposit is one of the world's leading deposits in terms of quality and consists of three parts: oxide, sulfur, and a mixed part. The oxide part is located in the highest part of the deposit, the sulfur part is located in the lowest part and the mixed part is located between the two [23]. The three-dimensional view and lithology of the Angouran deposit are shown in Figure 2.

Figure 2. 3D view and lithology of Angouran deposit
In general, the mineral and rock masses in the underground section of Angouran lead and zinc mine are classified in class II (medium to high-quality rock masses: 50 ≤ RMR < 60), class III (medium quality rock masses to low: 35MRR < 50) and class IV (poor quality rock masses: RMR < 35). The geomechanical parameters used in the modeling based on the results of drilling cores including RQD (rock quality index), GSI (geological strength index), compressive strength, and joint properties are shown in Table 1.

Numerical modeling and support system design
The primary purpose of installing a support system is to control the displacement and convergence of the walls of the excavated space and to prevent loosening of the rocks around the cavity, which may lead to the collapse of the excavated space. On the other hand, in the distance from without support to the tunnel face (excavation step), the ground releases some of its stresses; this reaction of the earth relative to the excavated space, depending on the type of mass-rock of the earth around the hole, is a function of time (earth movement). In modeling underground spaces, disregarding a proper ground stress release and closing the support system sooner or later will cause instability or excessive rigidity of the support system [24].

Convergence-confinement method
Since the installation of the support system is not implemented immediately after excavation operation and during this period a large amount of stress around the stope is released, to have a criterion for evaluating the stope during excavation and installation of the support system, the convergence-confinement method has been employed. First, the ground response curve (GRC), the radius of the plastic zone, and the maximum convergence of the tunnel are obtained according to the mentioned characteristics for the rock mass. According to the choice of the average drilling step of one meter, the degree of convergence at a one-meter distance from the tunnel face is also obtained from the longi-tudinal deformation profile of the tunnel (LDP). Then, in the GRC curve, according to the degree of convergence obtained, the amount of proportional stress release is determined. After applying the determined stress release, the support system which is simulated using a structural element (a combination of beam, concrete, and mesh) is installed. Afterward, the residual stress is released and finally, the whole model is mechanically analyzed and the results of the analysis are examined for the support system in the form of axial forcebending moment and axial force-shear force diagrams [24], [25]. The obtained GRC curve for the roof, floor, and wall of the stope along with its longitudinal displacement profile curve is shown in Figure 3.

Equivalent cross-section of the composite support system
In tunnels where the design of the support system is a combination of steel frames and shotcrete, it is necessary to design the components of the consolidation, the role and share of each of these support systems to be calculated separately. One of the methods to simulate the combined frame support system incorporating steel frame with shotcrete is the equivalent cross-sectional method, i.e., to simulate the support system, which consists of two different types of materials (steel frame and concrete) and with different crosssections, employing the relevant relations, an equivalent cross-section is considered and finally after calculating the loads subjected to the equivalent cross-section, the share of each steel frame and shotcrete support in the loads entering the support system is recalculated separately. In Figures 4 and 5 the cross-sections of each of the two different types of materials (steel frame and shotcrete) for the composite support system and its equivalent cross-section are shown.

Figure 7. Redistribution of loads over each of the cross-sections: (a) bending moment (M); (b) axial force (N); (c) shear force (Q)
The thickness (teq) and Young's modulus (Eeq) of the equivalent cross-section are calculated using the following equations: where: k1, k2, D1 and D2 values according to Young's modulus (Eeq), Poisson's coefficient (ϑ), cross-sectional area (A), and moment of inertia (I) of each of the cross-sections are obtained from the following relations.
Related to the axial force (N): The final stability of each component of the composite support system is obtained via assuming that the beam element for each of the structural support elements due to the relation between axial force and the following bending moment and based on the maximum stress created in the supporting component: where: Pthe axial force; Mthe bending moment created in each element (or their maximum value); Athe cross-section; Sthe cross-section of the element (beam). According to the presented materials, cross-section No. 1 is related to the steel frame and cross-section No. 2 is related to shotcrete. The characteristics of each cross-section (for example, for IPE160 frame, shotcrete, and 6 mm mesh) are given in Table 2.

Simulation assumptions and model geometry
The distance between the excavations and the boundaries of the model is more than 5 times the width of the excavations to eliminate the boundary effects of the model. The 2D model has 6-node triangular elements, which to study the stresses and displacements in more detail, the elements around the stope are considered to be smaller than other parts of the model. Horizontal displacement along the vertical boundaries and vertical and horizontal displacement along the lower horizontal boundary is prevented. The upper boundary is also modeled according to the overburden depth and its horizontal and vertical displacements are considered to be completely free. Vertical and horizontal stresses according to the depth of the stope, the density of the rock mass, and the ratio of horizontal to vertical stress according to ϑ = 0.3 from Equation 12, K = 0.42 is considered, which is applied to the upper and lateral boundaries of the model. The behavior of the surrounding rock mass is assumed to be elastoplastic and the Mohr-Coulomb behavior model is used for the simulation.
The lowest floor level of the open-pit mine is currently 2830 meters and the lowest floor level of the stope is 2719 meters. Therefore, the vertical in situ stress is 2.6 MPa and the horizontal in situ stress is 1.1 MPa, which is considered by linear distribution with a slope equal to the specific gravity of the mass at height. In numerical models, the ratio of horizontal to vertical stress is 0.42. One-meter steps have been considered for excavation in the rock mass. i.e., the stress release related to the non-supported space is considered to be one meter in each progression step. The cross-section of all stopes and preparatory spaces is considered per reality.

Stability analysis based on the concept of critical strain
Sakurai in 1983 suggested that the stability of tunnels could be assessed based on continuous strain in the surrounding rock mass. This strain is defined as the ratio of the tunnel convergence to its diameter. If the uniaxial compressive strength of the rock mass and its elastic modulus are assumed to be σc and E, respectively. Based on the stress-strain elastic relations, the critical strain (εc) can be calculated by the following equation: In these relations, εcr is the critical strain per percentage (displacement ratio to tunnel radius) and E is the elastic modulus per kg/cm 2 . In risk alert level 1, the tunnel has an instability problem. Strain resulting from risk alert level 2 has been proposed as the basis for tunnel support design and risk alert level 3 indicates short-term stability. Via allowable strain and using Equation 17, the allowable displacement is determined [26], [27]:

Modeling and analysis of batching room stability
According to the available geological information, the rock surrounding the batching is limestone with high geomechanical characteristics, the characteristics of which are given in Table 1, but due to the layering and joints that exist in this area, some local instabilities of batching have been observed. To investigate this situation, two sets of joints are considered in the worst possible case to predict the worst possible conditions (Table 3). Due to the effects of open-pit mine related explosions and the low vertical and horizontal distance between the final pit of the open-pit mine and batching room, it is necessary to design a proper support system so that it can be stable against the waves caused by explosions. For correct modeling of blasting effects, a seismic load coefficient of 0.2 has been considered. Figure 8 shows the model geometry, zoning, rock mass classification, and batching joint. It should be noted that the direction of batching drilling in modeling is in stages and each stage a part of batching is excavated and there has been an effort to make modeling to be very consistent with reality.

Figure 8. Model geometry, zoning, rock mass classification, and batching joint at the end of tunnel No. 2
To investigate the displacements and determine the appropriate support system, first, the batching was excavated and modeled without using the support system, and the displacements and the surrounding plastic zone were measured. The results of these models are given in Figures 9-11, respectively. The contours and displacement vectors around the stope are shown in Figure 9, as well as the maximum plastic area created around the excavated stope in Figure 10. As it is known, the maximum displacement around the batching room is in the sidewall and about 9.5 cm. To investigate the effect of the support system on the batching stability, two types of support systems were considered. Bolt and shotcrete support system is the first type of selective support system in which bolts (according to the specifications of Table 4) with a length of 4 meters and a diameter of 25 mm with shotcrete and mesh with a thickness of 15 cm were selected, and the second support system is shotcrete and a steel frame incorporating IPE140 beam with a one-meter distancing.  The results of modeling these two support systems are shown in Figures 12-17. The amount of displacement, the plastic area around the tunnel after the installation of the shotcrete and mesh support system are shown in Figures 12 and  13, respectively. It is observed that after the installation of the support system, the amount of displacements reaches about half of the initial value and the maximum displacements are about 2.5 cm, which takes place in the right wall and stope toe.
To evaluate the efficiency and safety of the mentioned support system, two diagrams of the axial force-bending moment and axial force-shear force which are known as load capacity diagrams are utilized.

Figure 14. Diagrams of: (a) axial force-shear force of the steel frame; (b) axial force-bending moment of shotcrete batching
These two diagrams are widely used in the design of concrete, engineering structures as well as tunnel support systems. Figure 14 shows the diagrams of the axial forcebending moment and axial force-shear force according to the parameters mentioned for the shotcrete and mesh support system and the loads applied to them in confidence levels of 1 and 1.5. The number of displacements, the plastic area around the tunnel after installation of the steel frame, and the shotcrete support system are shown in Figures 15 and 16, respectively. It can be inferred that after the installation of the support system, the amount of displacements has decreased and the maximum displacements are about 6.4 cm, which occurs in the right wall and stope toe. Figure 17 shows diagrams of the axial force-bending moment and axial force-shear force for the steel frame and shotcrete support system and the loads applied to them in confidence levels of 1 and 1.5. In these diagrams, the value of safety factors represents the range that, if the points drawn for the internal forces are in this range, indicates the stability of the support system. According to the diagrams shown, it can be seen that all points related to both composite support systems are within the range of 1.5 safety factors, which indicates the suitability of the support system. As can be seen from the results, there is not much difference in the results between the shotcrete and bolt and the steel frame support systems, and both support systems control the displacement range and the plastic zone around the batching and place it in the safety range. Therefore, the best support system for batching should be selected according to the priority and considering the technical and economic conditions. Thus, according to the existing conditions, the shotcrete and bolt were considered as the preferred support system.

Stability analysis of batching room
In this paper, to investigate the stability of the underground space, the displacement values obtained from the numerical modeling are compared with the critical displacement values obtained from Sakurai relations and criteria. According to the modeling performed and the results obtained from the Sakurai criterion presented in Table 5, it was found that batching drilling without a support system causes its instability and therefore, the utilization of a support system is warranted. In this research, two different support systems were used to stabilize the excavation space. A comparison between the two selected support systems based on the Sakurai criterion and a comparison of the displacements that occurred around the stopes with the number of critical strains and allowable strains created in the stope shows that there is not much difference between these two support systems and in both cases, the batching stability can be ensured. Thus, according to the stope conditions and available facilities, rock bolt and shotcrete have been selected as the optimal support system.

Conclusions
In this research, finite element modeling has been utilized to analyze the stability and design of the support system of Angouran mine batching room. The results of the analysis show that the excavation of the Batching room without a support system causes its instability and failure, though it is imperative to design a suitable support system according to the geomechanical and stress conditions of the area. The effect of the support system on the stability of the batching room was evaluated by selecting two types of composite support systems of bolt and shotcrete and steel frame with a one-meter distancing and shotcrete based on Sakurai criterion. A comparison between the two selected support systems showed that there was not much difference between the two in terms of performance. Therefore, according to the facilities available in the plant, rock bolt and shotcrete were selected as the final support system. These results are suitable for the analysis of stability and batching design by considering the assumptions about the seismic load coefficient of the blast wave. Due to the proximity of the batching room to the open-pit mine and the direct impact of open-pit minerelated blast waves on its stability due to its high importance and its applicability throughout the life of the mine, to study as closely as possible in this section, it is recommended to install accelerometers, record and study the waves caused by open-pit mine explosions. With comprehensive dynamic modeling and the application of surveying results, more accurate and transparent results of the impact of open-pit mining activities can be achieved.