Numerical Simulation of Surface Movement and Deformation Caused by Underground Mining with Complex Stratigraphic Boundary

/e interface program of finite element software based on surface spline interpolation is developed by using MATLAB software. /e controllable 3D surface modeling based on CAD contour is realized. Taking a mine as an example, the method of establishing the 3D numerical calculation model including complex stratum boundary is studied. /e influence of underground mining on surface movement and deformation under complex stratum conditions by using the FLAC3D software further was discussed./e research results show that the developed interface program of finite element software can well realize the numerical modeling of complex formation and provide help for subsequent numerical simulation. /e calculated subsidence value is in good agreement with the measured value./e values of surface tilt, surface curvature, and surface horizontal deformation are less than the relevant regulations. /e mining method of the filling method has no obvious effect on the surface structures. /e research results can provide reference for similar numerical simulation.


Introduction
With the rapid development of economy, people's demand for mineral products is increasing, which makes the mining industry develop vigorously. e surface subsidence disaster brings hidden danger to people's safety [1][2][3]. e underground mining can cause deformation and damage to surrounding rock mass, surface structures, railways, and water bodies. erefore, the surface movement and deformation caused by underground mining and the damage of surface buildings have become one of the global manmade environmental disasters and a hot issue in the mining field [4][5][6].
With the continuous improvement of mining related theory and technology, many achievements have been made in the study of surface subsidence caused by underground mining. At present, the prediction methods of surface movement and deformation can be divided into three categories. ey are the theoretical analysis method, similarity model, and numerical simulation method. In terms of the theoretical analysis method and the similar model method of surface movement and deformation caused by underground mining, Sasaoka et al. [7], Villegas et al. [8], Fazio et al. [9], Wu et al. [10], and He and Kang [11] have carried out a very meaningful discussion and practice on the theory of surface movement and deformation caused by underground mining. e theoretical analysis method can analyze and restore the surface movement and deformation curve by introducing related variables. However, the analysis of complex strata and terrain related problems is slightly insufficient. Nikolaos et al. [12], Guarino et al. [13], Scotto di Santolo [14], Yin et al. [15], and Li et al. [16] analyzed the surface movement caused by underground mining by simulation experiment method. However, the model experiment is established based on the similarity theory. erefore, the material properties, stress factors, and boundary conditions are difficult to fully meet the similar requirements. Moreover, the model test process is complex.
At present, with the development of computer hardware and the progress of related numerical model software, the numerical simulation method of surface movement and deformation of underground mining has made rapid progress [17,18]. e influence of dynamic load, groundwater, and other comprehensive factors on the surface movement and deformation can be considered in the numerical simulation [19]. erefore, the numerical simulation has become a common analysis method [20][21][22][23]. Hadi et al. [24] used the FLAC3D software to simulate the ground surface subsidence due to sublevel caving method. en, the longitudinal and transverse profiles of subsidence and the influence zone were studied. Zuo et al. [25] used a new program MDDA (mining discontinuous deformation analysis) which has been developed to simulate the continuous excavation process in mining engineering based on the existing discontinuous deformation analysis (DDA). Both the real-time stress distribution and evolution of rock strata movement during the mining process could be effectively obtained. Nick et al. [26] used a distinct element numerical model PFC2D to establish a large-scale avalanche of earth material. A general stretching and thinning of the avalanche are observed.
However, it is difficult to convert complex strata and complex terrain into the 3D numerical model. e establishment process of 3D stratigraphic boundary is complex. At present, the simplified stratigraphic boundary modeling is generally adopted by making a plane with consistent occurrence along a single direction, although some software can directly use the complex 3D stratum boundary to construct the model. However, it is very difficult to control the model element size at the 3D stratigraphic boundary, due to the complexity of CAD isoline spacing and nodes. e high-density grid is often generated in the area with large height difference (dense contour lines), which makes the total number of model units increase sharply. e existence of high-density mesh seriously affects the calculation speed of the 3D model.
Based on the self-developed CAD and finite element software interface program, this paper realizes the controllable grid of complex strata. Taking a mine as an example, a mine numerical model with complex geological strata is established by using the FLAC3D numerical simulation method. e relevant horizontal and vertical sections, geological boundary contour map, ore body distribution map, and other data are used in the model. rough the field engineering geological survey and rock mechanics, the mine numerical model is established. e actual rock mass mechanical parameters are obtained from the test. en, the calculation and analysis are carried out. In this way, not only the influence of complex strata is considered but also the influence of joint and other structural planes is reflected in the process of rock mass mechanical parameters. e surface subsidence is in good agreement with the measured values. e calculation results of surface movement and deformation can provide reliable guidance for mining.

Mesh Construction Method of 3D
Complex Surface e MATLAB program is used to realize the interface between CAD and finite element software. e 3D complex surface model can be directly generated from CAD data, which can be used by finite element software. e features of the interface program are as follows: (1) e interface program can automatically extract the coordinates of isoline points and corresponding elevation values and coexist in the text file. (2) e interface program can automatically set the horizontal and vertical dimensions of the 3D surface mesh. e difference of the grid near the actual coordinate point can be found automatically. Finally, the regular and uniform 3D surface is obtained. (3) e interface program can realize the automatic extension of isoline to nonisoline data area. e automatic construction of area surface model without isoline data can be realized.
If the density of isoline points extracted by the program does not meet the needs of horizontal and vertical grid size of 3D surface, the surface spline interpolation method is used for encryption. e coordinates and corresponding elevation values of the points before encryption are where X i and Y i are the X and Y coordinates of elevation point I, respectively, and H i is the elevation value of point i. Formula (1) is a bivariate single-valued list function, which is fitted. e expression of bivariate spline function is as follows: where r 2 � (x − x i ) 2 + (y − y i ) 2 , c 1 , c 2 , . . . , c 3+n is the coefficient to be solved, and ε is the adjustment coefficient, ranging from 0.01 to 1.00 for flat areas and 10 −6 to 10 −5 for singular surfaces. e coefficient to be calculated can be determined by the following formula: 2 and h j is the weighting coefficient of the j node, which can be taken as 0. Equation where C m×1 is a symmetric matrix composed of node coordinate values and weighted coefficients. If A m×m is not a singular matrix, then the equation can be solved and the coefficient matrix is en, equation (2) can determine that if the model area is divided into k quadrilateral grids, the corresponding isoline elevation of each grid node is expressed as e partial derivative of elevation function of isoline at any grid node is expressed as c 3+i is the coefficient to be solved. k is the number of quadrilateral meshes. i is the specific grid node. ε is the adjustment coefficient, ranging from 0.01 to 1.00 for flat areas and 10 −6 to 10 −5 for singular surfaces.
rough formula (7) and grid node coordinates, the corresponding coordinates of the mesh nodes of the fitting surface can be obtained. e three-dimensional coordinates of the interpolated points can be obtained.
Taking the CAD contour, as shown in Figure 1, as an example, the grid length in horizontal and vertical directions (25 m in horizontal direction and 20 m in vertical direction) is set independently in the interface program. e local part of CAD contour line (such as the lower left corner of CAD drawing) cannot meet the grid density. e program carries out automatic interpolation calculation to generate the 3D surface model. e different grayscale in Figure 2(a) represents different elevation. e model is consistent with the isoline in Figure 1.

General Situation of Mining Area.
e surface of the mining area is gentle, slightly inclined to the east, with a gradient of about 7‰. e surface of the mining area is approximately horizontal. e strata in the mining area are composed of quaternary loose rocks, neogene clastic rocks, metamorphic rocks of Daqueshan formation of Middle Proterozoic Zhujiashan group, and ultrabasic rocks. ere are mainly two ore belts K1 and K2 in the mining area. e K1 ore belt occurs in the top alteration of the rock mass. e K2 mineralization belt occurs in the alteration of the rock bottom. e occurrence of K1 and K2 ore belts is consistent with the top and bottom of the rock mass, respectively.
A large number of buildings and structures are distributed on the surface of the mining area.
e several villages and a provincial highway are mainly distributed within the scope of the mining right. If the surface is destructively deformed, it will have serious consequences. is paper will focus on the analysis of the surface movement and deformation of underground mining. e geological boundary of strata in the mining area is complex. e contour line of typical rock stratum boundary is shown in Figure 3. Because the theoretical analysis method has strict requirements on the occurrence of strata, it is very difficult to establish the complex geological boundary in the similar model test.
erefore, it is very practical to use numerical simulation analysis for the mine with complex stratum boundary.

Construction of Numerical Model.
According to the actual situation of the mine, the developed interface program is used to realize the three-dimensional surface modeling of multiple complex strata, which provides the basis for the subsequent surface movement and deformation analysis. According to the geological exploration data, a three-dimensional grid of quaternary stratigraphic boundary and surface topography is established, as shown in Figure 4. According to the contour map of stratigraphic boundary, the boundary 3D grid of upper K1 ore body and lower K2 ore body is established, as shown in Figure 5. e final 3D grid of stratigraphic boundary is shown in Figure 6.
According to the established three-dimensional grid of stratum boundary, the numerical calculation model is completed and the units are divided. According to the Advances in Civil Engineering horizontal projection of ore body, the high-and low-grade ore bodies are divided. e three-dimensional model of the ore body and surrounding rock three-dimensional model is shown in Figure 7. e different colors represent different grades of ore bodies. e two-layer ore bodies are composed of four three-dimensional curved surfaces. e internal units of the calculation model are shown in Figure 8. In Figure 8, group 9 is the quaternary system, group 11 is the Neogene rock mass, group 10 is the roof rock mass of K1 seam, group 7 is the roof rock mass of K1 seam and K2 seam, group 8 is the floor rock mass of K2 seam, and group 3, group 5, and group 6 are different grade ore bodies.   In the numerical analysis of underground engineering, the key to successful calculation is not only the accurate calculation model but also the initial geo-stress field [27,28].
is time, the load is applied according to the self-weight stress of the rock mass. e z-direction constraint is applied at the bottom of the model, and the X-direction and Ydirection constraints are applied around the model. e surface is a free surface.

Numerical Simulation Scheme and Parameter Value.
is numerical simulation mainly analyzes the influence law and scope of mining on the surface. According to the mine design, the simulated mining steps are determined. e filling is carried out at the end of each mining step. e specific numerical simulation scheme is shown in Figure 9.
According to the results of rock physical and mechanical test and rock mass quality evaluation, the mechanical parameters of surrounding rock mass are calculated according to Hoek-Brown criterion. e typical Hoek-Brown calculation curve is shown in Figure 10. e final mechanical parameters of rock mass are shown in Table 1.

Analysis of Numerical Simulation Results.
After the highgrade ore above 280 m in K1 ore body is mined and filled, the displacement of the upper rock mass next to the filling body is the largest, which is 0.15 m. e surface displacement above the mining area is the largest, about 0.05 m. After the mining and filling of high-grade ore below −280 m of K1 ore body and −730 m∼−590 m of K2 ore body, the maximum displacement of rock stratum appears above the filling body of K2 ore body, which is about 0.3 m. e surface displacement is about 0.1 m. After the high-grade ore above −730 m and below −590 m of K2 ore body is mined and filled, the maximum displacement of strata is located above the filling body of K2 ore body, and the increase is not obvious, which is 0.32 m. e surface displacement reaches 0.2 m. e mining of all low-grade ore bodies has a significant impact on strata movement. e maximum displacement inside the strata reaches 0.37 m, which is located in the upper part of K2 ore body and relatively close to the surface. At this time, the maximum displacement of the surface reaches 0.3 m. e cloud chart of strata movement in different mining stages is shown in Figure 11. e results of FLAC3D are imported into surfer software, and the final deformation of the surface is shown in Figure 12.

Analysis of eoretical Calculation
Results. At present, the mine has arranged the surface subsidence observation points along the 1-1 survey line. e mining of the lower ore body of Miaozhuang has been completed, and the mining area is far away from here. e measured results show that the surface subsidence has become stable. e maximum measured value of subsidence near village 4 is 0.07 m, which can be regarded as the final subsidence value of this point.
ere is no obvious damage sign of surface structures. FLAC3D can directly export the surface subsidence isoline map and the surface horizontal movement isoline map. Surface slope contour map and surface curvature contour map are transformed into surface horizontal deformation contour map. It can be realized by the derived surface subsidence and surface horizontal movement data by the MATLAB, according to the following calculation formula. e formula for calculating the inclination of the earth's surface in x and y directions is where W(x, y) is the surface subsidence value. e formula for calculating the curvature of the earth's surface in x and y directions is where i (x, y) is the surface inclination value in the x and y directions. e formula for horizontal deformation in the direction of x and y of the surface is where U (x, y) is the horizontal surface movement in x and y directions. e surface subsidence isoline map, surface horizontal movement contour map, and corresponding survey line profile map obtained by numerical simulation are shown in Figures 13 and 14. e surface slope contour map, surface curvature contour map, surface horizontal deformation contour map, and corresponding survey line profile map obtained by formulas (8)- (10) are shown in  e final surface movement deformation value obtained is shown in Table 2.
e damage degree of buildings affected by mining depends on the size of surface deformation and the ability of buildings to resist deformation. As there is no specification for this aspect in noncoal mines, the influence range of surface deformation of underground mining is delineated according to horizontal deformation ε > 2.0 mm/m or curvature k > 0.2 × 10 −3 or inclination i > 3.0 mm/m with Advances in Civil Engineering reference to the regulations for retaining coal pillars in buildings, water bodies, railways, and main shafts [29].
According to the calculation, the underground mining has little influence on the surface of the mine. e calculated subsidence value at village 4 is 0.075 m, which is close to the measured value of 0.07 m. e numerical simulation results are consistent with the actual situation. e regional subsidence and horizontal movement of the ground surface will not affect the surface buildings and structures. erefore, in the three regulations and other specifications, only three surface deformation indexes, such as horizontal deformation, curvature, and inclination, are used to reflect the influence degree of surface buildings and structures. Horizontal surface deformation can cause tension cracks in buildings, surface curvature will cause uneven deformation of buildings and structures, and surface tilt may cause construction buildings overturn, which affects the safety of surface buildings and structures. According to the calculation results, the maximum inclination, maximum curvature, and horizontal deformation of the mine are −0.5∼0.3 mm/m, K2 ore body -730m and below -590m high-grade ore Low-grade ore body K1 ore body -280m highgrade ore K1 ore body below -280 m and K2 ore body -730 m~-590 m high

Conclusion
In this paper, the interface program between CAD and finite element software based on surface spline interpolation developed by MATLAB software is used to realize the controllable mesh of complex stratum. Based on a mine engineering example, the mine numerical model including multiple complex stratum boundaries is established. e surface movement and deformation of underground mining are analyzed emphatically.
(1) Based on the self-developed interface program of CAD and finite element software, it can realize the controllable gridding of complex strata and provide a good foundation for the establishment of numerical calculation model. (2) According to the results of numerical simulation, the surface subsidence caused by underground mining is in good agreement with the measured value. e horizontal deformation, curvature, and inclination of the ground surface do not reach the critical value. e mining method of the filling method has no obvious effect on the surface structures. e filling mining can effectively control surface subsidence. Mining steps have an impact on surface subsidence. e results of the previous mining are affected by the next mining.
(3) e developed interface program overcomes the difficulty of 3D complex formation boundary modeling. e program has the advantages of controllable unit size and node number. e program can realize area interpolation and automatic expansion of data without contour. e program effectively avoids the appearance of local high-density grid on the stratum boundary. e program improves the calculation speed of the model. e research results can provide reference for similar numerical simulation.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper. Advances in Civil Engineering 13