At present, the coal seam exploited in large open-pit mines in China is mostly near-horizontal coal seam, and the internal dumping technology in the mining area is widely adopted (Chen and Wu 2016; Wang et al. 2019; Wu et al. 2020; Zhao et al. 2020; Lelek and Kulczycka 2021). In the traditional open-pit mining design, the transportation lines discharged by mining haul trucks in these large-scale open-pit mines are set on the two ends of the wall, and the formed end slope angle is less than the stable slope angle. As a result, a large number of coal resources are occupied by the end slope (Ross et al. 2019; Kuznetsova and Anfyorov 2019). With the continuous advancement of the inner dump and working slope, the coal occupied by the end slope will be re-buried by the inner dump. Due to the limitations of the current mining technology and mining equipment in China, it is difficult to conduct secondary mining in these mines (Elmouttie and Karekal 2017; Pan et al. 2021; Li et al. 2020).

Generally, a double lane inner row line of heavy-duty and no-load large mining haul trucks is set on the haulage bench of the end slope of the open-pit mine. For heavy-duty trucks, the truck weight plus load can reach 600 t, which brings a strong dynamic load and a great impact on the stability of the lower rib pillar (Wu et al. 2020; Chen et al. 2020). Therefore, when the highwall miner (Sasaoka et al. 2016; Bykadorov et al. 2017; Hood et al. 2020) is used to mine the retained coal in the end slope, the load on the reserved rib pillar is the composite load superimposed by the strong dynamic load generated by mining haul trucks and the non-uniform static load generated by end slopes.

Wilson's two-zone constraint theory has been widely used in the load calculation of rib pillars in strip mining in China (Wilson 1973). In this method, the ultimate bearing capacity and actual load of rib pillars are calculated, and the safety coefficient of the rib pillar is calculated to evaluate the stability of the rib pillar. Lin'Kov (2001) calculated the strength at different positions in the core area by combing the strength of the core area in a rib pillar with the actual stress, and put forward a general equation for calculating the failure envelope of the long rib pillar. Considering the strike section of the goaf as a flat elliptical orifice in an infinite plate with uniformly distributed load acting on the boundary, Pietruszczak and Mroz (1981) derived the stress calculation equation of any distance point between the rib pillar at the end of the orifice and the coal wall using the elastic fracture theory. From the 1950s to 1960s, the ideal elastic–plastic analytical solution of surrounding rock of axisymmetric circular roadway, namely Kastner equation, has been derived (Obert and Duvall 1967; Jaeger and Cook 1978). Based on the stress balance theory of loose media and the stress differential equilibrium equation, Hou and Ma (1989) obtained the coal seam interface stress and the width of the stress limit equilibrium zone (plastic zone) of the coal. Based on the elastic–plastic mechanics, Li et al. (2004) deduced the width of stress limit equilibrium zone of retained strip rib pillar using the stress balance differential equation and Coulomb criterion. Besides, the plastic zone width of the strip rib pillar was also deduced by using the complex function model of elastic theory, and then different theoretical equations were obtained.

Wang et al. (2021) explored the instability mechanism of the pillar based on cusp catastrophe theory. Aiming at the deformation characteristics of the surrounding rock, Zhang et al. (2019) proposed a calculation method based on cusp catastrophe theory. Yuan et al. (2019) investigated the surrounding rock system of the water-rich roadway in sandy-gravel stratum based on catastrophe theory, and put forward the instability criterion of the roadway in sandy-gravel stratum under seepage. Considering the cusp catastrophe theory and elastic thin plate theory, Ma et al. (2019) discussed the necessary and sufficient conditions for the instability of strip-supported rib pillars under two different constraints. Combining the Hoek–Brown failure criterion with the Mohr–Coulomb failure criterion, Wang et al. (2019) obtained the equation for calculating the width of strip-supported rib pillar, and analyzed the failure mechanism of strip-supported rib pillar and the factors affecting the width of strip-supported rib pillar in detail. Through the FLAC3D software, Porathur et al. (2014) and Wang et al. (2019) researched the slope stability under the conditions of different pillar widths and different adit widths in highwall mining of open-pit coal mines. Mo et al. (2018) quantitatively studied the influence of backfilling on rib pillar strength in highwall mining by using the FLAC2D software. Nikumbh and Chandar (2020) researched the slope stability by using the ANSYS/LS-DYNA. Woo et al. (Woo et al. 2018) studied interaction behavior between soil and a post based on Winkler’s spring theory. Sanches et al. (2020) carried out nonlinear finite element analysis on the dynamic characteristics of nonlinear Winkler foundation beam under moving load, calculated the critical velocity of moving load, and analyzed the influence of physical and geometric nonlinear characteristics of the system on the critical velocity.

Highwall mining has great practical value in terms of industrial policy and mining efficiency. However, there is little research on highwall mining in China. Due to the different application conditions of the highwall mining techniques, different coal rock occurrence and geological conditions in China and foreign countries, it is difficult to ensure the rationality and reliability of rib pillars designed by empirical equation via analyzing the stability of supporting rib pillars. In this study, the instability mechanical model of strip-supported rib pillar in highwall mining is established, which provides a basis for efficient highwall mining in an open-pit mine.

2.1 Mechanical model of the rib pillar load-bearing

The external load of the strip-supported rib pillar mainly comes from the static load of overlying strata and the dynamic load caused by the driving of mining haul trucks. Due to the action of external load, a yield zone is formed on both sides of the rib pillar. According to the two-zone constraint theory proposed by Wilson (1973), in the plastic zone (i.e. from the peak value of the rib pillar stress to the boundary of the rib pillar), the stress of the rib pillar exceeds the yield stress and flows to the goaf.

As shown in Fig. 1, the load borne by the rib pillar includes the static load of the overlying strata, and the dynamic load caused by the driving of mining haul trucks. The stress of the coal in the plastic zone of the rib pillar does not exceed the yield stress, which conforms to the elastic law. This area is surrounded by the plastic zone and constrained by the plastic zone, which is called the elastic core zone of the rib pillar. According to the effective area theory, the elastic core zone of the rib pillar carries three kinds of load, including the static load of rock strata above the upper elastic core zone of the rib pillar, the load of the overlying strata evenly divided by adjacent rib pillars and the dynamic load caused by trucks.

Load-bearing of the rib pillar in highwall mining

For the dynamic loads generated by mining haul trucks, the following four factors should be considered:

1. (1)

   The running speed of the mining haul truck on the haulage bench of an open-pit mine is generally not more than 30 km/h. Therefore, the dynamic loads caused by the horizontal movement of the mining haul truck is ignored, and only the impact load caused by the uneven road surface of the mining haul truck is considered;

2. (2)

   Because the wheel diameter of a large mining haul truck is larger than the width of the rib pillar, the front and rear wheels of the mining haul truck will not act on the same rib pillar at the same time, that is, each rib pillar can only carry the front or rear wheels of the mining haul truck;

3. (3)

   Mining haul trucks generally have 2 front wheels and 4 rear wheels. As shown in Fig. 2, under the full load, the load distribution of front and rear axles is approximately 1:2, that is, the specific pressure of wheel-to-ground is approximately equal and the wheel-to-ground contact area is approximately the same;

4. (4)

   When the rear wheel of a fully-loaded mining haul truck acts directly above the rib pillar, the rib pillar bears the maximum external load. At this time, the wheel-to-ground contact area of the mining haul truck is calculated by a × b, where a is the length of the wheel-to-ground contact area and b is the wheel width. Combined with factor (1), the wheel-to-ground load can be simplified as a uniformly distributed static load, with a width of a, the load value of q₀ and the dynamic load coefficient of k_d.

Calculation diagram of wheel-to-ground contact area during mining haul truck driving

When the cutting width of the highwall miner is 2.9–3.7 m (Chen et al. 2013; Yin and Dong 2015; Kuznetsova and Anfyorov 2019), no caving occurs in the roof of the rib pillar roof. As shown in Fig. 1, it is assumed that the average thickness of the overlying strata is H and the unit weight is γ, the width of the rib pillar is W_p and the width of the mining opening is W_m, then the load borne by the rib pillar (in plane) is:

The relationship of the rib pillar stress σ, strain ε and damage parameter D can be expressed as follows (Xie et al. 2009; Xie 1990):

where, ε₀ is the strain of the rib pillar under load; E is the initial elastic modulus of the rib pillar.

Assuming that the width of the plastic zone on one side of the rib pillar is Y and the thickness of the coal seam is H_p, then in the plastic zone of the rib pillar with a total width of 2Y or in the elastic core zone of the rib pillar with a width of W_p − 2Y, the relationship between plastic zone load P_s, elastic core zone load P_e and deformation u can be expressed as follows:

where u₀ is the compression of the rib pillar under load.

2.2 The rib pillar failure and instability model based on the cusp catastrophe theory

The strain energy V₁ of the plastic zone of the strip-supported rib pillar, the elastic potential energy V₂ of the elastic core zone of the strip-supported rib pillar and the gravity potential energy V₃ of the overlying strata and the mining haul truck can be obtained as follows:

The total potential energy of the roof-rib pillar system composed of rib pillar, overlying strata and mining haul truck can be expressed as:

Taking u as the state variable, the catastrophe theory is used for analysis. The first derivative of V(u) is taken and V(u)' = 0, then the equation of the equilibrium surface can be obtained:

Equation (8) is derived further as:

After solving Eq. (9), then:

The equilibrium surface equation V(u)' is expanded into a power series according to Taylor equation at u = u₁ = 2u₀, and is taken to the cubic term. Then it can be simplified as:

The dimensionless quantity ζ is introduced as the state variable, p and q as the control variables, and let:

where, k₀ is the stiffness ratio, k₀ = k_e/k_s; k_e is the stiffness of medium in the elastic core zone of the rib pillar, k_e = E(W_p − 2Y)/H_p; k_s is the stiffness of the medium in the plastic zone of rib pillar, k_s = 2EYe⁻²/H_p; t is defined as a parameter related to mining conditions, i.e. mining depth H, mining opening width W_m, rib pillar width W_p and unit weight γ, the coal seam thickness H_p and relevant mechanical parameters of coal samples. Then Eq. (13) can be obtained as follows:

By substituting Eq. (12) and Eq. (13) into Eq. (11), an equilibrium equation (with ζ as a state variable, p and q as control variables) can be obtained, which conforms to the standard form of cusp catastrophe:

By solving the simultaneous Eq. (14) and inflection point equation 3ζ² + p = 0, the bifurcation set equation can be obtained:

By substituting Eq. (12) and Eq. (13) into the bifurcation set Eq. (15), the bifurcation set Eq. (15) of the cusp catastrophe (Chen et al. 2021; Dong et al. 2021) of the system can be obtained:

2.3 Criterion equation of the rib pillar stability

The stress of the rib pillar σ_s can be obtained as:

According to the ultimate strength theory, the safety factor f of the rib pillar can be obtained from the ultimate stress σ_p of the rib pillar and stress σ_s borne by the rib pillar.

Due to the difference in actual mining conditions and geological conditions of the highwall mining, the selection of the safety factor of the rib pillar is different. According to the previous studies of highwall mining, the safety factor of the rib pillar is generally greater than 1.3 (Chen et al. 2013; Zipf and Mark 2005).

As described in Sect. 2.2, to maintain the stability of the rib pillar in highwall mining, the following two conditions must be met simultaneously for the stability of the rib pillar under highwall mining:

According to Eq. (19), there are three main variables in the mechanical criterion of the rib pillar instability, including the plastic zone width Y of the rib pillar, the ultimate stress σ_p borne by rib pillar, and the compression u₀ of the rib pillar under external load.

3.1 Calculation of plastic zone width and ultimate stress of the rib pillar

Based on the limit equilibrium theory, it is assumed that the roof and floor of coal seams are consistent, and the compressive strength is greater than that of coal. Without considering the physical strength, the mechanical analysis model of the lower mining opening and surrounding rock in highwall mining is established (Gu et al. 2014). As shown in Fig. 3, in highwall mining, the mining opening width is W_m, the rib pillar width is W_p, the mining opening height is H_p, and the coal seam buried depth is H, the vertical stress of σ₀ is γH + [k_dq₀a/(W_m + W_p)], the horizontal stress at infinity from the mining opening is λσ₀, side pressure coefficient λ = μ/(1 − μ), and μ is the Poisson's ratio of coal.

Mechanical model of mining opening and surrounding rock

The model is bilaterally symmetrical along the central line of the mining opening and longitudinally symmetrical along the x-axis. A micro element in the plastic zone of the rib pillar is taken, as shown in Fig. 4. Because the plastic zone is the limit equilibrium zone, the average shear stress at the interface between roof and floor and coal seam τ and vertical load σ_y meet the function:

where c₀ and φ₀ are the cohesion and internal friction angle of the interface between coal seams and the roof and floor.

Micro element mechanical model of the plastic zone

According to the Mohr–Coulomb yield criterion, the plastic zone of coal is in the limit equilibrium state, and the following equation can be obtained:

The vertical stress and horizontal stress of micro elements are the main stress. In Eq. (21), c and φ are cohesion and internal friction angle of coal, and σ₁ = σ_y, σ₃ = σ_x.

In Eq. (21), let A = (1 + sin φ)/(1 − sin φ), then boundary conditions \(\left. {\sigma_{x} } \right|_{x = 0} = 0\) are substituted. Given that N₁ = (2ccos φ)/(1 + sin φ), if x \(\in\)[0,Y], the vertical stress \(\sigma_{y}^{{\text{p}}}\) and horizontal stress \(\sigma_{x}^{{\text{p}}}\) of micro elements in the plastic zone can be obtained:

The vertical ultimate stress σ_p of the rib pillar can be obtained from the Eq. (22):

Combined with the principle of stress balance, after the mining opening is excavated, the load originally borne by the coal in the mining opening section is transferred to the rib pillars on both sides of the mining opening. Assuming that the vertical stress and horizontal stress of micro elements are equal to the average stress distributed along with the coal seam thickness, and let X = (2YAtan φ₀)/H_p, then Eq. (24) can be obtained:

The expression of plastic zone width Y of the rib pillar on both sides of the mining opening is finally obtained by solving the equation:

3.2 Calculation of the rib pillar compression under the dynamic load of mining haul trucks

3.2.1 Theoretical model of the beam on Winkler elastic foundation

To establish the differential equation of the deflection curve satisfied by the beam (in Fig. 5), a section of the beam with a length of l, a width of b₁ and a height of h is taken, and the foundation coefficient of the elastic foundation is presented as k and the width as b₂. Under the load q(x), the compression of the foundation is y(x), the reaction of the foundation to the beam is q_f(x), a micro section dx is intercepted in the beam, the bending moment is M and the shear force is Q.

Elastic foundation beam model composed of the rib pillar and the roof

According to Winkler's assumption (Soon et al. 2021; Li et al. 2021), the base reaction force is:

If the influence of shear force on beam deflection is not considered, according to the balance condition of the force, α = (kb₂/4E₁I)^1/4 and αx can be used to replace x, then:

The additional term caused by uniformly distributed load q(x) = q₀ is introduced, Eq. (28) can be obtained:

where

3.2.2 Rib pillar compression under dynamic loads of the mining haul truck

Since the compression of the rib pillar under the dynamic loads of a fully-load truck is greater than that of a no-load truck, only the compression of the rib pillar under the dynamic load of a fully-loaded truck is discussed. As shown in Fig. 6, the rib pillar group between two permanent rib pillars is taken as the research object along the direction of the haulage bench in the end slope.

Dynamic load of front and rear wheels of a mining haul truck under full load

When calculating the compression of dynamic loads of a mining haul truck on the rib pillar, the rib pillar-roof model can be regarded as a plane strain model. In this model, displacement boundary conditions are adopted, and fixed boundary conditions are employed by the left and right sides. Combined with the Winkler elastic foundation beam theory, the coal seam after highwall mining is simplified into countless parallel springs, and the rock stratum above the rib pillar is simplified into a beam. Therefore, this problem can be simplified to the solution of the compression of the rib pillar in the middle of the rib pillar group (near y = 0) under the dynamic loads of the truck.

As shown in Fig. 7, the coordinate system yoz is established, the y-axis is established along the horizontal direction, and the z-axis is established along the vertical direction. The boundary conditions of the calculation model are:

Schematic diagram of the coordinate system calculated by the model

The dynamic loads of a mining haul truck are further simplified. Assuming that the compression generated on the rib pillar in the middle of the rib pillar group (near y = 0) does not decay when the truck is driving, the compression generated by dynamic loads of a truck on a certain rib pillar can be simplified as the sum of the compression generated by the static load between the truck and the rib pillar at different distances, and this compression is defined as u₂.

As shown in Figs. 6 and 7, the width of the mining opening is W_m, the height is H_p, the width of the rib pillar is W_p, the elastic modulus is E, the width of the plastic zone is Y, the width of the elastic foundation is unit one, and the foundation coefficient is k = [E(W_p − 2Y)]/[H_p(W_m + W_p)]; the height of the roof is h, the length is 2l₁, the width of the roof is unit one, the elastic modulus is E₁, and the elastic index value is α = [3kb₂/(E₁b₁h³)]^1/4. The interval of the uniformly distributed load is y \(\in\)[y₁,y₂], y \(\in\)[y₃,y₄]. By substituting the boundary conditions (30) into Eq. (28), the initial parameters z₀, θ₀, M₀ and Q₀ can be obtained. Subsequently, the calculation equation of coal seam compression z(y) along the y-axis can be obtained.

The obtained z(y) curve is integrated on the interval of y \(\in\)[− l₁,l₁]. Finally, the compression u₂ of the rib pillar caused by the dynamic loads of the fully-loaded truck is obtained:

3.3 Calculation of the rib pillar compression under the combined static load of overlying strata and mining haul truck

The loads of overlying strata that deform the elastic foundation are q₁, q₂, …, q_n(n = \(\in\) N⁺), corresponding to the coordinate intervals x \(\in\)[x₁,x₂], x \(\in\)[x₂,x₃], …, x \(\in\)[x_n,x_n+1]. The loads of the mining haul truck that deform the elastic foundation are q_n+1, q_n+2, …, q_n+m, (m = \(\in\) N⁺ and m is a multiple of 4), corresponding to the coordinate interval x \(\in\)[x_n+1,x_n+2], x \(\in\)[x_n+2,x_n+3], …, x \(\in\)[x_2m+n,x_2m+n+1], where q_n+1, q_n+3, q_n+5, …, q_n+m-3, q_n+m-1 are the static loads of the rear wheels of the truck under no-load; q_n+2, q_n+4, q_n+6, …, q_n+m-2, q_n+m are the static load of the rear wheel of the truck when it is fully loaded. Figure 8 shows the established coordinate system xoz.

Calculation model of the rib pillar compression under the static load of overlying strata and mining haul truck

As shown in Figs. 6 and 8, the width of the roof is b₁ = W_m + W_p, and the width of the elastic core area of the rib pillar is b₂ = W_p − 2Y. The foundation coefficient is k = E/H_p, and the elastic index value of roof (beam) is α = {3 k(W_p − 2Y)/[E₁h³(W_m + W_p)]}^1/4.

As shown in Fig. 8, the boundary conditions of the model are:

By substituting the boundary conditions into Eq. (28), the initial parameter z₀ and θ₀ can be obtained by solving the equations. Thus, the coal seam compression z(x) along the x-direction of the coordinate axis can be obtained:

When i = 1, 2, 3,…, n, x_i and x_i+1 are the left and right endpoints of the interval where q_i is located; When i = 1, 2, 3,…, m, x_2i+n+1 and x_2i+n are the left and right endpoints of the interval where q_i+n is located.

By analyzing the z(x) curve, the maximum compression u₃ and corresponding coordinate x of the rib pillar under the static load of overlying strata and mining haul truck can be obtained:

To sum up, the total compression u₀ of the rib pillar under dynamic and static loads superposition is equal to the sum of the rib pillar compression u₂ generated by dynamic loads of a truck and the rib pillar compression u₃ generated by the static load of overlying strata and truck, that is:

Figure 9 shows the schematic diagram of the end slope of an open-pit mine in highwall mining. The average thickness of the coal seam closest to the surface is 3 m and the buried depth is 122 m. After open-pit mining, the width of the compressed coal seam in the end slope is about 190 m, the width of each haulage bench is about 40 m and that of each safety bench is about 10 m, the slope angle of the loess layer is 55° and that of each rock layer is 70°. The physical and mechanical parameters of coal and rock strata of end slope are shown in Table 1.

Section view of end slope in an open-pit mine

When the mining haul truck is fully loaded (600 t), according to the calculation of design parameters, the bearing capacity of each wheel is 100 t and wheel-to-ground contact area is 0.440 m × 1.250 m. When the mining haul truck is unloaded (237 t), each wheel of the front wheel carries 54.5 t, each wheel of the rear wheel carries 32 t and the wheel-to-ground contact area is 0.357 m × 1.250 m.

When the highwall miner is used to recover the retained coal in the end slope of the open-pit mine, the mining height and mining width are 3 m. The width of the rib pillar is designed by using the stability model of the lower rib pillar of highwall mining established in Sec. 2.

According to the established rib pillar instability model based on the cusp catastrophe theory, the plastic zone width of the rib pillar, the ultimate stress of the rib pillar and the compression of the rib pillar are calculated respectively. The reasonable rib pillar width W_p is determined by analyzing the criterion values Δ of catastrophe instability of rib pillars with different widths and the safety factor value f of the rib pillar.

4.1 Calculation of plastic zone width of the rib pillar

The plastic zone width of the rib pillar includes the sum of the plastic zone width caused by static load (slope) and the plastic zone width caused by dynamic load (truck). The rib pillar under the No. 1 haulage bench is taken as an example. According to Table 1, the cohesion c of the coal and the cohesion c₀ of the interface between the coal seam and the roof and floor are set as 1.62 MPa; the friction angle φ in the coal and the internal friction angle φ₀ at the interface with roof and floor are set as 36°; the mining opening width is W_m = 3 m, and the mining opening height is H_p = 3 m; the dynamic load coefficient is k_d = 1.2, and the static load is q₀ = 1.8 MN/m when a mining haul truck is fully loaded; the width of the uniformly distributed static load is a = 0.44 m; the vertical stress for No. 1 haulage bench is γH = 2491.6 kPa. The plastic zone width Y of the rib pillar with different pillar widths W_p beneath the No. 1 haulage bench can be obtained by solving Eq. (24) and Eq. (25), as shown in Fig. 10.

Plastic zone width of the rib pillar with different pillar widths beneath the No. 1 haulage bench

As shown in Fig. 10, for the rib pillar directly below the No. 1 haulage bench, the plastic zone width Y on both sides of the rib pillar gradually decreases with the increase of the width W_p of the rib pillar. When the pillar width W_p increases from 0.54 to 3.00 m, the pillar plastic zone width Y decreases from 0.269 to 0.260 m, decreasing by 3.3%. The plastic zone width Y of the rib pillar decreases linearly.

4.2 Calculation of the rib pillar compression

4.2.1 Compression of the rib pillar under dynamic loads of the mining haul truck

Along the direction of a mining opening, the rib pillar is equivalent to the elastic foundation. As the result of the highwall mining, the load above the elastic foundation is redistributed. When the fully-loaded truck and unloaded truck meet on each two-lane haulage bench at the same time, then the maximum static load is generated by trucks.

As shown in Fig. 7, one barrier pillar is reserved in every ten rib pillars, and the width of the barrier pillar is 12 m; then half the length of the model l₁ = (22.5 + 5W_p) m is taken, the uniformly distributed load q₀ = 1.8 MPa, the width of the uniformly distributed load is a = 0.44 m, and the abscissa of the uniformly distributed load is y₁ = −3.61 m, y₂ = −3.17 m, y₃ = 3.17 m and y₄ = 3.61 m.

For the No. 1 haulage bench, the roof height is h = 107 m, the width is unit one, and the elastic modulus is E₁ = 3.03 GPa. The elastic modulus of the rib pillar is E = 1 GPa, the width of the elastic foundation is unit one, and the foundation coefficient is k = [E(W_p − 2Y)]/[H_p(W_m + W_p)].

The rib pillar plastic zone width Y obtained with different rib pillar reserved width W_p (W_p \(\in\)[0.54,3.00] m) is substituted into the foundation coefficient k; a series of obtained foundation coefficients k is substituted into the Eq. (28), so as to obtain the coal seam compression z(y) along the coordinate axis y under different rib pillar reserved width W_p. The obtained curve z(y) is integrated on the interval y \(\in\)[− l₁,l₁], and then the compression u₂ of the rib pillar due to the dynamic load generated by the fully-loaded truck is obtained, as shown in Fig. 11.

Compression of the rib pillar under the dynamic load of truck at full load under the No. 1 haulage bench

4.2.2 Compression of the rib pillar under the combined action of static loads of overlying strata and mining haul truck

The loads of overlying strata that deform the elastic foundation are q₁, q₂, …, q₇, and the loads of the truck that deform the elastic foundation are q₈, q₉, …, q₂₃, therefore, n = 7, m = 16 in the Eq. (33). The overlying strata load at all levels can be expressed as follows:

The static load of the unloaded truck is q₈ = … = q₂₂ = 307.2 kN/m, corresponding to interval x \(\in\)[5.675,8.175] m,…, x \(\in\)[161.825,164.325] m.

The static load of a fully-loaded truck is q₉ = … = q₂₃ = 960.0 kN/m, corresponding to interval x \(\in\)[25.675,28.175] m,…, x \(\in\)[181.825,184.325] m.

The parameter setting is as follows: the roof length is l₂ = 500 m, the roof height is h = 12 m, the roof elastic modulus is E₁ = 4 GPa, the mining opening height is H_p = 3 m, the mining width is W_m = 3 m, and the foundation coefficient is k = 1/3 GPa. As shown in Fig. 8, the boundary conditions of the mechanical model are shown in Eq. (32). By substituting the boundary conditions (32) into the Eq. (33), the coal seam compression z(x) along the x-direction, the maximum compression u₃ of the rib pillar under the static load of overlying strata and truck and the corresponding coordinate x is obtained.

As shown in Fig. 12, with the increase of the rib pillar width, the compression of the rib pillar under static load approximately decreases nonlinearly. When the rib pillar width W_p increases from 0.54 to 0.80 m, the position of the maximum compression u₃ of the rib pillar decreases sharply from x = 12,942.73 mm to x = 99.20 mm. When the rib pillar width W_p continues to increase to 3.00 m, the position of u₃ decreases gently to x = 11.16 mm.

The maximum compression of the rib pillar under the static load of overlying strata and mining haul truck

As shown in Fig. 13, with the increase of the rib pillar width W_p, the position of the maximum compression u₃ of the rib pillar gradually shifts to the position of maximum buried depth of the rib pillar. When the rib pillar width W_p increases from 0.54 m to 0.70 m, the position of maximum compression u₃of the rib pillar shifts sharply from x = 129.00 m (below the No. 2 haulage bench) to x = 157.25 m (below the No. 1 haulage bench); When the rib pillar width W_p continues to increase to 3.00 m, the position of maximum compression u₃ of the rib pillar slowly shifts from x = 157.25 m to x = 162.75 m.

Position of the maximum compression of the rib pillar

4.2.3 Compression of the rib pillar under superposition of dynamic and static loads

Through the analysis of the compression curve z(x) of the rib pillar below the No. 1 haulage bench (x \(\in\)[150,190] m) under static load, it can be obtained that: when the rib pillar width W_p \(\in\)[0.54,0.58] m, the curve z(x) decreases monotonically in the interval, then the maximum compression of the rib pillar below the No. 1 haulage bench under static load occurs at the position of x = 150 m; when the rib pillar width W_p \(\in\)[0.54,3.00] m, the curve z(x) first increases in the interval, reaches the peak and then decreases.

As mentioned above, the total compression is u₀ = u₂ + u₃. Moreover, the maximum compression of the rib pillar generated by dynamic loads of the truck is located directly below the contact between the wheel and the ground, that is, for the No. 1 haulage bench, the rib pillar compression u₂ generated by dynamic loads of the truck has a superposition effect with the rib pillar compression u₃ generated by static loads of overlying strata and the truck only at x \(\in\)[155.675,158.175] m, x \(\in\)[161.825,164.325] m, x \(\in\)[175.675,178.175] m, and x \(\in\)[181.825,184.325] m The results are shown in Fig. 14.

Compression of the rib pillar under dynamic and static load superposition of the No. 1 haulage bench

From Fig. 14, it is concluded that with the increase of the rib pillar width, the total compression u₀ of the rib pillar approximately decreases nonlinearly. When the rib pillar width W_p increases from 0.54 m to 0.80 m, the total compression u₀ under the No. 1 haulage bench decreases sharply from 12,471.42 mm to 99.20 mm. When the rib pillar width W_p continues to increase to 3.00 m, the total compression u₀ decreases gently to 11.16 mm.

The curve shown in Fig. 14 is consistent with that in Fig. 12. For the rib pillar directly below the No. 1 haulage bench, the rib pillar compression u₂ generated by dynamic loads of the truck is much smaller than the compression u₃ generated by static loads of overlying strata and the truck. The rib pillar compression u₃ generated by static loads of overlying strata and the truck plays a decisive role in the total rib pillar compression u₀ below the No.1 haulage bench.

4.3 Verification of the rib pillar stability criterion equations

4.3.1 Criterion value of the rib pillar stability

When the rib pillar width W_p takes different values (W_p \(\in\)[0.54,3.00] m), the calculated plastic zone width Y of the rib pillar and the total compression u₀ of the rib pillar are substituted into Eq. (16), then criterion values Δ under different rib pillar widths can be obtained. The calculation results are shown in Fig. 15.

Mechanical criterion value Δ of the rib pillar instability under dynamic and static loads of the No. 1 haulage bench

As shown in Fig. 15, the criterion value Δ sharply increases with the gradual increase of the rib pillar width W_p. When the rib pillar width W_p increases from 0.54 m to 3.00 m, the criterion value Δ increases from 7 to 80,065 in an approximate power function.

4.3.2 Ultimate stress and safety factor of the rib pillar

When the rib pillar width takes different values, the calculated plastic zone width Y of the rib pillar and the total compression u₀ of the rib pillar can be substituted into Eqs. (17) and (18). The ultimate stress σ_p and safety factor f under different rib pillar widths can be obtained, and the calculation results are shown in Figs. 16 and 17.

Ultimate stress of the rib pillar under the No. 1 haulage bench

Safety factor of the rib pillar under the No. 1 haulage bench

As shown in Figs. 16 and 17, with the increase of the rib pillar width, the ultimate stress of the rib pillar under the No. 1 haulage bench approximately decreases linearly, and the safety factor of the rib pillar approximately increases linearly. When the rib pillar width W_p increases from 0.54 to 3.00 m, the ultimate stress of the rib pillar gradually decreases from 11.965 to 11.733 MPa, and the safety factor of the rib pillar increases from 0.66 to 2.21. The ultimate stress of the rib pillar decreases slightly, with a reduced rate of 1.9%, indicating that the ultimate stress is less affected by the change of the rib pillar width. The safety factor of the rib pillar increases greatly, with a growth rate of 234.8%, which is greatly affected by the change of the rib pillar width.

4.4 Analysis on reasonable reserved width of the rib pillar

Combined with Figs. 15 and 17, the criterion value Δ of the rib pillar under the No. 1 haulage bench and statistical chart of the safety factor f can be obtained. According to the crust catastrophe instability model of rib pillars in highwall mining, if the rib pillar remains stable, criteria value Δ > 0 is required. When the safety factor of the rib pillar f > 1.3 is taken, Eq. (37) can be obtained by solving the simultaneous equations of Eq. (19), Eq. (23), Eq. (24) and Eq. (25):

Based on the above calculation, the rib pillar width under the No. 1 haulage bench is at least 1.3 m.

Based on the calculation of the rib pillar width under the No. 1 haulage bench, when the width of rib pillar is W_p = 1.3 m, the safety factor f and criterion value Δ of the No. 2, 3 and 4 haulage benches are checked, and the results are shown in Table 2.

Based on the calculation of the No. 1–4 haulage benches, the width of rib pillar under this mining condition is at least 1.3 m.

4.5 Verification of numerical simulation by the ANSYS/LS-DYNA

Referring to the dimensions and physical and mechanical parameters of the end slope model (in Fig. 9 and Table 1), the end slope model in Fig. 18 is established in the ANSYS/LS-DYNA. A permanent barrier pillar is reserved in every ten rib pillars. The mining opening is 3 m high, 3 m wide and 190 m long. The rib pillar is 190 m long, 1.3 m wide and 3 m high. The cross-section is stretched 58 m forward along the z-axis, and 11 cuboids as mining openings are reserved.

Global view of end slope

To eliminate the influence of boundary conditions and reduce the calculation amount of the model, two rib pillars in the middle of the rib pillar group are taken as the main research objects and the mesh refinement is carried out. After the initial in-situ stress is balanced, the life-and-death unit command is used to delete the grid of the mining opening (Fig. 19).

Meshing after highwall mining

As shown in Fig. 20, a moving load on each haulage bench is applied, and the load of a mining haul truck is simplified into four-point loads. The load parameters are set referring to Chapter 4, and the load moving speed is 30 km/h.

Schematic diagram of mesh division and dynamic load application of haulage bench

After entering the ANSYS post-processor, the calculation results of the last step are read, and the unit after end slope mining is selected. As shown in Fig. 21, the end slope does not produce an obvious plastic zone, so the slope is in a stable state under this mining condition. Due to the symmetry of the model, the unit belonging to the rib pillar group in the most middle part is further selected. As shown in Fig. 22, the rib pillar does not produce an obvious plastic zone at the shallow buried depth, but an obvious plastic zone can be seen at the deeply buried position (below the No. 1 haulage bench). The section diagram at the plastic strain peak of the rib pillar node is obtained by enlarging the section in the middle red frame area of Fig. 22. The plastic area of the rib pillar is not connected, and the width of the plastic area is about 0.26 m. The results are in good agreement with the theoretical calculation, which verifies the feasibility of the proposed theoretical calculation.

Cloud diagram of the plastic strain of the slope global nodes

Cloud diagram of plastic strain at the node of the middle rib pillar of rib pillar group

Based on the catastrophe theory, energy theory and rock mass mechanics, the mechanical analysis and numerical simulation are performed, the instability mechanism of the rib pillar under dynamic and static loads in a highwall mining system is comprehensively studied, and the design method of the rib pillar width in highwall mining is put forward. The main research results are as follows:

1. (1)

   According to the stress characteristics of the rib pillar in highwall mining, the cusp catastrophe theory and the requirements of rib pillar safety coefficient are considered, a comprehensive criterion for predicting the rib pillar stability is established. The comprehensive criterion is mainly determined by three factors: the plastic zone width of the rib pillar, the ultimate stress of the rib pillar and the corresponding compression of the rib pillar under external load.

2. (2)

   According to the limit equilibrium theory, the limit stress of the rib pillar is analyzed, and the calculation equation of the plastic zone width of the rib pillar on both sides of the mining opening is obtained. The results show that the plastic zone width of the rib pillar can be determined by three factors: the mining height, the ratio of vertical stress peak to horizontal stress at the elastic–plastic boundary, and the internal friction angle between coal seam and roof and floor interface. Based on the Winkler foundation beam theory, an elastic foundation beam model composed of the rib pillar and roof under highwall mining is established, and the calculation equation of compression of the rib pillar under load is obtained.

3. (3)

   With the increase of the rib pillar width, the plastic zone width of the rib pillar decreases linearly; under the action of dynamic and static loads, the total compression of the rib pillar decreases nonlinearly. In the total compression of the rib pillar, the compression of the rib pillar caused by static loads of overlying strata and the mining haul truck plays a decisive role.

4. (4)

   Taking the actual mining conditions of the open-pit mine as the engineering background, the 3 m thick coal seam with a buried depth of 122 m is mined by the highwall miner through the numerical simulation. According to the established rib pillar instability model of the highwall mining system, a theoretical calculation is performed. The results show that when the mining opening width is 3 m, the reasonable width of the rib pillar is at least 1.3 m, and the safety factor of the rib pillar is 1.3. Through the verification of the numerical simulation test via the ANSYS/LS-DYNA, the theoretical calculation results are in good agreement with the numerical simulation results.