Nội dung

Journal of Science and Technology in Civil Engineering, NUCE 2020. 14 (3): 67–74 NUMERICAL INVESTIGATION ON THE TUNNELING AND MINING INDUCED GEO-HAZARDS: CASE STUDY IN QUANG NINH, VIETNAM Nguyen Cong Gianga,∗, Nguyen Van Manhb , Nguyen Quang Phichc a Faculty of Civil Engineering, Hanoi University of Architecture, Km 10, Nguyen Trai road, Thanh Xuan district, Hanoi, Vietnam b Faculty of Civil Engineering, Hanoi University of Mining and Geology, No. 18 Vien street, Bac Tu Liem district, Hanoi, Vietnam c Faculty of Civil Construction, Van Lang University, 401C, Campus 45 Nguyen Khac Nhu street, District 1, Ho Chi Minh city, Vietnam Article history: Received 19/05/2020, Revised 29/06/2020, Accepted 29/06/2020 Abstract In the field of rock mechanics, underground construction and mining, there have been many proposed methods for studying geo-hazards and also many research results that have been published in the world. In Vietnam, the numerical method is mainly used for analysis and design but not going deeply to predict the possible causes that lead to geo-hazards due to complex geological conditions. On the other hand, underground constructions and exploitation projects are often designed based on standards, regulations and experiences. The physical mechanism as well as the possibility of geo-hazards occurred when constructing underground structures and mining can take on various forms, depending on geological conditions and construction technology. Therefore, using numerical methods to simulate and analyze the possible geo-hazards is essential. This article presents a number of specific analysis cases, taking into account geological conditions and boundary conditions, and from that, raising a number of issues to note when using numerical methods. Keywords: underground mining; numerical method; geo-hazards; rock mechanics; FLAC2D. https://doi.org/10.31814/stce.nuce2020-14(3)-06 c 2020 National University of Civil Engineering 1. Introduction Geo-hazards are types of disaster-related to geological processes induced by natural or human activities. In recent years, various geo-hazard or geo-risks have occurred in tunneling and mining in Vietnam. In the mining field, the geo-hazard is due to exploitation of natural resources from mine including subsidence, slope stability, landslides and other related damage have been reported by numerous authors [1–4]. The prediction and management of geo-hazard are of great importance in the mining industry [5–7]. Understanding the behaviour of rock masses has always been difficult for mining and underground engineers because of the presence of discontinuities, anisotropic and heterogeneity. Empirical, analytical and numerical methods have been widely used for modeling the behavior of rock mass [8–10]. ∗ Corresponding author. E-mail address: gianglientca@gmail.com (Giang, N. C.) 67 Giang, N. C., et al. / Journal of Science and Technology in Civil Engineering In recent years, numerical methods have been used in design of underground openings in the world, Giang, N. C., of Science andand Technology in Civil Engineering Giang, N. et C.,al./Journal et al./Journal of Science Technology in Civil Engineering however, in Vietnam, this issue has received little attention. Numerical modeling application in mining engineering aims to provide a better understanding of the mining and rock mechanics engineers for solving problems related to the design of support systems [7, 11]. The numerical methods are conthis article, rock mass layers andfor underground aditadit in Quang Ninh province of of Incostly this article, the rock mass layers and underground in Quang Ninh province venient,Inless andtheless time-consuming the analysismine ofmine stress redistribution and their effects Vietnam were selected to investigate the influence of underground mine adit location and rock layers Vietnam wereofselected to investigate the influence of underground mine and rock layers on the behavior rock mass and designing of support system within the adit rocklocation mass environment. positions to stress states, yielded zonezone and and displacement of the rockrock mass surrounding aditadit by using positions to stress states, yielded displacement of the mass surrounding using In this article, the rock mass layers and underground mine adit in Quang Ninh province ofby Vietnam FLAC2D [12].[12]. TheThe results suggest thatthat the the subsidence of the surface could be triggered duedue to to FLAC2D results suggest subsidence of the surface could be triggered were selectedcollapse. to collapse. investigate the influence of underground mine adit location and rock layers positions underground underground to stress states, yielded zone and displacement of the rock mass surrounding adit by using FLAC2D 2. Mechanical and Simulation Diagram 2.The Mechanical Parameter and Simulation [12]. resultsParameter suggest that the subsidence of Diagram the surface could be triggered due to underground collapse.TheThe rockrock massmass in the coalcoal mining in Quang Ninh province, Vietnam consists normally of 4of 4 in the mining in Quang Ninh province, Vietnam consists normally layers: sandstone, siltstone, argillite and and coalcoal lying inclined withwith mechanical parameters as in layers: sandstone, siltstone, argillite lying inclined mechanical parameters as the in the following Table 1. following Table 1. 2. Mechanical parameter and simulation diagram Table 1. Mechanical parameters of rock mass Table 1. Mechanical parameters of rock mass The rock mass in the coal mining in Quang Ninh province, Vietnam consists normally of 4 layers: Layer Density Angle Modulus Shear Modulus Layersiltstone, Density Cohesion cFriction Friction Angle Bulk Modulus Shear Modulus sandstone, clayCohesion and coalc lying inclined with Bulk mechanical parameters as in G the Gfollowing K Table 1. The constitutive model of Mohr–Coulomb is used for modeling the behaviour of the rock K r r j j (MPa) (GPa) (MPa) (GPa) mass. 3 3 (g/cm ) ) (g/cm Sandstone Sandstone 2.612.61 Coal Coal Layer Slitstone Slitstone Sandstone Clay Clay Coal Slitstone Clay Density 1.301.30ρ (g/cm3 ) 2.502.50 2.61 2.60 2.60 1.30 2.50 2.60 (GPa) (GPa) (degrees) (degrees) Table 1. Mechanical parameters of rock mass 1.001.00 40 40 11.60 11.60 Cohesion 0.010.01 c (MPa) 1.001.00 8.708.70 Friction Bulk Shear modulus G 35 35 angle ϕ 2.60 2.60modulus K1.301.30 (degrees) (GPa) (GPa) 1.00 0.100.01 0.10 1.00 0.10 25 25 40 30 3035 25 30 10.00 10.00 11.60 9.609.60 2.60 10.00 9.60 7.007.00 2.702.70 8.70 1.30 7.00 2.70 TwoTwo cases are investigated withwith different orders of rock layers: cases are investigated different orders of rock layers: - Case 1: from upper to lower layers are clay-coal-clay-siltstone layers - Case 1: from upper to lower layers are clay-coal-clay-siltstone layers - Case 2: cases from upper to lower layers are sandstone-clay-coal-siltstone-sandstone -Two Case 2: from to lower layers are sandstone-clay-coal-siltstone-sandstone areupper investigated with different orders of rock layers: TheThe tunnel has a semi-circular and straight-walled shape or D-shape with a width andand height of 4m tunnel has a semi-circular and straight-walled shape or D-shape with a width height of 4m - Case 1: from upper to lower layers are clay-coal-clay-siltstone layers. each, excavation in coal layer. TheThe height fromfrom top top of the aditadit to the surface is nearly 30m. TheThe each, excavation in coal layer. height of the to the surface is nearly 30m. analysis model is shown in Fig. 1. 1. analysis model is shown in Fig. (a) Case 1 (b) Case 2 (a) Case (a) Case 1 1 (b) Case (b) Case 2 2 Figure 1.1.The model cases study Figure The model ofthe thetwo two cases study Figure 1. The model of of the two cases study 3. Simulation results discussions 3. Simulation results andand discussions 68 Redistribution of stress induced to excavation opening a complex subject in actual Redistribution of stress induced due due to excavation opening is aiscomplex subject in actual mining conditions because of the influence of rock mass layers. results of the numerical mining conditions because of the influence of rock mass layers. TheThe results of the numerical simulation show all information on the of mechanical changes occurring in the mass simulation can can show all information on the lawslaws of mechanical changes occurring in the rockrock mass Giang, N. C., et al. / Journal of Science and Technology in Civil Engineering - Case 2: from upper to lower layers are sandstone-clay-coal-siltstone-sandstone. The tunnel has a semi-circular and straight-walled shape or D-shape with a width and height of 4 m each, excavation in coal layer. The height from top of the adit to the surface is nearly 15 m. The model was built with size of 30 × 30 m. The left and right boundary of model are fixed at horizontal direction; the bottom boundary is fixed in both vertical and horizontal direction and the top boundary of model is free. The analysis model is shown in Fig. 1. The initial boundary condition of the model is in-situ rock mass stress state. 3. Simulation results and discussions Redistribution of stress induced due to excavation opening is a complex subject in actual mining conditions because of the influence of rock mass layers. The results of the numerical simulation can Giang, N. C., et al./Journal of Science Technology in the Civil Engineering Giang, N.on C.,the et al./Journal of Science andand Technology in in Civil Engineering show all information laws of mechanical changes occurring rock mass surrounding the adit, including the stress redistribution, displacement, deformation and the failure zones. Based on that information the designers can analyze and choose the possibilities of support systems for reinforcing rock mass keep the stability oftothekeep opening. By introducing the simulation results, the systems for reinforcing rock mass tounderground keep stability of the underground opening. introducing systems forto reinforcing rock mass thethe stability of the underground opening. ByBy introducing thethe advantage of numerical simulation in general as simulation well as in be demonstrated. simulation results, advantage of numerical in general asanalysis well ascould in geohazards analysis simulation results, thethe advantage of numerical simulation in geo-hazard general as well as in geohazards analysis could demonstrated. could be be demonstrated. (a) Case 1 (b) Case 2 Case (b)(b) Case 2 2 Case (a) (a) Case 1 1 Figure 2. Major principal stress distribution (Pa) Figure 2. Major principal stress distribution Figure Major principal stress distribution 23and 3 show the distribution of the major minor principal stresses in the the rock mass 2 and 3 show distribution the major and minor principal stresses in in the rock mass Figs.Figs. 2Figs. and show thethedistribution ofofthe major andand minor principal stresses rock mass surrounding tunnel. is easy to see in Fig. 2 that general redistribution principles surrounding thethe tunnel. It isIt easy to see in Fig. 2 that thethe general redistribution principles areare thatthat thethe surrounding the tunnel. It is easy to see in Fig. 2 that the general redistribution principles are that major principal stresses peak value occur hard rock layers, however order major principal stresses of of peak value occur in in thethe hard rock layers, however duedue to to thethe order of of thedifferent majorrock principal stresses ofbe peak value in the hard layers, however due to the order of rock layers, it can seen that the average value of principal stress forming different different layers, it can be seen that theoccur average value of rock thethe principal stress forming different different rock layers, it can be seen that the average value of the principal stress forming different shapes. case 2, as siltstone under layer, which is mechanically stronger than shapes. In In case 2, as thethe siltstone lieslies under thethe coalcoal layer, which is mechanically stronger than thethe shapes. In the case 2, as the thethe siltstone lies under coal layer, which isMPa mechanically stronger the clay layer, the area of principal stress with an average value of 0.25 MPa spreads deeper in than the clay layer, area of principal stress with anthe average value of 0.25 spreads deeper in the coalcoal than in case than in case 1. 1. clay layer, the area of the principal stress with an average value of 0.25 MPa spreads deeper in the coal than case 1. thethe By comparing minimum principal stress distribution principles in Fig. 3, the findings show Byin comparing minimum principal stress distribution principles in Fig. 3, the findings show By comparing the minimum principal distribution principles 3,with the findings show quite similarities in the distribution areas instress both cases. However, in case 2, Fig. the with the smallest quite similarities in the distribution areas in both cases. However, in case 2, in the areaarea the smallest minor stress components with values ranging from toHowever, 0.05 MPa iscase more than insmallest case minor stress components with values ranging from 0 to0 0.05 MPa isinmore widespread in 1 1 quite similarities in the distribution areas in both cases. 2, widespread the area than with thecase with the stress fluctuating the range 0.15 It smaller is is smaller lower part of the with the stress fluctuating in in thevalues range of of 0.15 to to 0.200.2 MPa. is in in thethe lower part of in thecase minor stress components with ranging from to MPa. 0.05It MPa more widespread than modeling. modeling. 1 with the stress fluctuating in the range of 0.15 to 0.2 MPa. It is smaller in the lower part of the Therefore, results showed effect layers very clear stress Therefore, thethe results showed thatthat thethe effect of of thethe layers is is very clear to to thethe stress modeling. redistribution, which very different from results obtained analytical methods with redistribution, which is is very different from thethe results obtained by by analytical methods with thethe "averagation" or “homogenization” of the model rock mass [13]-[15]. "averagation" or “homogenization” of the model on on thethe rock mass [13]-[15]. 69 clay layer, area of the principal stress an average value of 0.25 spreads deeper in coal the coal clay layer, thethe area of the principal stress withwith an average value of 0.25 MPaMPa spreads deeper in the than in case 1. than in case 1. comparing minimum principal stress distribution principles in Fig. 3, findings the findings show ByBy comparing thethe minimum principal stress distribution principles in Fig. 3, the show quite similarities in the distribution areas in both cases. However, in case 2, the the smallest quite similarities in the distribution areas in both cases. However, in case 2, the areaarea withwith the smallest minor stress components with values ranging from to 0.05 is more widespread in case minor stress components with values ranging from 0 to0 0.05 MPaMPa is more widespread thanthan in case 1 1 Giang, N. C., et al. / Journal of Science and Technology in Civil Engineering with stress fluctuating in the range of 0.15 to 0.2 MPa. is smaller in the lower of the with thethe stress fluctuating in the range of 0.15 to 0.2 MPa. It isIt smaller in the lower partpart of the modeling. modeling. Therefore, the results showed that the effect of the layers is very clear to the stress redistribution, which isTherefore, very different from the results by analytical methods the “averagation” or Therefore, results showed that effect layers is with very clear to the stress thethe results showed thatobtained the the effect of of the the layers is very clear to the stress redistribution, which very different from results obtained analytical methods redistribution, which very different the the results obtained by by analytical methods withwith the the “homogenization” of istheis model on the from rock mass [13–15]. "averagation" or “homogenization” of the model on the mass [13]-[15]. "averagation" or “homogenization” of the model on the rockrock mass [13]-[15]. (a) (a) Case 11 1 (a) Case Case (b) Case Case (b) Case 2 22 (b) 3. Minimum principalprincipal stress distribution (Pa) Figure 3.Figure The principle of minimum stress distribution Figure 3. The principle of minimum principal stress distribution Similarly, the results obtained with the principle of movement show that due to the influence of the Giang, N. et al./Journal ofmass Science Technology in Civil Engineering rock mass layers, theN.movement in the rock around the underground opening is not symmetrical Giang, C.,C., et al./Journal of Science andand Technology in Civil Engineering but dependent on the specific geological structures. Fig. 4 shows the displacement on the boundary of the opening, reflected across the boundary of3the3 opening after the displacement. (a) Case 1 Case (a) (a) Case 1 1 (b) Case 2 Case (b)(b) Case 2 2 Figure 4. Displacement of the opening boundary after excavation Figure 4. Displacement of the opening boundary after excavation Figure 4. Displacement of the opening boundary after excavation Figs. 5 and 6 show the formation of failure zone (area with symbol) in the rock mass around the Similarly, results obtained with principle show to the opening. The failure zone in obtained both cases develops mainly inof themovement coal layer andthat tothat the surface ofinfluence the hard Similarly, thethe results with thethe principle of movement show duedue to the influence of the rock mass layers, movement in rock mass around underground opening is of themass. rock mass layers, thethe movement thethe rock mass around underground opening notnot rock However, comparing the twoincases with different orderthe ofthe rock mass layers showsis that in symmetrical but dependent on the specific geological structures. Fig. 4 shows the displacement on the symmetrical but dependent on the specific geological structures. Fig. 4 shows the displacement on the case 2, the failure zone is wider. The results clearly show the influence of the stratification structure boundary the of the opening, reflected across boundary of the opening after displacement. boundary opening, reflected across thethe boundary of the opening after thethe displacement. as well asofthe order of the rock mass layers on the formation of geohazards. The failure state occur Figs. 5stress and 6 isshow formation astrength failure zone (area with symbol) inThe the rock mass around when the shear more than the of shear of(area rock mass element. symbol of “*” in Figs. 5 and 6 show thethe formation aoffailure zone with symbol) in the rock mass around the opening. The failure zone in both cases develops mainly in the coal layer and develop the The failure bothelement cases develops mainly in the and coalthe layer and develop to to thethe Fig.opening. 5 is indicated that thezone rockinmass was failed by shearing symbol of “x” indicating surface of the hard rock mass. However, comparing the two cases with different order of rock mass surface the hard rockwas mass. However, comparing the rockofmass element failed in elasticity state. the two cases with different order of rock mass layers shows that in case 2, the failure zone is wider. The results clearly show the influence of layers shows that in case 2, the failure zone is wider. The results clearly show the influence of thethe stratification structure as well as the order of the rock mass layers formation of geohazards. 70mass stratification structure as well as the order of the rock layers on on thethe formation of geohazards. Figs. 5 and 6 show the formation of a failure zone (area with symbol) in the rock mass around Figs. 5 and 6 show the formation of a failure zone (area with symbol) in the rock mass around the opening. The failure zone in both cases develops mainly in the coal layer and develop to the the opening. The failure zone in both cases develops mainly in the coal layer and develop to the surface of the hard rock mass. However, comparing the two cases with different order of rock mass surface of the hard rock mass. However, comparing the two cases with different order of rock mass layers shows that in case 2, the failure zone is wider. The results clearly show the influence of the layers shows that in case 2, the failure zone is wider. The results clearly show the influence of the stratification structure as N. well the/ Journal order of the rockandmass layers in on the Engineering formation of geohazards. Giang, C., as et al. of Science Technology stratification structure as well as the order of the rock mass layers onCivil the formation of geohazards. 5. Failure zone around Casex:1elastic failure) Figure 5. FailureFigure zone around the opening, case 1the (*:opening: shear failure; Figure 5. Failure zone around the opening: Case 1 Figure 6. FailureFigure zone around the opening, case 2the (*:opening: shear failure; 6. Failure zone around Casex:2elastic failure) Figure 6. Failure zone around the opening: Case 2 Several comments can be drawn from the numerical results: 4 4 principal stress redistribution, displacement of - When the rock mass has a layered structure, tunnel boundary and the formation of failure zones are complex. Rock formation does not behave as homogeneous material, it requires therefore advanced numerical model to solve the problem; - It is clear that the mechanical behavior of rock mass depends not only on the location of the opening, but also on the order and distribution of the rock mass layers, which are clearly shown in the numerical results; - In the second model, the displacement and deformation processes achieve relatively larger values, although in case 2 there are both sandstone layers in pillars and cliffs; - The change in the position of the layers clearly affects the processes of stress redistribution and movement in the rock masses; - The development of the failure zone in the latter case is stronger; - In both cases, as the distance from the top of the structure to the surface is not wide the failure zone is developed to the surface of hard rock mass layers. In this case, it may cause landslide or land subsidence, with varying intensity. 4. Influence of tunnel shape on geomechanical process To study the effect of the cross-section shape of underground opening to the redistribution of stress states and failure zone, simulations were performed with the case of the circular shape which 71 subsidence, with varying intensity. 4. Influence of tunnel shape on geomechanical process To study the effect of the cross-section shape of underground opening to the redistribution of stress states and failureGiang, zone,N. simulations were performed with in the case of the circular shape which C., et al. / Journal of Science and Technology Civil Engineering has a radius of 2m, with similar mechanical parameters and order of distribution of rock mass layers as has a radius of 2 m, with similar mechanical parameters and order of distribution of rock mass layers in Section 2. The obtained results show that, when the opening is a circular shape, the rules of stress as in Section 2. The obtained results show that, when the opening is a circular shape, the rules of stress redistribution and displacement also show dependence on the layering of the rock mass. However, the redistribution and displacement also show dependence on the layering of the rock mass. However, the failure zone does not develop to the surface. It also means that it is not likely to lead to landslides or failure zone in this case does not develop to the surface. It also means that it is not likely to lead to land subsidence under the investigated conditions. landslides or land subsidence under the investigated conditions. 7–9 show the the rulerule of the and minor stress redistribution as well as the Figs. Figs. 7, 8 and 9 show ofmajor the major and principal minor principal stress redistribution as failure well as the of of rock mass layers are sandstone, claystone, coal andcoal claystone. failure zone, zone,for forthe thecase case rock mass layers are sandstone, claystone, and claystone. Giang, N. et C.,al./Journal et al./Journal of Science Technology in Civil Engineering Giang, N. C., of Science and and Technology in Civil Engineering Figure 7. Modeling of a circular opening Figure 7. Modeling of a circular opening The numerical modeling results reveal that when the opening is a circular cross-section shape, the rules of principal stress redistribution, displacement and deformation strongly depend on the shape of the opening compared to the two cases analyzed above (Figs. 3 and 8). The failure zone is formed within the coal layer, although there are also some local failure points on the hard rock mass, which are not symmetrical, due to the inclined rock mass layers. Especially, when paying attention to ground subsidence and landslides, it shows that when selecting a circular cross-section of opening the land subsidence decreases, and it seems difficult that the failure reaches the surface (Fig. 9). 5 (a) (b) Figure 8. The principle of maximum stress redistribution (left)(left) and and minimum (right) withwith the circular Figure 8. The principle of maximum stress redistribution minimum (right) the circular Figure 8. The principle of (a) maximum and (b) minimum stress (Pa) redistribution with the circular cross-section of opening cross-section of opening cross-section of opening The numerical modeling results reveal that when the opening is a circular cross-section shape, the rules of principal stress redistribution, displacement and deformation strongly depend on the shape of the opening compared to the two cases analyzed above (Figs. 3 and 8). The failure zone is formed 72 Giang, N. C., et al. / Journal of Science and Technology in Civil Engineering within the coal layer, although there are also some local failure points on the hard rock mass, which symmetrical, due to the inclined rock mass layers. Especially, when paying attention to ground Figureare8.not The principle of maximum stress redistribution (left) and minimum (right) with the circular subsidence and landslides, it shows that when selecting a circular cross-section of opening the land cross-section ofthe opening subsidence decreases, and it is difficult to collapse to surface (Fig. 9). Figure Figure 9. The failure zone aroundzone the circular (*: shear failure; x:tunnel elastic failure) 9. The failure (slash)tunnel around the circular 5. Conclusion 5. Conclusions The numerical modeling results show that the numerical methods can solve the complex numerical modeling resultsinshow the numerical methods can solve complexbehavior prob- of problems The of the tunneling and mining rockthat masses by taking into account thethe complex lems of thesuch tunneling and mining in rock masses byheterogeneities. taking into account therules complex behavior of rock rock formation as discontinuities, anisotropic, The for stress redistribution, formation as discontinuities, and as heterogeneities. rules for stress redistribution deformation, thesuch development of the anisotropic failure zone, well as theirThe magnitude, depend clearly on the and characteristics, deformation, the development of theof failure zone as welllayers, as theircross-section magnitude, depend clearly on structural the arrangement the rock mass shapes of opening. the structural characteristics, arrangement of the rock mass layers and cross-section shapes of openObviously, in order to obtain accurately, the mechanical behavior in rock mass with complex ing. Obviously, to obtain accurately, mechanical behavior in rock complex geological structures,initorder is necessary to analysethe specificallytion case by case.mass On with the one hand, by geological structures, it is necessary to analyse specifically case by case. On the one hand, by using using numerical method, it is possible to analyze the possibility and type of development of specific numerical method, it islead possible to analyzeand theaccidents", possibility and type ofthat development of specific ge- the geological conditions that to "incidents meaning it is possible to identify conditions that lead to “incidents and accidents”, meaning that On it is the possible identify type ofological "geological disaster" which can be caused by human factors. othertohand, bythe clearly type of “geological disaster” which can be caused by human factors. On the other hand, by clearly understanding the rock masses behavior after excavated opening, it can be helpful for the designer to the rock masses after method, excavatedsupport opening,system it can befor helpful for the designer to the select understanding suitable cross-section shape,behavior excavation the opening to prevent select suitable cross-section shape, excavation method and support system for the opening to prevent possibility of incidents and accidents, i.e. limiting the geological disaster. the possibility of incidents and accidents, i.e. limiting the geological disaster. References [1] Marschalko, M., Yilmaz, I., Bednárik, M., Kubečka, K. (2011). Variations in the building site categories in the underground mining region of Doubrava (Czech Republic) for land use planning. Engineering Geology, 122(3-4):169–178. 6 73 Giang, N. C., et al. / Journal of Science and Technology in Civil Engineering [2] Marschalko, M., Yilmaz, I., Křístková, V., Fuka, M., Bednarik, M., Kubečka, K. (2012). Determination of actual limit angles to the surface and their comparison with the empirical values in the Upper Silesian Basin (Czech Republic). Engineering Geology, 124:130–138. [3] Yilmaz, I., Marschalko, M. (2012). A leaning historical monument formed by underground mining effect: an example from Czech Republic. Engineering Geology, 133:43–48. [4] Swift, G. (2014). Relationship between joint movement and mining subsidence. Bulletin of Engineering Geology and the Environment, 73(1):163–176. [5] Marschalko, M., Yilmaz, I., Bednárik, M., Kubečka, K. (2012). Influence of underground mining activities on the slope deformation genesis: Doubrava Vrchovec, Doubrava Ujala and Staric case studies from Czech Republic. Engineering Geology, 147:37–51. [6] Yang, Y.-Y., Xu, Y.-S., Shen, S.-L., Yuan, Y., Yin, Z.-Y. (2015). Mining-induced geo-hazards with environmental protection measures in Yunnan, China: an overview. Bulletin of Engineering Geology and the Environment, 74(1):141–150. [7] Singh, G. S. P., Singh, U. K., Murthy, V. M. S. R. (2010). Applications of numerical modelling for strata control in mines. Geotechnical and Geological Engineering, 28(4):513–524. [8] Jing, L. (2003). A review of techniques, advances and outstanding issues in numerical modelling for rock mechanics and rock engineering. International Journal of Rock Mechanics and Mining Sciences, 40(3): 283–353. [9] Harrison, J. A., Hudson, J. P. (2000). Engineering Rock Mechanics: Part 2: Illustrative Worked Examples. Pergamon, Oxford, UK. [10] Kien, N. T. (2020). Multi-scale modeling of geomechanics problems using coupled finite-discreteelement method. Journal of Science and Technology in Civil Engineering (STCE)-NUCE, 14(1V):93–103. (in Vietnamese). [11] Hussain, S., Ur Rehman, Z., Mohammad, N., Tahir, M., Shahzada, K., Wali Khan, S., Salman, M., Khan, M., Gul, A. (2018). Numerical modeling for engineering analysis and designing of optimum support systems for headrace tunnel. Advances in Civil Engineering, 2018. [12] Itasca Consulting Group, Inc. FLAC2D User’ guide. [13] Mindlin, R. D. (1940). Stress distribution around a tunnel. Transactions American Society of Civil Engineers, 105. [14] Mindlin, R. D. (1948). Stress distribution around a hole near the edge of a plate under tension. Proceedings of the Society of Experimental Stress Analysis, 5(2):56–68. [15] Phich, N. Q. (2007). Rock mechanics. Construction Publishing, Hanoi. (in Vietnamese). 74