Abstract

The interface program of finite element software based on surface spline interpolation is developed by using MATLAB software. The 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. The influence of underground mining on surface movement and deformation under complex stratum conditions by using the FLAC3D software further was discussed. The 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. The calculated subsidence value is in good agreement with the measured value. The values of surface tilt, surface curvature, and surface horizontal deformation are less than the relevant regulations. The mining method of the filling method has no obvious effect on the surface structures. The research results can provide reference for similar numerical simulation.

1. Introduction

With the rapid development of economy, people’s demand for mineral products is increasing, which makes the mining industry develop vigorously. The surface subsidence disaster brings hidden danger to people’s safety [13]. The underground mining can cause deformation and damage to surrounding rock mass, surface structures, railways, and water bodies. Therefore, the surface movement and deformation caused by underground mining and the damage of surface buildings have become one of the global man-made environmental disasters and a hot issue in the mining field [46].

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. They 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. The 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. Therefore, 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]. The influence of dynamic load, groundwater, and other comprehensive factors on the surface movement and deformation can be considered in the numerical simulation [19]. Therefore, the numerical simulation has become a common analysis method [2023]. Hadi et al. [24] used the FLAC3D software to simulate the ground surface subsidence due to sublevel caving method. Then, 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. The 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. The 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. The 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. The relevant horizontal and vertical sections, geological boundary contour map, ore body distribution map, and other data are used in the model. Through the field engineering geological survey and rock mechanics, the mine numerical model is established. The actual rock mass mechanical parameters are obtained from the test. Then, 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. The surface subsidence is in good agreement with the measured values. The calculation results of surface movement and deformation can provide reliable guidance for mining.

2. Mesh Construction Method of 3D Complex Surface

The MATLAB program is used to realize the interface between CAD and finite element software. The 3D complex surface model can be directly generated from CAD data, which can be used by finite element software.

The features of the interface program are as follows:(1)The interface program can automatically extract the coordinates of isoline points and corresponding elevation values and coexist in the text file.(2)The interface program can automatically set the horizontal and vertical dimensions of the 3D surface mesh. The difference of the grid near the actual coordinate point can be found automatically. Finally, the regular and uniform 3D surface is obtained.(3)The interface program can realize the automatic extension of isoline to nonisoline data area. The 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. The coordinates and corresponding elevation values of the points before encryption arewhere Xi and Yi are the X and Y coordinates of elevation point I, respectively, and Hi is the elevation value of point i.

Formula (1) is a bivariate single-valued list function, which is fitted. The expression of bivariate spline function is as follows:where 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. The coefficient to be calculated can be determined by the following formula:where and is the weighting coefficient of the node, which can be taken as 0. Equation (3) is expressed as matrixwhere Cm×1 is a symmetric matrix composed of node coordinate values and weighted coefficients. If Am×m is not a singular matrix, then the equation can be solved and the coefficient matrix is

Then, 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

The partial derivative of elevation function of isoline at any grid node is expressed as

c3+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.

Through formula (7) and grid node coordinates, the corresponding coordinates of the mesh nodes of the fitting surface can be obtained. The 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. The local part of CAD contour line (such as the lower left corner of CAD drawing) cannot meet the grid density. The program carries out automatic interpolation calculation to generate the 3D surface model. The different grayscale in Figure 2(a) represents different elevation. The model is consistent with the isoline in Figure 1.

3. Engineering Case Analysis

3.1. General Situation of Mining Area

The surface of the mining area is gentle, slightly inclined to the east, with a gradient of about 7‰. The surface of the mining area is approximately horizontal. The 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. There are mainly two ore belts K1 and K2 in the mining area. The K1 ore belt occurs in the top alteration of the rock mass. The K2 mineralization belt occurs in the alteration of the rock bottom. The 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. The 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. This paper will focus on the analysis of the surface movement and deformation of underground mining. The geological boundary of strata in the mining area is complex. The 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. Therefore, it is very practical to use numerical simulation analysis for the mine with complex stratum boundary.

3.2. 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. The 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 horizontal projection of ore body, the high- and low-grade ore bodies are divided. The three-dimensional model of the ore body and surrounding rock three-dimensional model is shown in Figure 7. The different colors represent different grades of ore bodies. The two-layer ore bodies are composed of four three-dimensional curved surfaces. The 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]. This time, the load is applied according to the self-weight stress of the rock mass. The z-direction constraint is applied at the bottom of the model, and the X-direction and Y-direction constraints are applied around the model. The surface is a free surface.

3.3. Numerical Simulation Scheme and Parameter Value

This 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. The filling is carried out at the end of each mining step. The 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. The typical Hoek–Brown calculation curve is shown in Figure 10. The final mechanical parameters of rock mass are shown in Table 1.

3.4. Analysis of Numerical Simulation Results

After the high-grade 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. The 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. The 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. The surface displacement reaches 0.2 m. The mining of all low-grade ore bodies has a significant impact on strata movement. The 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. The cloud chart of strata movement in different mining stages is shown in Figure 11. The results of FLAC3D are imported into surfer software, and the final deformation of the surface is shown in Figure 12.

3.5. Analysis of Theoretical Calculation Results

At present, the mine has arranged the surface subsidence observation points along the 1-1 survey line. The mining of the lower ore body of Miaozhuang has been completed, and the mining area is far away from here. The measured results show that the surface subsidence has become stable. The maximum measured value of subsidence near village 4 is 0.07 m, which can be regarded as the final subsidence value of this point. There 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.

The formula for calculating the inclination of the earth’s surface in x and y directions iswhere is the surface subsidence value.

The formula for calculating the curvature of the earth’s surface in x and y directions iswhere i (x, y) is the surface inclination value in the x and y directions.

The formula for horizontal deformation in the direction of x and y of the surface iswhere U (x, y) is the horizontal surface movement in x and y directions.

The 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. The 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 Figures 1517. The final surface movement deformation value obtained is shown in Table 2.

The 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 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. The calculated subsidence value at village 4 is 0.075 m, which is close to the measured value of 0.07 m. The numerical simulation results are consistent with the actual situation. The regional subsidence and horizontal movement of the ground surface will not affect the surface buildings and structures. Therefore, 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, −1.6∼4.0 × 10−6/m, and −0.38∼0.18 mm/m, which are not within the scope of influence specified in the three regulations. The mining method of the filling method has no obvious impact on the surface structures.

4. 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. The 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. The horizontal deformation, curvature, and inclination of the ground surface do not reach the critical value. The mining method of the filling method has no obvious effect on the surface structures. The filling mining can effectively control surface subsidence. Mining steps have an impact on surface subsidence. The results of the previous mining are affected by the next mining.(3)The developed interface program overcomes the difficulty of 3D complex formation boundary modeling. The program has the advantages of controllable unit size and node number. The program can realize area interpolation and automatic expansion of data without contour. The program effectively avoids the appearance of local high-density grid on the stratum boundary. The program improves the calculation speed of the model. The research results can provide reference for similar numerical simulation.

Data Availability

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

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

The authors would like to thank Shijiazhuang Tiedao University for the support for the project research funding.