Controlled Roof Collapse during Secondary Mining in Coal Mines

Abstract

The problem considered is an investigation of the possible collapse of the roof between the pillar next to be mined in secondary coal mining and the first line of pillar remnants called snooks. The roof rock between the pillar, which is the working face, and the snook is modelled as an Euler-Bernoulli beam acted on at each end by a horizontal force and by its weight per unit length. The beam is clamped at the pillar and simply supported (hinged) at the snook. The dimensionless differential equation for the beam and the boundary conditions depend on one dimensionless number 𝐵. We consider the range of values of 𝐵 before the displacement and curvature first become singular at 𝐵=𝐵₁. The model predicts that for all practical purposes, the beam will break at the clamped end at the pillar. The failure of the beam for values of 𝐵 greater than 𝐵₁ is investigated computationally.

1. Introduction

We consider the challenge posed by coal mine pillar extraction [1, 2]. Secondary mining involves revisiting a mine and extracting coal from the pillars. The mining of these pillars commences from the area furthest away from the point of entry of the mine. This exercise involves cutting the existing pillars into smaller pillars called snooks. As each section is mined, the roof must collapse in a controlled manner in order to pose no safety risk to those miners operating underground. We analyse the behaviour of the roof of the mine between the pillar next to be mined and the first line of snooks. This is the work area and must be safe for the miners.

2. Model

In Figure 1, a mining panel in shown prior to pillar extraction. The tunnels are excavated in coal which are approximately 5 m to 7 m wide. They are excavated in a fixed pattern crossing at right angles creating a checker board layout. The coal between the tunnels forms the pillars which support the overburden rock. The width of the pillars is approximately 10 m to 20 m wide and is a function of the depth of the mine. The height of the tunnels ranges from 3 m to 4 m. Secondary mining is carried out in two stages. In the initial stage, approximately 5 to 10 pillars are removed and the roof is left to collapse. This stage is modelled in [2]. Following this, adjacent pillars are mined and smaller sections are left to collapse. The purpose of this paper is to model the second stage in the extraction process. Figure 2 shows the snooks after pillar extraction. The pillars are cut to leave four snooks, approximately 2 m, one at each corner. The snooks have to be small enough to fail when the miners are a safe distance (about the width of a pillar) from the working face but they have to be large enough to be stable right next to the unmined pillars.

Figure 1 

A mining panel showing the pillars before pillar extraction. Reproduced with permission from N. van der Merwe.

Figure 2 

A mining panel showing the snooks after pillar extraction. The pillar is replaced by four snooks, one at each corner. Reproduced with permission from N. van der Merwe.

The roof consists of horizontal layers of rock of approximate thickness 0.5 m to 20 m, as shown in Figure 3. The Euler-Bernoulli beam equation can be used to describe the horizontal layers of rock in the roof. The use of the Euler-Bernoulli beam equation assumes that the roof is thin compared with its horizontal extent and that only the horizontal direction is important. The horizontal extent of the beam is the distance from the next pillar to be mined to the first line of snooks which is the width of the tunnel and is approximately 6 m. The ratio of the thickness of the beam to its length ranges from about 0.1 to about 3 and thus for the theory to apply the thickness of the beam should not exceed 2 m. The width of the mining panel ranges from about 130 m to 200 m. If we take the width of the mining panel as the width of the beam then the ratio of the length of the beam to its width varies from about 0.05 to 0.03. Dependence of the variables in the direction of the width of the beam can therefore be neglected. The use of the Euler-Bernoulli beam is therefore justified for a beam of thickness less than about 2 m.

Figure 3 

A cross-section of the mine showing failed snooks and goaf, a stable snook, the next pillar to be mined, and the overburden which consists of horizontal layers of rock. Reproduced with permission from N. van der Merwe.

In this paper we will investigate if roof collapse can occur between the next pillar to be mined and the first line of snooks when these snooks are stable and do not fail. In order to achieve this, we consider the roof of the mine to be clamped at a pillar while at the adjacent snook, the roof is simply supported or hinged. We consider the roof to be simply supported at the snook since as secondary mining takes place, disturbances in the rock mass and also the roof collapse due to the failure of the neighbouring snooks could change the roof structure in the region where the snooks support the beam [3]. We model this by assuming that at the snook the beam is no longer clamped and use instead that the beam is simply supported or hinged at the snook. We also consider the behaviour of this small section of the roof when a disturbance, such as a seismic event, causes a sudden increase in the horizontal force acting at each end of the section of the roof.

An analysis of the problem where both ends of the beam are clamped is presented in [2]. The beam number 𝐵, which occurred in the dimensionless Euler-Bernoulli beam equation, was defined in [2] as follows: 𝑃𝐵=𝐿𝐸𝐼^1/2, (2.1) where 𝑃 is the horizontal axial force applied to the ends of the beam, 𝐿 is the length, 𝐸 is the Young’s modulus of the beam, and 𝐼 is a second moment of the cross-sectional area of the beam. The beam number was the only dimensionless parameter in the problem. This number has arisen before in the literature, for example in [4], but no name was assigned to it. The displacement became singular when 𝐵=2𝑛𝜋 where 𝑛=1,2,3…. The magnitude of the displacement was greatest at the centre of the beam for 0<𝐵<2𝜋 due to the symmetry of the problem. For 0<𝐵<2𝜋 the magnitude of the curvature was greatest at the endpoints of the beam and thus the beam collapsed at these points when the tensile strength of the beam was exceeded. The problem of one end clamped and one simply supported is not symmetric. Our task is to solve and analyse this problem and to compare it with the problem with both ends clamped.

Previous work on roof failure due to the failure of snooks is reviewed in [5]. Useful texts are [6–9].

3. Derivation of the Differential Equation

The combined beam and strut is shown in Figure 4. We will use the notation and conventions of Segal and Handelman [10]. The coordinate axes are defined in terms of the undeformed beam. The 𝑥₁- and 𝑥₂-axes are along the axes of principle moments of inertia of the cross-section of the beam with the 𝑥₁-axis vertically downwards. The 𝑥₃-axis is horizontal and passes through the centroid of each cross-section. The origin of the coordinate system is at the centroid of the cross-section of the left end of the beam. Unit vectors, i, j, and k are directed along each coordinate axis. For simplicity, we denote 𝑥₃ by 𝑥.

Figure 4 

Beam with end at 𝑥=0 clamped and end 𝑥=𝐿 simply supported (hinged).

An outline of the derivation of the differential equation for the Euler-Bernoulli beam when both ends are clamped is given in [2, 10]. The potential energy 𝑉 of an elastic beam of length 𝐿 and Young’s modulus 𝐸 is given by the following: 𝑉(𝑤)=^𝐿₀12𝑑𝐸𝐼²𝑤𝑑𝑥²²−12𝑃𝑑𝑤𝑑𝑥²−𝑞(𝑥)+𝑠(𝑥)𝑤(𝑥)𝑑𝑥, (3.1) where 𝑤(𝑥) is the displacement of the beam from the horizontal position 𝑥 vertically downwards in the direction of i, 𝐼 is the second moment of area about the 𝑥₂-axis, 𝑞(𝑥) is the magnitude of the body force per unit length, 𝑠(𝑥) is the magnitude of the applied surface traction per unit length in the direction of 𝑖, and 𝑃 is the horizontal force acting at each end of the beam. We assume that the simple support or hinge can oppose an axial force thus disallowing any axial motion. The nonlinear strain tensor is used in part of the derivation of 𝑉(𝑤) in [10].

The derivation of the beam equation depends on the boundary conditions. We show that the boundary conditions for the present problem yield the same beam equation as in [2] and thus the only difference between the two problems is their boundary conditions. The boundary conditions with one end clamped and the other end simply supported are as follows: 𝑤(0)=0,𝑑𝑤(0)(𝐿)𝑑𝑑𝑥=0,𝑤=0,²𝑤𝑑𝑥²(𝐿)=0. (3.2)

At equilibrium, the potential energy is at an extremum. In order to impose this condition we let 𝑤(𝑥)=𝑤₀(𝑥)+𝜖𝑤₁(𝑥), (3.3) where 𝜖 is a constant parameter. Since 𝑤(𝑥) satisfies boundary conditions (3.2) for all 𝜖 it implies that 𝑤₀(𝑥) and 𝑤₁(𝑥) must separately satisfy (3.2). Thus we have, 𝑑𝑉𝑤𝑑𝜖⁰+𝜖𝑤₁|_𝜖=0=^𝐿₀𝑑𝐸𝐼²𝑤₀𝑑𝑥²𝑑²𝑤₁𝑑𝑥²−𝑃𝑑𝑤₀𝑑𝑥𝑑𝑤₁−𝑞(𝑥)(𝑥)𝑤𝑑𝑥+𝑠₁(𝑥)𝑑𝑥, (3.4) giving that, at equilibrium, 0=^𝐿₀𝑑𝐸𝐼²𝑤₀𝑑𝑥²𝑑²𝑤₁𝑑𝑥²−𝑃𝑑𝑤₀𝑑𝑥𝑑𝑤₁−𝑞(𝑥)(𝑥)𝑤𝑑𝑥+𝑠₁(𝑥)𝑑𝑥. (3.5)

We use integration by parts and the boundary conditions to deduce the following results: ^𝐿₀𝑑²𝑤₀𝑑𝑥²𝑑²𝑤₁𝑑𝑥²𝑑𝑥=^𝐿₀𝑑⁴𝑤₀𝑑𝑥⁴𝑤₁(𝑥)𝑑𝑥,_𝐿0𝑑𝑤₀𝑑𝑥𝑑𝑤₁𝑑𝑥𝑑𝑥=−^𝐿₀𝑑²𝑤₀𝑑𝑥²𝑤₁(𝑥)𝑑𝑥. (3.6) Using (3.6) we can rewrite (3.5) as 0=^𝐿₀𝑑𝐸𝐼⁴𝑤₀𝑑𝑥⁴𝑑+𝑃²𝑤₀𝑑𝑥²−𝑞(𝑥)−𝑠(𝑥)𝑤₁(𝑥)𝑑𝑥. (3.7) Since (3.7) holds for arbitrary 𝑤₁(𝑥), we can deduce that 𝑑𝐸𝐼⁴𝑤₀𝑑𝑥⁴𝑑+𝑃²𝑤₀𝑑𝑥²=𝑞(𝑥)+𝑠(𝑥). (3.8)

In our model the roof is made of horizontal layers of rock each acting as a beam. The horizontal force 𝑃 acting on each end of the beam arises as a result of the compressive stresses due to the rock mass above. The quantity 𝑞 is the weight per unit length of the beam which we assume is constant. The quantity 𝑠 is the magnitude of the applied normal surface traction per unit length due to the transfer of stresses from the adjoining layers. We assume that 𝑠 is constant. The displacement depends on the weight and on the applied normal surface traction in the same way. We will therefore denote the combined forces, 𝑠 and 𝑞, simply by 𝑞. The above derivation and in [2] differ from that of Segal and Handelman [10] by the inclusion of 𝑞(𝑥) and 𝑠(𝑥) in the analysis. Also, in [10], the boundary conditions were as follows: 𝑤(0)=0,𝑤^′′(0)=0,𝑤(1)=0,𝑤^′′(1)=0, (3.9) while in [2], the boundary conditions were as follows: 𝑤(0)=0,𝑤^′(0)=0,𝑤(1)=0,𝑤^′(1)=0. (3.10)

We see that the Euler-Bernoulli beam equation remains valid for the boundary conditions used in this paper. Equation (3.8) is now written in dimensionless form. Define [2]: 𝑥𝑥=𝐿,𝑤𝑤=𝑆,𝑆=𝑞𝐿⁴,𝐸𝐼 (3.11) where 𝑆 is the characteristic displacement. Equation (3.8) becomes 𝑑⁴𝑤𝑑𝑥⁴+𝐵²𝑑²𝑤𝑑𝑥²=1, (3.12) where the beam number 𝐵 is defined by (2.1). The boundary conditions, (3.2), when expressed in dimensionless variables become 𝑤(0)𝑑=0,𝑤𝑑𝑥(0)=0,𝑤(1)𝑑=0,²𝑤𝑑𝑥²(1)=0. (3.13)

The bending moment 𝑀 is [10] as follows: 𝑑𝑀=𝐸𝐼²𝑤𝑑𝑥², (3.14) where, since it is assumed that the displacement is sufficiently small that linear theory applies [10], 𝑑²𝑤𝑑𝑥²=curvatureoftheneutralaxisofthebeam. (3.15)

Define [2]: 𝑀𝑀=𝑞𝐿². (3.16)

Then 𝑑𝑀=²𝑤𝑑𝑥², (3.17) and we will refer to 𝑀 as both the bending moment and curvature of the beam. The overhead bar will be suppressed in the rest of the paper to keep the notation simple.

4. Mathematical Solution

Consider the model of the roof rock between a pillar at 𝑥=0 which is the working face and a snook at 𝑥=1 described by an Euler-Bernoulli beam with end 𝑥=0 clamped and end 𝑥=1 simply supported or hinged. The displacement 𝑤(𝑥) satisfies the differential equation 𝑑⁴𝑤𝑑𝑥⁴+𝐵²𝑑²𝑤𝑑𝑥²=1, (4.1) subject to the boundary conditions 𝑤(0)=0,𝑑𝑤(0)(1)𝑑𝑑𝑥=0,𝑤=0,²𝑤𝑑𝑥²(1)=0. (4.2)

The solution for 𝑤(𝑥) is 𝑤(𝑥)=𝑥²2𝐵²−1𝐵⁴()𝐵sin𝐵−𝐵cos𝐵1+²2(+𝐵sin𝐵−𝐵(1−cos𝐵𝑥))1−1+²2cos𝐵(𝐵𝑥−sin(𝐵𝑥)), (4.3) provided that 𝐵 does not satisfy 𝐹(𝐵)=sin𝐵−𝐵cos𝐵=0. (4.4)

We can rewrite (4.4) as follows: tan𝐵=𝐵. (4.5)

Equations (4.4) and (4.5) are plotted in Figures 5 and 6. The first five roots of (4.4) and (4.5) are 𝐵₀=0,𝐵₁=4.4934,𝐵₂=7.7253,𝐵₃=10.9041,𝐵₄=14.0662. (4.6)

Figure 5 

Graph of 𝐹(𝐵)=sin𝐵−𝐵cos𝐵 against 𝐵 for 𝐵 in the range [0,10].

Figure 6 

Graph of 𝑦=tan𝐵 and 𝑦=𝐵 against 𝐵 for 𝐵 in the range [0,10].

At these values of 𝐵, the displacement becomes infinite. Since 𝐵=4.4934 is the first nonzero value of 𝐵 for which the displacement becomes singular, our primary concern is in the interval 0<𝐵<4.4934. In the next section, we discuss the solution (4.3) for the displacement and calculate the curvature of the beam which determines the location at which the beam will break.

For a beam with both ends clamped [2], the displacement becomes singular for 𝐵=2𝑛𝜋 where 𝑛=1,2,3,…. Comparing the first points at which the displacement becomes singular in the two models, we note that 4.4934<2𝜋=6.2832. Small displacements and small derivatives are used in the derivation of the Euler-Bernoulli beam equation. The theory therefore breaks down in the neighbourhood of the points 𝑥=𝐵_𝑛 where the curvature has singular behaviour. A full nonlinear theory would need to be used in these regions. However, the beam will break when its tensile strength is exceeded which could be well before the singularities in the curvature are reached.

5. Analysis of the Results

Graphs of the displacement 𝑤(𝑥) for values of the beam number, 𝐵, in the range 0<𝐵<4.4934 are shown in Figure 7. The displacement has two stationary points which are located at 𝑥=0 and at an interior point.

Figure 7 

Graphs of the displacement 𝑤(𝑥) for 𝐵=0.5,1.5,2.5,3.5,4.

Consider first the beam for small values of 𝐵. The asymptotic expansion of 𝑤(𝑥) as 𝐵→0 is given by the following: 𝑤(𝑥)=𝑥²(1−𝑥)𝐵483−2𝑥+²604𝑥³−11𝑥²+4𝑥+6+𝑂𝐵⁴. (5.1)

The displacement is nonzero when 𝐵=0 because of the weight per unit length, 𝑞, acting on the beam. Graphs of the displacement 𝑤(𝑥) for small values of 𝐵 are presented in Figure 8. From the graphs, we can see that (5.1) is a good approximation for the displacement for 𝐵≤0.9.

Figure 8 

Comparison of the asymptotic expansion (5.1) correct to order 𝐵2 (- - -) with the exact solution (4.3) (—) of the deflection 𝑤(𝑥) for 𝐵=0.1,0.25,0.5,0.9. The curves for the exact solution and asymptotic solution overlap.

Denote by 𝑥₀ the point of maximum displacement of the beam. In order to estimate 𝑥₀ for small 𝐵, consider, from (5.1), 𝑑𝑤=𝑥𝑑𝑥488𝑥²𝐵−15𝑥+6−²208𝑥⁴−25𝑥³+20𝑥²+2𝑥−4+𝑂𝐵⁴=0,as𝐵⟶0. (5.2)

The root of the quadratic equation 8𝑥²−15𝑥+6=0 (5.3) in the range 0<𝑥<1 is 𝑥=0.578 and therefore, 𝑥₀=0.578+𝑂𝐵²,as𝐵⟶0. (5.4)

The root of 8𝑥⁴−25𝑥³+20𝑥²+2𝑥−4=0 (5.5) in the range 0<𝑥<1 is 𝑥=0.594 which is only 2.77% larger than 0.578. The maximum turning points of the curves in Figure 10 are all close to the zero order in 𝐵 value. This shows that (5.4) is a good approximation of 𝑥₀ for small values of 𝐵. Substituting (5.4) into (5.1) gives 𝑤𝑥₀=0.0054+𝑂𝐵²,as𝐵⟶0, (5.6) which from Figures 7 and 8 is a good approximation for 𝐵≤0.5. In comparison, for a beam with clamped ends, from symmetry, the magnitude of the deflection is a maximum at 𝑥₀=0.5 for 0<𝐵<2𝜋.

In order to gain insight into the possible failure of the beam, we need to determine the point at which the beam is under maximum stress. We assume that this is the point at which the magnitude of the curvature is greatest. The dimensionless curvature of the beam is given by the following: 𝑑²𝑤𝑑𝑥²=1𝐵²−1𝐵²()𝐵sin𝐵−𝐵cos𝐵1+²2()+𝐵sin𝐵−𝐵cos𝐵𝑥1−1+²2cos𝐵sin(𝐵𝑥). (5.7)

Graphs of the magnitude of the curvature for a range of values 𝐵<𝐵₁ are given in Figure 9.

Figure 9 

Graphs of the magnitude of the curvature for 𝐵=0.5,1.5,2.5,3.5,4.

Figure 10 

Graph of the magnitude of the curvature for 𝐵=4.4930.

The curvature vanishes at 𝑥=1 because the end 𝑥=1 is simply supported. Denote by 𝑥₁ the position of the local maximum of the magnitude of the curvature when 0<𝑥<1. Since the sign of the curvature at 𝑥=0 is opposite to that at 𝑥=𝑥₁, the curvature must vanish at a point, say 𝑥₂, where 0<𝑥₂<𝑥₁. We can deduce from Figure 9 that 𝑥₂ does not change significantly as 𝐵 is increased. For the range of values used in Figure 9, the magnitude of the curvature at 𝑥=0 is greater than at 𝑥=𝑥₁. We will investigate later as to whether this is always the case.

Consider first the curvature for small values of 𝐵. The asymptotic expansion of (5.7) as 𝐵→0 is given by the following: 𝑑²𝑤𝑑𝑥²=18(𝐵1−𝑥)1−4𝑥+²3010𝑥³−15𝑥²+1+𝑂𝐵⁴. (5.8) To zero order in 𝐵, the zeros of the curvature occur at 𝑥=1 and 𝑥=1/4 and therefore, 𝑥₂=0.25+𝑂𝐵²,as𝐵⟶0. (5.9)

The root of the cubic equation 10𝑥³−15𝑥²+1=0 (5.10) in the range 0<𝑥<1 is 𝑥=0.287 which explains why 𝑥₂ does not greatly depend on 𝐵. The expansions are in good agreement with Figure 9.

We now consider the turning point 𝑥=𝑥₁ of the curvature. From (5.7), 𝑑³𝑤𝑑𝑥³=1𝐵²()𝐵sin𝐵−𝐵cos𝐵1+²2()−𝐵sin𝐵−𝐵sin𝐵𝑥1−1+²2cos𝐵cos(𝐵𝑥). (5.11) which vanishes for 𝑥=𝑥₁, where tan𝐵𝑥₁=𝑃(𝐵)𝑄(𝐵),𝑃(𝐵)𝐵=1−1+²2𝑄(𝐵)=𝐵cos𝐵,1+²2sin𝐵−𝐵. (5.12)

We now examine 𝑅, the absolute value of the ratio of the curvature at 𝑥=0 to the curvature at 𝑥=𝑥₁ as follows: 𝑅(𝐵)=𝑤^′′(0)𝑤^′′𝑥₁. (5.13)

For the range of values of 𝐵 considered in Figure 9, 𝑅(𝐵)>1. However, as 𝐵 approaches the first singular value 𝐵=𝐵₁=4.4934, the ratio 𝑅 approaches and may exceed unity. This is illustrated in Figure 10 where 𝐵=4.4930. We now investigate analytically the ratio 𝑅(𝐵) for 0≤𝐵≤𝐵₁.

Consider first the asymptotic behaviour as 𝐵→0. Now 𝑃(𝐵)=5𝐵24⁴+𝑂𝐵⁶,𝑄(𝐵)=13𝐵³+𝑂𝐵⁵, (5.14) and therefore, from (5.12) 𝑥₁=58+𝑂𝐵²,as𝐵⟶0. (5.15)

The value 𝑥₁=0.625 for small 𝐵 is consistent with the graphs in Figure 9. Also, from (5.8), as 𝐵→0, 𝑤^′′(0)=181+𝑂𝐵²,𝑤^′′𝑥₁9=−𝐵1281+𝑂², (5.16) and therefore, 𝑅(𝐵)=169+𝑂𝐵²,as𝐵⟶0. (5.17)

A more detailed calculation yields the following result 𝑅(𝐵)=1691−47𝐵2304²+𝑂𝐵⁴, (5.18) so that 𝑅(𝐵) decreases initially as 𝐵 increases.

Consider next the limit as 𝐵→𝐵₁. The curvature 𝑤^′′(𝑥), given by (5.7), may be expressed in terms of 𝐹(𝐵), 𝑃(𝐵), and 𝑄(𝐵) as follows:𝑑²𝑤𝑑𝑥²=1𝐵²𝐹(𝐵)[𝐹(𝐵)−𝑄(𝐵)cos(𝐵𝑥)−𝑃(𝐵)sin(𝐵𝑥)]. (5.19)

Then, 𝑤^′′(0)=1𝐵²𝐹(𝐵)[𝐹(𝐵)−𝑄(𝐵)]. (5.20)

Also, since 𝑥₁ is defined by (5.12), cos𝐵𝑥₁=𝑄(𝐵)𝑃²(𝐵)+𝑄²(𝐵)^1/2,sin𝐵𝑥₁=𝑃(𝐵)𝑃²(𝐵)+𝑄²(𝐵)^1/2, (5.21) and therefore, 𝑤^′′𝑥₁=1𝐵²𝐹(𝐵)𝐹(𝐵)−𝑃²(𝐵)+𝑄²(𝐵)^1/2. (5.22)

Hence, 𝑅(𝐵)=𝑤^′′(0)𝑤^′′𝑥₁=𝐹(𝐵)−𝑄(𝐵)𝐹(𝐵)−𝑃²(𝐵)+𝑄²(𝐵)^1/2. (5.23)

Now at 𝐵=𝐵₁, 𝐹(𝐵₁)=0 and therefore, sin𝐵₁=𝐵₁cos𝐵₁. Thus, 𝑄𝐵₁=−𝐵₁𝑃𝐵₁, (5.24) and hence 𝑅𝐵₁=𝐵₁1+𝐵²₁^1/2. (5.25)

Thus, 𝑅(𝐵₁)<1. Also, using (5.12) and (5.24), we find that in the limit 𝐵=𝐵₁, 𝑥₁ is given by the following: tan𝐵₁𝑥₁1=−𝐵₁. (5.26)

Since 𝐵₁=4.4934, it follows that 𝑅𝐵₁=0.9761, (5.27) and that 𝑥₁𝐵₁=0.66. (5.28)

The graph of 𝑅(𝐵) against 𝐵 for 0≤𝐵<𝐵₁ is presented in Figure 11. The analytical results for 𝐵→0 and 𝐵→𝐵₁ agree with the graph. In comparison when the two ends of the beam are clamped, 𝑅(𝐵)>1 for 0≤𝐵<𝐵₁ and 𝑅(𝐵₁)=1, where 𝐵₁=2𝜋.

Figure 11 

Graph of the ratio 𝑅(𝐵) for 0≤𝐵≤𝐵1, where 𝐵1=4.4934. The point of intersection of the two curves is 𝐵=4.4324.

6. Numerical Estimates

Consider first the beam number 𝐵 defined by (2.1). The total horizontal force 𝑃 acting on each end section is given by [2, 9, 10] the following: 𝑃=𝑘𝜌𝑔𝐻𝑏ℎ, (6.1) where 𝐻 is the depth of the mine below the surface of the earth, 𝜌 is the average density of the rock from the surface of the earth to the depth 𝐻, 𝑏 is the breadth of the roof beam, ℎ is the thickness of the beam, and 𝑘 is the lateral stress coefficient. The lateral stress coefficient is a function of the rock properties. The value 𝑘=0 corresponds to a material that is completely solid while 𝑘=1 corresponds to a fluid in which the pressure is isotropic. Some models can predict values of 𝑘>1 and that 𝑘 decreases with depth [7]. For the shallow coal mines which we will consider, we will take 𝑘=2 [11]. The second moment of area about the 𝑥₂-axis is [2] 𝐼=𝑏ℎ³.12 (6.2)

The beam number (2.1) becomes 𝐿𝐵=3.46ℎ𝑘𝜌𝑔𝐻𝐸^1/2. (6.3)

For a beam with one end clamped and the other end simply supported, the displacement and curvature become infinite first at 𝐵=𝐵₁=4.4934. Thus we obtain the upper limit 𝐿_𝑐 for the length that a beam can have without collapsing: 𝐿<𝐿_𝑐𝐸=1.2971ℎ𝑘𝜌𝑔𝐻^1/2. (6.4)

As the beam becomes more fractured with time the Young’s modulus 𝐸 will decrease and 𝐿_𝑐 will decrease.

When both ends of the beam are clamped, 𝐵₁=2𝜋 and the ratio of the maximum length 𝐿_𝑐(2) when one end is clamped and one is simply supported, to the maximum length 𝐿_𝑐(1) when both ends are clamped is𝐿_𝑐(2)𝐿_𝑐(1)=𝐵₁(2)𝐵₁(1)=0.7151. (6.5) The critical length 𝐿_𝑐(1) describes the initial stage of the process of pillar extraction when several pillars are removed and the roof is left to collapse, while 𝐿_𝑐(2) models the second stage when adjacent pillars are mined and smaller sections of the roof are left to collapse. We see that 𝐿_𝑐(2)<𝐿_𝑐(1) which is consistent with the two models.

We consider a beam made of sandstone. We use the following estimates: 𝐸=3×10¹⁰Pa,ℎ=0.1m,0.2m,…,1m,𝑘=2,𝜌=2.5×10³kg/m³,𝑔=9.8m/s²,𝐻=120m,500m,1000m. (6.6)

Table 1 summarizes numerical estimates for 𝐿_𝑐 for a beam with one end clamped and the other simply supported. In a coal mine the distance between the pillars ranges from 5 m and 7 m. When the bending moment or curvature exceeds the tensile strength of the beam the roof will collapse. This may occur for beam lengths 𝐿 less than 𝐿_𝑐 since 𝐿_𝑐 provides only an upper limit on the length of the beam for collapse.

Table 1 

Critical length 𝐿_𝑐.

The axial force 𝑃 may experience a sudden increase due to a seismic event which could last for a short time. Using (2.1), the upper limit for the length, 𝐿_𝑐, can be written as follows: 𝐿_𝑐=𝐵₁𝐸𝐼𝑃^1/2. (6.7)

If 𝐿_𝑐 is reduced below about 6 m due to an increase in 𝑃 then a roof collapse may occur.

We now compare the effect of hinged and clamped supports at the snook at 𝑥=1 on the curvature at the pillar at 𝑥=0. When both ends of the beam are clamped [2], 𝑤₁^′′(0)=1𝐵²𝐵1−(),𝑤2tan𝐵/2_1′′(0)=1𝐵12+𝑂²,as𝐵⟶0. (6.8)

When the end 𝑥=0 is clamped and the end 𝑥=1 is simply supported (hinged), from (5.7) and (5.8), 𝑤₂^′′(0)=1−cos𝐵−(1/2)𝐵sin𝐵𝐵(),𝑤sin𝐵−𝐵cos𝐵_2′′(0)=18+𝑂𝐵²,as𝐵⟶0. (6.9) Thus, 𝑤₂^′′(0)𝑤₂^′′(0)|_𝐵=0=32. (6.10)

In Figure 12 the graphs of 𝑤₁^′′(0) and 𝑤₂^′′(0) are plotted against 𝐵. We see that when the end 𝑥=1 is simply supported the curvature at 𝑥=0 is greater for a given value of 𝐵 than when it is clamped. The tensile strength of the beam will be exceeded at the end 𝑥=0 for lower values of 𝐵 when the end 𝑥=1 is simply supported than when it is clamped. The simply supported boundary condition has the effect of increasing the bending moment at the end 𝑥=0 and causing the beam to break at 𝑥=0 for lower values of 𝐵.

Figure 12 

Graphs of the curvature 𝑤′′(0) at 𝑥=0 plotted against 𝐵 for a beam with: (1) 𝑥=0 clamped and 𝑥=1 clamped, (2) 𝑥=0 clamped and 𝑥=1 simply supported.

7. Values of the Beam Number 𝐵 Greater Than 𝐵₁

A graph of the magnitude of the curvature of the beam at the end 𝑥=0, plotted against 𝐵, is given in Figure 13. It divides the values of 𝐵 into the intervals, 𝐼₁=0,𝐵₁,𝐼₂=𝐵₁,𝐵₂,𝐼₃=𝐵₂,𝐵₃,…,𝐼_𝑛=𝐵_𝑛−1,𝐵_𝑛, (7.1) where 𝐵₁,𝐵₂,𝐵₃,… are the values at which the displacement and curvature become infinite. We have only considered the first interval 𝐼₁. Since, from (2.1), we see that the beam number 𝐵 is proportional to 𝐿 and 𝑃^(1/2) the value of 𝐵 would increase if either 𝐿 or 𝑃 were to increase. The length of the beam increases by a finite amount when a snook fails. In the second stage of the pillar extraction process, small sections of the mine collapse but if the snooks are too strong they will support a longer section of the roof which will form a beam and collapse when the snooks fail. Another way in which 𝐵 could increase suddenly is due to a seismic event which may produce a discontinuous increase in 𝑃 which could last for a short period of time.

Figure 13 

Graph of the magnitude of the curvature at the end 𝑥=0 plotted against 𝐵.

Consider first the displacement. The displacement for values of 𝐵 in the first interval, 𝐼₁, was considered in Figure 7. In Figure 14 graphs of the displacement for representative values of 𝐵 in the intervals, 𝐼₁ to 𝐼₆, are presented. We see that as 𝐵 increases through the intervals the number of turning points increases and that the displacement can take negative values beyond interval 𝐼₁. The amplitude of the displacement will depend on how close 𝐵 is to the singular end points of the interval.

Figure 14 

Graphs of the displacement for representative values of the beam number 𝐵 in the intervals 𝐼1 to 𝐼6.

Consider next the magnitude of the curvature. In Figure 9 graphs of the magnitude of the curvature were plotted against 𝑥 for values of 𝐵 in the first interval 𝐼₁. In Figure 15, graphs of the magnitude of the curvature are plotted for the same representative values used to plot the displacement in Figure 14. For the values of 𝐵 and 𝑛 considered, the number of local maxima of the magnitude of the curvature in the 𝑛th interval 𝐼_𝑛 is 𝑛. The greatest local maximum is not at the end 𝑥=0 but at interior points. There may be several points for which the magnitude of the curvature has the maximum value. If the bending moment exceeds the tensile strength the beam will break at these interior points. The magnitude of the curvature depends on how close 𝐵 is to the singular endpoints of the interval. Since 𝐵 is proportional to the length of the beam 𝐿 we see from Figure 13 that if the value of 𝐵 is in the range outside of 𝐵₁ then a longer beam could be less susceptible to failure than a shorter beam. This could be associated with the beam taking on a higher mode of bending. We see from Figure 14 that the displacement can be negative. For this to be possible in practice the beam would have to be detached sufficiently from the layers above.

Figure 15 

Graphs of the magnitude of the curvature for representative values of the beam number 𝐵 in the intervals 𝐼1 to 𝐼6.

8. Conclusions

We investigated the possible roof collapse between the next pillar to be mined and the first line of snooks by modelling the roof as an Euler-Bernoulli beam. The beam was simply supported (hinged) at the snook and clamped at the pillar which was the working face. The model contained one dimensionless number—the beam number 𝐵.

Numerical estimates obtained for the critical length 𝐿_𝑐 are comparable to the expected distance between pillars which in a coal mine is 5–7 m thus making the model credible. The model also predicts that the roof may collapse at the clamped end at the pillar for 0<𝐵<4.4324. For the range 4.4324<𝐵<𝐵₁=4.4934, the model predicts that the roof may collapse at an interior point closer to the snook than to the pillar where the magnitude of the curvature attains its maximum value. However, this range contributes only 1.36 percent of the range [0,𝐵₁] and since it is likely that the threshold of the stress would have been exceeded for values of the beam number below 𝐵=4.4324, for all practical purposes, the beam will break at the clamped end if the threshold of its stress is exceeded. It will therefore break at the pillar which is the working face. However, the model showed that it is not necessary for the beam to break at the clamped end. The beam may break at interior points. In 𝐼₁, the interior point is unique for a given value of 𝐵. Other boundary conditions would need to be considered and analyzed in order to determine whether the beam can collapse at an interior point for practical values of 𝐵.

From Figure 12 we deduced that a beam which is simply supported (hinged) at the snook produces a larger bending moment at the pillar (where it is clamped) than a beam which is clamped at the snook. The beam will break at a lower value of the beam number when it is hinged at the snook.

The displacement and curvature become infinite at the zeros of 𝐹(𝐵) defined by (4.4) which divides the value of 𝐵 into intervals. We considered mainly the first interval 0≤𝐵≤𝐵₁. However, as the snooks fail 𝐵 can increase discontinuously by finite amounts and may take values in the higher intervals. A preliminary computational investigation was undertaken of the displacement and curvature for values of 𝐵 in these intervals. It was found that the displacement can take negative values. In the 𝑛th interval the magnitude of the curvature had 𝑛 local maxima. The maximum value of the magnitude of the curvature did not in general occur at the pillar but occurred instead at interior points. If the tensile strength of the beam is exceeded the beam could break at several interior points.

In practice, roof bolts are installed [12].

Acknowledgments

The problem of pillar extraction in coal mines was submitted to the Mathematics in Industry Study Group 2011 (MISG) at the University of the Witwatersrand, Johannesburg. The author thanks the referees for their valuable comments which improved this paper greatly. The author thanks Matthew Woolway, University of the Witwatersrand, Johannesburg, for assisting in creating the diagrams. The author would like to express gratitude and appreciation to Professor Nielen van der Merwe and Dr. Halil Yilmaz of the School of Mining Engineering, University of the Witwatersrand, for submitting the problem. The author would particularly like to acknowledge that she has profited greatly from the discussions held with Professor Nielen van der Merwe. The author would like to acknowledge Professor Colin Please of the University of Southhampton, England, who formulated the model for pillar extraction and roof collapse at the MISG. The author thanks Professor David Mason for his valuable comments, advice, and encouragement while she was writing this paper.