Mechanical Model for Calculating Surface Movement Related to Open-Pit and Underground Caving Combined Mining

Abstract

Surface movement under the combined mining of open-pit and underground caving is the main problem affecting safe and efficient mining. Taking the combined mining of Dahongshan Iron Mine as an example, a mechanical model for calculating surface movement range was established to reveal the movement mechanism. The topography, the caving zone, the gravel filling of the caving zone and mining depth were considered, and the surface principal stress was calculated. Based on the mechanical model, the mechanical relationship between the subsidence angle and the surface cracking range was established and the definition of subsidence angle was supplemented: the subsidence angle is the angle measured from the horizontal of the moving boundary line, a straight line with zero displacement from the edge of the deepest extraction level to the surface and the largest surface cracking range. A trial calculation method for calculating the subsidence angle was proposed: by presetting different subsidence angles, the subsidence angle corresponding to the maximum surface prone cracking range calculated by mechanical model was found as the real subsidence angle. The subsidence angle calculated by the mechanical model was 60°, consistent with the observed subsidence angle of 64° in the Dahongshan Iron Mine, which verifies the reliability of the model.

1. Introduction

The combined mining method of open-pit and underground caving can be divided in two ways. One is the transformation from open-pit to underground caving mining, and the other is open-pit and underground caving combined mining [1,2,3,4,5]. The transformation from open-pit to underground mining is mainly a combined mining method in chronological order, and is the only way to extend the service life of mines and improve the recovery rate. In the early stage, open-pit mining is adopted, and in the later stage, underground mining is mainly used. At present, research on the transformation from open-pit to underground caving mining is relatively mature, including the maintenance of the open-pit slope, the management and control of ground pressure in the underground stope, the structural parameters of stope, the development and mining method in the transition period, and the thickness of overburden and roof pillar [6,7,8,9,10,11,12]. However, some mines are affected by the occurrence of ore bodies, and there can be many ore veins or ore belts in the deep and shallow at the same time, and they can be at a certain distance in the horizontal direction. In order to maximize the efficient mining of mineral resources and avoid a large amount of waste of resources caused by the destruction of the mineability of adjacent ore bodies due to the separate mining of an ore body, shallow ore bodies have been mined by the open-pit method and the deep ore bodies have been mined by the underground caving method through optimized design. Therefore, efficient open-pit and underground caving combined mining modes, mainly reflected in spatial coordination, have been formed, such as the combined mining of the upper and lower ore belts in the southeast of Gongchangling Iron Mine in China, and the combined mining of Dahongshan Iron Mine in Yuxi in China [5,13,14]. It is clearly stipulated in the “Safety regulation for metal and nonmetal mines” that, within the range of surface movement caused by underground caving mining, open-pit mining shall not be carried out at the same time unless effective technical measures are taken [15]. Therefore, the study of surface movement under the combined mining of open-pit and underground caving in metal mines has important practical significance for safe and efficient mining.

Research on the mechanisms of rock strata movement caused by coal mining is relatively mature, and many theoretical models with practical application value have been proposed, including the theory of voussoiring beams, the theory of transferring rock beams, the theory of key stratum and the theory of plate model [16,17,18,19,20]. However, these theoretical models are not suitable for describing the rock movement mechanisms in metal mines because of the great difference between the geological conditions and mining methods. In order to better understand the mechanisms and more accurately estimate the range of surface subsidence, researchers have also proposed some classic models and methods for the rock movement in metal mines. The arch theory has been widely recognized in the study of describing the caving mechanism of metal ore caving mining. It is believed that after the ore body is mined, the overlying rock mass continues to cave upward in an arched shape. Under the condition of large mining depth, the caving arch is in a natural balance state before the caving reaches the surface. After caving to the surface, due to the different mining and geological conditions, the surface has different degrees of deformation, fracture or no obvious deformation [21]. In addition to arch caving, rock caving in some mines also shows other forms, including cylindrical shapes and plunger shapes [22]. There are two types of “prism models” for the collapse of rock mass: straight line (Figure 1a) and broken line (Figure 1b) [23]. According to on-site mine observations, the broken-lines prism model is more appropriate for describing the fracture and shear slip processes of the rock mass. Based on the stress analysis of surrounding rock, the “arc plate model” considers that the boundary between the stress-reduced area and stress-increased area around the goaf is elliptical (Figure 1c). The stress-reduced area corresponds to the caving zone, and the stress-increased area corresponds to the fracture zone and the bending deformation zone. The model preliminarily considers the stress distribution and rock mass strength, and has been verified from centrifugal simulation tests, which is consistent with many mine cases [23]. Hoek regarded the hanging wall as a slope, and proposed a limit equilibrium method to predict the surface subsidence range, which has been widely applied [24,25]. However, most of these models are descriptive and fail to effectively reveal the essential mechanism of surface subsidence, and there are still significant discrepancies in predicting surface subsidence. In addition, these models fail to effectively consider the impact of open-pit mining on surface subsidence under combined mining.

With the development of numerical simulation technology, scholars have also adopted different numerical methods to study the strata movement and surface subsidence caused by underground mining [26]. Li et al. [27] reproduced the entire process of strata movement and surface deformation caused by underground caving mining using the RFPA numerical method. Zhang et al. [28] further analyzed the fracture development characteristics during rock mass movement using RFPA. Parmar et al. [29] explored the impact of ground surface subsidence due to underground mining on surface infrastructure. Takhanov et al. [30] analyzed the possible zones of the mass deformation under combined mining based on the Phase2 software. Imashev et al. [31] studied the stress–strain state of the mass under an open pit with a change in the open-pit bottom width using a numerical simulation method. Ghazdali et al. [32] used RS2 software to analyze the stability of the underground mining method used in shallow-dip vein deposits hosted in poor quality rock. It is generally believed that numerical simulation methods are low-cost and time-saving. However, in order to obtain reliable calculation results, numerical calculations must be combined with on-site monitoring data to repeatedly verify and correct parameters. Thus, it inevitably takes a long time to ensure prediction accuracy.

Many research works and engineering practices have shown that topography is an major factor affecting surface movement range caused by underground caving mining in metal mines. The range of surface movement is often different under different topography features [23,33,34]. For open-pit and underground caving combined mining, under static conditions, continuous open-pit mining essentially affects the range of surface movement by changing the original topography features. Thus, in order to reveal the surface movement mechanism under the open-pit and underground caving combined mining, it is necessary to establish a mechanical model that can consider the topography. In this study, taking Dahongshan Iron Mine as the background, a mechanical model was established to calculate the surface subsidence range under combined mining. The change of topography caused by open-pit mining, the shape of the caving zone, the gravel filling of caving zone and mining depth were considered in the model. The definition of the subsidence angle was supplemented, and the calculating method was proposed by the mechanical model. The reliability of the mechanical model was verified by the field measured data, which can provide a new theoretical method for calculating surface movement range for metal mines used for open-pit and underground caving combined mining.

2. Engineering Background

Dahongshan Iron Mine is located in Yuxi City, Yunnan Province, China. There are two ore veins in the mining area, shallow and deep, which are mined simultaneously by open-pit mining and underground block caving, separately. The strike of the underground ore body is almost east–west. The shape of the underground ore body is high in the east and low in the west, thick in the middle, thin at the edges, and upturned in the north and south, generally resembling a ship shape. The deep ore vein is 1969 m from east to west, 440~640 m wide from north to south and 362.48~988.31 m deep, and the occurrence elevation is 25.72~945.00 m. The occurrence elevation of the shallow ore vein is 680~1012 m, the strike is northwest, the dip is southwest, and the average thickness is 11~18 m. The shallow vein is open-pit mined by the process of “steep slope stripping and slow slope mining”. The upper dimension of the stope is 1130 m long and 724 m wide, and the bottom dimension is 523 m long and 172 m wide. The production capacity of open-pit mining after optimized design is 5.5 million t/a. The production capacity of the underground mining is 4 million t/a, and it is divided into four sections according to the elevation, namely the head section, the middle section, the lower section and the deep section. Except for the head section, which is mined by the open-stope method, the other sections are mined by block caving, as shown in Figure 2.

The main problem restricting safe and efficient mining in Dahongshan Iron Mine is the mutual influence between open-pit and underground caving mining. The minimum plane distance between the open-pit final boundary and the underground main mining area is only about 250 m. With the continuous increase of the caving scope and the depth of the underground ore body, the range of surface subsidence will continue to expand, which will have an important impact on the stability of the open-pit slope. In addition, since the grade of the deep ore body is twice that of the shallow ore body, it is difficult to slow down the underground mining plan to ensure the productivity and economic benefits. However, there are still about 30 million tons of ore bodies in a shallow ore vein. The shutdown of the open pit would cause a huge waste of resources. Therefore, determining the range of surface movement is the key to evaluate whether open-pit and underground caving combined mining can be used.

3. Characteristics and Angle of Rock Movement

The rock strata movement process caused by metal mine caving mining can be summarized as follows: with underground caving mining, when the span and volume of the goaf reach the limit, the overlying rock mass will collapse. Under continuous mining, the caving of the rock mass will develop to the surface, causing surface collapse and deformation. Cheng et al. [35] divided the surface deformation process into two temporal stages, the chimney caving development stage and the post-chimney deformation stage, based on the type of driving force. Van As [36] subdivided the mining-induced surface influence area into a caved zone, fractured zone and continuous subsidence zone. Combining these two methods, the division of surface movement caused by underground caving mining is shown in Figure 3. The surface influence area is also generally divided into caved zone, fractured zone and continuous subsidence zone. However, the caved zone is subdivided into the initial caved zone and the caving extension zone. When the rock caving reaches the surface, the initial caved zone will be formed. After caving to the surface, due to the loss of support, the surrounding rock mass will fully collapse to form the caving extension zone. In the surface within a certain range from the edge of the caving extension zone, only different lengths of tension cracks occur, forming the fractured zone. No cracking occurs for the surface within a certain range outside the fractured zone, but it presents continuous subsidence, which is the continuous subsidence zone.

According to Figure 3, the initial caving angle is the angle measured from the horizontal of the line extending from the edge of the deepest extraction level to the position where the initial caved zone developed to the surface, which is set as α. The second caving angle is the angle measured from the horizontal of the line extending from the edge of the deepest extraction level to the edge of the caving extension zone, which is set as α₁. The fracture initiation angle is the angle measured from the horizontal of the line extending from the edge of the deepest extraction level to the edge of the fractured zone, which is set as α₂. The subsidence angle is the angle measured from the horizontal of the line extending from the edge of the deepest extraction level to the edge of the continuous subsidence zone, which is set as α₃, as shown in Figure 3.

4. Establishment of the Mechanical Model

Caving mining of the lower section (elevation 400~500 m) in the deep ore vein has the greatest impact on the surface subsidence in Dahongshan Iron Mine. The lower section is divided into several subsections with a mining height of 20 m, named 480, 460, 440, 420 and 400 subsections according to the elevation of the floor of each subsection, and caving mining is carried out from top to bottom successively. The initial mining subsection was the 480 subsection. By the end of 2013, the 420 subsection had been mined, and then the 400 subsection was mainly mined, as shown in Figure 4. In reference [14], the shape of the initial caved zone is determined based on the field monitoring data; it is roughly “arched”, as shown in Figure 2.

The open-pit in Dahongshan Iron Mine is a hillside open pit. Among the geological profiles of Dahongshan Iron Mine, the A38 profile passes through both the open-pit and the center of the underground caving mining at the same time, and can be used as the main profile to determine the surface movement. Therefore, the A38 exploration profile was selected to establish a two-dimensional mechanical model, as shown in Figure 5. The average slope angle on the A38 section is β = 45°, the elevation of the bottom of the slope is about 800 m, and the top of the slope is a near-horizontal topography with an elevation of 1170 m. The initial caving point C on the surface is basically located in the middle of the surface projection of the 420 subsection. For the convenience of calculation on the A38 profile, it is assumed that the shape of the initial caved zone is approximately an isosceles triangle. The horizontal distance from point C to point D is 460 m, and the maximum elevation of the initial mining subsection (480 subsection) is 500 m. The model adopts the rectangular coordinate system, and takes the boundary point A of the 420 subsection as the origin. The span of the 420 subsection is 390 m (2l = 390 m), and the depth is 750 m (H₂ = 750 m). Thus, the initial caving angle (α) is 75°.

With the mining of the ore body at the elevation of 420~500 m, the overburden continued to fall, and the caving gravel continued to fill the initial caved zone. After the goaf was filled with the caving gravel, the initial caved zone was not fully filled, and the unfilled space (CFG) is shown in Figure 5. In the model, the coordinate of point F at the maximum elevation of gravel filling in the initial caved zone can be estimated from the initial caved zone volume and the mined ore volume based on triangle similarity. Assuming that the initial caved zone is a cone formed by the rotation of the isosceles triangle ACI around the central axis, the total volume of the filled gravel can be set as V₁ in consideration of the crushing expansion of the rock mass. The mined volume of the ore body was set as V₂. The volume of the cone formed by rotating the isosceles triangle BCH around the central axis in Figure 5 was set as V₃. Assuming BH = 2R and FG = 2r, the following similarity relationship can be obtained.

$$\frac{V_3-\left(V_1-V_2\right)}{V_3}=\frac{\frac{1}{3}{r}^2\left(H-H_1\right)\pi }{\frac{1}{3}{R}^{2}H\pi }=\frac{r^3}{R^3}$$

(1)

It is assumed that the goaf is completely filled by caving gravel. Since Formula (1) does not involve the size of the goaf, it is applicable to any irregular goaf formed by caving mining. According to the field data, V₂ = 5,790,000 m³. Thus, R, V₁ and V₃ can be calculated with the following formula.

$$\{\begin{array}{l}R=l-\frac{500-420}{\tan \alpha }\\ V_1=\frac{1}{3}{R}^{2}H\pi K_{\mathrm{s}}\\ V_3=\frac{1}{3}{R}^{2}H\pi \end{array}$$

(2)

where K_s is the crushing expansion coefficient (K_s = 1.15). Thus, the coordinate of point F is calculated to be (108, 400); that is, the filling height of the initial caved zone is 400 m.

In Figure 5, the rock mass outside the moving boundary line (AE) will not move, so the AE line can be set as a fixed boundary. Therefore, considering the gravity G, the normal support force P, the topography (CDE), the moving boundary line (AE), and the initial caved line (AC), a mechanical model was established with the quadrangle (ACDE) as the research object, as shown in Figure 6. In the mechanical model, the moving boundary line was preset and the moving angle α₃ was unknown. The equations of straight lines AC, CD and DE, the coordinates of point E and the normal support force P can be expressed as [37]:

$$\{\begin{array}{l}{x}_E=-\frac{1015}{1+\tan {\alpha }_3}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{x}_{\mathrm{E}}\le -265\\ y_{\mathrm{E}}=\frac{1015\times \tan {\alpha }_3}{1+\tan {\alpha }_3}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{x}_{\mathrm{E}}\le -265\\ x_E=-\frac{750}{\tan {\alpha }_3}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{x}_{\mathrm{E}}\ge -265\\ y_E=750\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{x}_{\mathrm{E}}\ge -265\\ y_{\left(x\right)\mathrm{AC}}=\tan \alpha \times x=3.73x\\ y_{\left(x\right)\mathrm{CD}}=750\\ y_{\left(x\right)\mathrm{DE}}=x+1015\\ P={\gamma }_1\left(y_{\mathrm{F}}-y\right)\left(\sin \alpha \cot \theta -\cos \alpha \right)={\gamma }_1\left(400-3.73x\right)\left(0.97\cot \theta -0.26\right)\hspace{1em}\hspace{1em}0\le x\le 108\end{array}$$

(3)

where γ₁ is the average bulk density of the filling gravel, and θ is the maximum slope angle of the gravel bulk structure. From field experience, θ ≈ 55°.

5. The Solution of the Mechanical Model

5.1. Coordinate Transformation

The mechanical model (Figure 6) can be regarded as a cantilevered beam model with variable cross-sections placed at an angle (α₃). To facilitate the solution, the rectangular coordinate system was rotated counterclockwise around point A by an angle (α₄), and α₄ = π/2 − α₃, as shown in Figure 7. The coordinates of each point after transformation were obtained by the following formula:

$$\{\begin{array}{l}{x}_1=x\cos {\alpha }_4+y\sin {\alpha }_4\\ y_1=y\cos {\alpha }_4-x\sin {\alpha }_4\end{array}$$

(4)

where x₁ is the abscissa of any point after transformation, y₁ is the ordinate of any point after transformation, x is the original abscissa of any point, and y is the original ordinate of any point.

The equations of straight lines and the normal support force P₁ after coordinate transformation are expressed as:

$$\{\begin{array}{l}{y}_{\left(x_1\right)\mathrm{AC}1}=\tan \left(\alpha -{\alpha }_4\right)\times x_1\\ y_{\left(x_1\right)\mathrm{CD}1}=-\tan {\alpha }_4\times x_1+750\left(\cos {\alpha }_4+\sin {\alpha }_4\tan {\alpha }_4\right)\\ y_{\left(x_1\right)\mathrm{DE}1}=\tan \left(\beta -{\alpha }_4\right)\times \left(x_1+265\cos {\alpha }_4-750\sin {\alpha }_4\right)+750\cos {\alpha }_4+265\sin {\alpha }_4\\ P_1={\gamma }_1\left\{400-\left[\sin {\alpha }_4+\tan \left(\alpha -{\alpha }_4\right)\cos {\alpha }_4\right]x_1\right\}\left(\sin \alpha \cot \theta -\cos \alpha \right)\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}0\le x_1\le x_{\mathrm{F}1}\end{array}$$

(5)

5.2. Piecewise Solution

The mechanical model (Figure 7) can be regarded as a cantilevered beam model with variable cross-sections. According to the mechanics of materials, the stress distribution on any section parallel to AE can be solved in three parts. The first part is the FC-segment, which cannot act towards the normal support force P₁; that is, x_F1 ≤ x₁ ≤ x_C1. The second part is the DF-segment that acted by P₁; that is, x_D1 ≤ x₁ < x_F1. The third part is the ED-segment that acted by P₁; that is, 0 ≤ x₁ < x_D1.

5.2.1. Solving the First Part (x_F1 ≤ x₁ ≤ x_C1)

When x_F1 ≤ x₁ ≤ x_C1, the shear force F_s generated by gravity G on any section parallel to AE can be calculated by the following formula:

$$F_{\mathrm{s}}={\displaystyle {\int }_{x_1}^{x_{\mathrm{C}1}}\left(y_{\left(x_1\right)\mathrm{CD}1}-y_{\left(x_1\right)\mathrm{AC}1}\right)}\mathrm{d}{x}_1\times {\gamma }_2\cos {\alpha }_4$$

(6)

where γ₂ is the average bulk density of the rock mass.

The mechanical model is a two-dimensional mechanical model with unit thickness, and any section is a rectangle with thickness b = 1. Thus, the shear stress τ_x1y1 can be calculated.

$$\{\begin{array}{l}{\tau }_{x_1{y}_1}=\frac{6F_{\mathrm{s}}\left(\frac{h_2^2}{4}-y_1^2\right)}{b{h}_2^3}=6F_{\mathrm{s}}\left[\frac{1}{4h_2}-\frac{y_1^2}{h_2^3}\right]\\ h_2=y_{\left(x_1\right)\mathrm{CD}1}-y_{\left(x_1\right)\mathrm{AC}1}\end{array}$$

(7)

where h₂ is the height of the rectangular section.

The axial force F_n generated by gravity G on any section parallel to AE can be calculated.

$$F_{\mathrm{n}}={\displaystyle {\int }_{x_1}^{x_{\mathrm{C}1}}\left(y_{\left(x_1\right)\mathrm{CD}1}-y_{\left(x_1\right)\mathrm{AC}1}\right)}\mathrm{d}{x}_1\times {\gamma }_2\sin {\alpha }_4$$

(8)

The bending moment M generated by gravity G on any section parallel to AE can be calculated.

$$M=F_{\mathrm{s}}\left(x_{\mathrm{G}}-x_1\right)={\displaystyle {\int }_{x_1}^{x_{\mathrm{C}1}}\left(y_{\left(x_1\right)\mathrm{CD}1}-y_{\left(x_1\right)\mathrm{AC}1}\right)}\mathrm{d}{x}_1\times {\gamma }_2\cos {\alpha }_4\left(x_{\mathrm{G}}-x_1\right)$$

(9)

where x_G is the abscissa of the gravity center of the right rock mass cut by calculated section. For the triangular, it can be calculated as:

$$x_{\mathrm{G}}=\frac{2x_1+x_{\mathrm{C}1}}{3}$$

(10)

Thus, the normal stress σ_x1 on the section is as follows:

$${\sigma }_{x_1}=\frac{F_{\mathrm{n}}}{h_2}+\frac{12M}{b{h}_2^3}\left(y_1-\frac{y_{\left(x_1\right)\mathrm{CD}1}+y_{\left(x_1\right)\mathrm{AC}1}}{2}\right)=\frac{F_{\mathrm{n}}}{h_2}+\frac{12M}{h_2^3}\left(y_1-\frac{y_{\left(x_1\right)\mathrm{CD}1}+y_{\left(x_1\right)\mathrm{AC}1}}{2}\right)$$

(11)

According to the “Code for the design of concrete structures” [38], the continuous variable cross-section beam can use material mechanics to calculate the principal tensile stress.

$${\sigma }_{{}_{1,3}}=\frac{{\sigma }_{x_1}+{\sigma }_{y_1}}{2}\pm \sqrt{{\left(\frac{{\sigma }_{x_1}-{\sigma }_{y_1}}{2}\right)}^2+{\tau }_{x_1{y}_1}{}^2}$$

(12)

where σ_y1 is the longitudinal normal stress. According to the assumption of no extrusion of longitudinal fibers in material mechanics, σ_y1 = 0.

In order to evaluate the surface cracking, it is necessary to calculate the principal stress σ_1CD on the surface (CD). It is specified that tension is positive and pressure is negative. Thus, σ_1CD can be calculated.

$${\sigma }_{{}_{1\mathrm{CD}}}={\sigma }_{x_1\mathrm{CD}}=-\frac{\left|F_{\mathrm{n}}\right|}{h_2}+\frac{6\left|M\right|}{h_2^2}$$

(13)

Surface cracking is usually manifested as tensile failure. Therefore, the place where σ_1CD is tensile stress and exceeds the tensile strength of rock mass σ_t, namely |σ_1CD| > |σ_t|, will crack.

5.2.2. Solving the Second Part (x_D1 ≤ x₁ < x_F1)

When x_D1 ≤ x₁ < x_F1, the shear force F_s2 generated by the normal support force P₁ on any section parallel to AE can be calculated.

$$F_{\mathrm{s}2}={\displaystyle {\int }_{x_1}^{x_{\mathrm{F}1}}\left(P_1\right)}\mathrm{d}{x}_1\times \cos \left(\alpha -{\alpha }_4\right)$$

(14)

The axial force F_n2 generated by P₁ on any section parallel to AE can be calculated.

$$F_{\mathrm{n}2}={\displaystyle {\int }_{x_1}^{x_{\mathrm{F}1}}\left(P_1\right)}\mathrm{d}{x}_1\times \sin \left(\alpha -{\alpha }_4\right)$$

(15)

From the resultant moment theorem of linear load, the bending moment M₂ generated by P₁ on any section parallel to AE can be obtained.

$$M_2=F_{\mathrm{s}2}\frac{x_{\mathrm{F}1}-x_1}{3}$$

(16)

Similarly, when x_D1 ≤ x₁ < x_F1, the principal stress σ_1CD2 on the surface (CD) can be calculated.

$${\sigma }_{{}_{1\mathrm{CD}2}}=-\frac{\left|F_{\mathrm{n}}\right|}{h_2}+\frac{6\left|M\right|}{h_2^2}-\frac{6\left|M_2\right|}{h_2^2}-\frac{\left|F_{\mathrm{n}2}\right|}{h_2}$$

(17)

5.2.3. Solving the Third Part (0 ≤ x₁ < x_D1)

When 0 ≤ x₁ < x_D1, the shear force F_s generated by gravity G on any section parallel to AE can be calculated.

$$F_{\mathrm{s}}={\displaystyle {\int }_{x_{\mathrm{D}1}}^{x_{\mathrm{C}1}}\left(y_{\left(x_1\right)\mathrm{CD}1}-y_{\left(x_1\right)\mathrm{AC}1}\right)}\mathrm{d}{x}_1\times {\gamma }_2\cos {\alpha }_4+{\displaystyle {\int }_{x_1}^{x_{\mathrm{D}1}}\left(y_{\left(x_1\right)\mathrm{DE}1}-y_{\left(x_1\right)\mathrm{AC}1}\right)}\mathrm{d}{x}_1\times {\gamma }_2\cos {\alpha }_4$$

(18)

The axial force F_n generated by gravity G on any section parallel to AE can be calculated.

$$F_{\mathrm{n}}={\displaystyle {\int }_{x_{\mathrm{D}1}}^{x_{\mathrm{C}1}}\left(y_{\left(x_1\right)\mathrm{CD}1}-y_{\left(x_1\right)\mathrm{AC}1}\right)}\mathrm{d}{x}_1\times {\gamma }_2\sin {\alpha }_4+{\displaystyle {\int }_{x_1}^{x_{\mathrm{D}1}}\left(y_{\left(x_1\right)\mathrm{DE}1}-y_{\left(x_1\right)\mathrm{AC}1}\right)}\mathrm{d}{x}_1\times {\gamma }_2\sin {\alpha }_4$$

(19)

The bending moment M generated by gravity G on any section parallel to AE can be calculated.

$$M=F_{\mathrm{s}}\left(x_{\mathrm{G}}-x_1\right)$$

(20)

where x_G can be calculated by the following formula under the quadrilateral condition.

$$x_{\mathrm{G}}=\frac{2x_1+x_{\mathrm{D}1}+x_{\mathrm{C}1}}{4}$$

(21)

When 0 ≤ x₁ < x_D1, the shear force, axial force and bending moment generated by P₁ on any section parallel to AE can be calculated by Formulas (14)–(16), respectively. Additionally, the principal stress σ_1DE on the surface (DE) can be calculated by Formula (17).

In the above calculation of surface tensile stress, the gravity, the topography, initial caved zone and the filling lateral support force were considered. In addition, the influence of mining depth and span could also be considered by geometric boundary Equation (3) in the model. Although the calculation of the variable section deep cantilever beam using material mechanics will produce an error, it ensures the calculability of the model and facilitates engineering application. According to mechanical research [39], the maximum error range of the stress solution of cantilever deep beam calculated by material mechanics is about 7%~14%. For mining engineering, such an error is acceptable.

6. The Solution of the Subsidence Angle

6.1. Supplementary Definition and Calculation Method of the Subsidence Angle

In the traditional definition, the subsidence angle and the fracture initiation angle are only determined according to the characteristics of surface movement, and they are independent of each other. Previous studies have shown that in the development time of surface movement, the surface cracking will go through the continuous deformation stage [35]. In terms of development space, the fractured zone is adjacent to the continuous subsidence zone (Figure 3). This indicates that the subsidence angle is related to the surface cracking range. In the mechanical model, the subsidence angle (α₃) is unknown. Thus, the surface principal stress and cracking range obtained by taking different subsidence angles are different. In fact, the moving boundary line determined by the subsidence angle is not only a fixed boundary with a displacement of 0, but also a supporting boundary for rock movement under the action of gravity.

Two exceptional cases can be analyzed by the mechanical model (α₃ = 0° and α₃ = 90°), as shown in Figure 8. When α₃ = 0°, the gravity action line is perpendicular to the moving boundary line, so the principal stress on the surface is all compressive stress, without cracking. When α₃ = 90°, the gravity action line is parallel to the moving boundary line, so tensile stress will only occur in the surface corresponding to the underground caving mining range, resulting in cracking. These two kinds of surface cracking are not consistent with the actual situation of most sites, and the surface cracking range often exceeds the underground mining scope [33]. Therefore, between 0° and 90°, there must be a critical subsidence angle which makes the surface cracking range the largest. The critical subsidence angle should be the true subsidence angle for the actual situation on site. Therefore, based on the mechanical model, the supplementary definition of the subsidence angle is that the subsidence angle is the angle measured from the horizontal of the moving boundary line, a straight line with zero displacement from the edge of the deepest extraction level to the surface and the largest surface cracking range. In this definition, the surface cracking range is the horizontal span of the fractured zone.

Based on the supplementary definition of the subsidence angle, the actual subsidence angle can be calculated by the mechanical model. Because the calculation formula for directly deriving the subsidence angle has a large amount of calculation and many high-order terms will appear, the trial calculation method can be used. First, by presetting different subsidence angles such as α₃ = 30°, α₃ = 45°, α₃ = 60° and α₃ = 75°, different mechanical models are established. Then, the surface tensile stress zone of each model is calculated, and the surface prone cracking range is determined by the tensile strength of rock mass. Finally, by comparing the surface prone cracking range of each model, the preset subsidence angle corresponding to the maximum surface prone cracking range is found as the real subsidence angle.

6.2. Calculation Result of the Subsidence Angle

Under the open-pit and underground caving combined mining, the impact of surface movement on the stability of the slope is the main problem, so point J is the boundary point of the mechanical model in Figure 5. Therefore, when the trial calculation method is used, the preset minimum subsidence angle α₃ is 30° (Figure 5). To ensure the calculation accuracy, the preset subsidence angle of each model is increased by 10° (α₃ = 30°, 40°, 50°, 60°, 70°, 80°). The parameters required for the calculation are indicated in

Table 1. The parameters of the mechanical model under open-pit and underground caving combined mining in the Dahongshan Iron Mine.

Average Bulk Density of the Filling Gravel (γ1)	Maximum Slope Angle of the Gravel Bulk Structure (θ)	Average Bulk Density of the Surrounding Rock (γ2)	Initial Caving Angle (α)
16.8 kN/m3	55°	27.5 kN/m3	75°

.

The trial calculation results of the surface principal stress distribution of mechanical models with different preset subsidence angles are shown in Figure 9. The surface principal stress increases from point C to point D. The principal stress will jump at point D, the compressive stress will increase significantly (Figure 9a,b), and the tensile stress will suddenly decrease (Figure 9c,d). From point D to point E, the compressive stress decreases rapidly and the tensile stress increases rapidly. When α₃ = 30°, the tensile stress zone is distributed on the slope. When α₃ = 40°, the surface principal stress changes into tensile stress, and the tensile stress area on the slope becomes larger, while the compressive stress area decreases. When α₃ ≥ 50°, the surface and slope are all in tension. The average tensile strength of the rock mass is σ_t = 0.6 MPa, from the Hoek–Brown criterion [40]. The critical point with the tensile stress of 0.6 MPa was set as the cracking point K. With the increase of α₃, the boundary point E continued to move to C, and the cracking point K also moved to C. The movement speed of point E was higher than that of point C, so the horizontal distance between two points changed from 654 m (α₃ = 30°) to 55 m (α₃ = 80°). The surface prone cracking range L₂ under different α₃ is calculated in Figure 9. When α₃ = 50°, the maximum surface prone cracking range is 497 m. When α₃ = 60°, the surface prone cracking range is 487 m, which is close to the maximum value—only 10 m smaller. According to the supplementary definition of the subsidence angle, the true subsidence angle is about 50°. However, compared with the subsidence angle of 50°, the moving boundary line of α₃ = 60° can produce effective support for the rock movement and prevent the subsidence of the rock mass outside the boundary line. Thus, the true subsidence angle should be taken to be 60°.

7. Analysis of On-Site Monitoring Data

In document [14], the initial caved zone in Dahongshan Iron Mine was determined according to the monitoring data (drilling data, laneway observation data and microseismic monitoring data). Its spatial shape was roughly a cone (Figure 10), consistent with the assumptions of the mechanical model. In addition, the microseismic monitoring system and the surface total station displacement monitoring system were arranged. Additionally, the development range of surface cracking was regularly observed, annually. The surface cracking range over the years and the cumulative surface subsidence distribution in 2021 on the A38 profile are shown in Figure 10.

In order to analyze the impact of surface movement on the open-pit slope and to maintain consistency with the mechanical model, the rock mass on the same side as the open pit was taken as the research object. Microseismic data reflecting rock mass damage near the A38 profile in 2012 and 2013 were extracted, respectively, as shown in Figure 11. There were frequent microseismic events in the rock mass at one side of the open pit, and there was obvious microseismic concentration in some parts. This shows that this part of the rock mass was in the influence range of combined mining and had damage and fracture. The rock mass with many microseismic events was within the range delineated by the subsidence angle of 60°. In addition, it can be seen that since the crack initiation point appeared on the surface in 2012, the surface cracking range continued to expand to the outside. Then, the surface cracking range increased by about 184 m on one side of the open pit in 2013, with the largest expansion seen (Figure 10). Until 2016, the cracking range only slightly expanded by 10 m on one side of the open pit, indicating that the expansion of surface cracking was sufficient in 2013. From 2013 to 2021, the 400 subsection was mined. The mining depth of the 400 subsection increased by 20 m, but the mining area was smaller than the 420 subsection (Figure 10). However, on the open-pit side, the surface fractured zone in 2021 remained within the scope of 2016 and had not expanded. This shows that when the mining area does not expand outward, the increase of the mining depth has little impact on the surface movement. The fracture initiation angle was maintained at 75° (Figure 10). This also shows that from 2013 to 2021, the surface deformation was always in the state of full expansion. From the surface subsidence distribution curve of the A38 profile in 2021 (Figure 10), the subsidence value in the center of the caved zone was the largest, reaching 1.1 m. In the continuous subsidence zone, the subsidence value at the 1150 bench of the slope was 0.1 m, and the subsidence value at the 1120 bench was close to 0.0 m. Therefore, in the A38 profile, the measured subsidence angle was 64°. This is consistent with the actual subsidence angle calculated by the mechanical model, which verifies the reliability of the established mechanical model.

8. Conclusions

(1)

Taking the combined mining of the Dahongshan Iron Mine as an example, a mechanical model for calculating the surface movement range was established according to the zoning characteristics of surface movement to reveal the movement mechanism. The change in topography caused by open-pit mining, the shape of the caving zone, the gravel filling of the caving zone and the mining depth were considered in the model.

(2)

Based on the mechanical model, the definition of the subsidence angle was supplemented: the subsidence angle is the angle measured from the horizontal of the moving boundary line, a straight line with zero displacement from the edge of the deepest extraction level to the surface and the largest surface cracking range.

(3)

According to the supplementary definition of the subsidence angle, a trial calculation method for calculating the subsidence angle was proposed: by presetting different subsidence angles, the subsidence angle corresponding to the maximum surface prone cracking range calculated by mechanical model can be found as the real subsidence angle.

(4)

The subsidence angle calculated by the mechanical model was 60°, consistent with the observed subsidence angle of 64° in Dahongshan Iron Mine. The mechanical model provides a new method for calculating the surface movement range for metal mines used for open-pit and underground caving combined mining.