**TRANSACTION PAPER**

**Numerical computation of average pillar stress and implications for pillar design **

**J.A.L Napier ^{I, II}; D.F. Malan^{II}**

^{I}CSIR, South Africa and Department of Mining Engineering, University of Pretoria, South Africa

^{II}Department of Mining Engineering, University of Pretoria, South Africa

]]>

**SYNOPSIS**

A number of issues relating to the computational aspects of pillar design are addressed in this paper. The computation of average pillar stress values is important when attempting to establish criteria for pillar design and in the analysis of the stability of tabular pillar layouts. One of the default 'classic' numerical methods that are used to determine pillar stresses is the displacement discontinuity method. In many instances it is not clearly understood that this approach does suffer from some limitations, particularly in relation to the fact that in coarse element simulations, the simulated average pillar stress (APS) can depend on the chosen mesh size. The nature of this error is highlighted in this paper and some strategies are suggested to bound the magnitude of these errors. It is demonstrated as well that the popular linear stiffness approximation to pillar or seam compressibility does appear to allow reasonably accurate estimates of the average pillar stress when either the pillar height is varied or when the seam modulus differs from the host rock modulus. A practical implication of this study is that if the seam modulus is noticeably lower than that of the host rock, such as for coal seams, it is important to use a linear stiffness constitutive model for the pillars rather than a 'rigid' pillar assumption. This added complexity seems unnecessary, however, when simulating hard rock pillars in mines where the seam modulus is very similar to that of the surrounding rock.

**Keywords:** pillar design, average pillar stress, numerical modelling, tabular layouts

**Introduction**

Pillars are ubiquitously used in hard rock underground mining operations for the control of regional layout stability, and in some cases in the form of crush pillars or yield pillars to control local roof conditions. Although pillar behaviour has been studied^{1-3} for many decades, the methodologies used to design pillars are nevertheless still imperfect and require additional research^{4}. A particular difficulty in pillar design is that pillar failure occurs according to a number of complex local mechanisms imposed by the heterogeneity of the pillar geological structure (Figure 1 illustrates one of these failure mechanisms) and by the deformations in the proximate roof and floor regions. It has also been understood for many years that the overall stability of a pillar layout system will depend on the relative stiffness characteristics of the surrounding host region and on the effective loaddeformation (strain-softening) characteristics of the pillar units^{5,6}. Owing to the many unknowns, empirical determination of pillar strength has remained a popular approach^{3}.

]]>

The empirical determination of pillar strengths^{7-9} usually depends on classifying case studies of selected visual observations and monitored deformations of actual pillars against an estimate of the load that is carried by these pillars. The usefulness of these analyses is strongly dependent on the accuracy of the estimated pillar load. Before numerical modelling programs were commonly available, tributary area theory was frequently used to estimate the pillar loads. Salamon and Munro7 and Hedley and Grant1 used this approach to estimate pillar loads and to derive their famous empirical pillar strength formulae.

A difficulty with tributary area theory is that it cannot take irregularities, such as solid abutments and large barrier pillars, into account. Ryder and Jager^{10} suggest that the deviations are small at shallow depth in regular pillar layouts having low extraction ratios, and can in fact be ignored in most aspects of coal mine design. In contrast, in hard rock mines, the deviations can be significant owing to the much higher extraction ratios. These layouts require numerical modelling to evaluate the average pillar loads^{11,9}. In estimating pillar strength, a typical approach is to identify pillars underground that are very small or are considered to have failed. By careful numerical simulation of these layouts and pillars, the strength of the pillars, or at least the minimum strength of intact pillars, can be estimated.

Although numerical techniques have become popular for estimating pillar stresses, these techniques may lead to erroneous results if the basic characteristics of these methods are not fully appreciated. As pillar load is such a fundamental building block in pillar strength studies, the objective of this paper is to create awareness regarding some of the pitfalls that arise when calculating average pillar stresses. In particular, the displacement discontinuity boundary element method (DDM) is a numerical technique that is very well suited to the analysis of pillar systems in large-scale tabular layout problems in which the plan projection of the mined area may cover several square kilometres. This paper discusses some of the numerical errors that arise when using the DDM to calculate average pillar stresses. It must be noted that while local failure and deformation mechanisms may be analysed using nonlinear continuum plasticity finite element or finite difference models, these methods are themselves difficult to apply to the interactive analysis of extensive rock mass regions surrounding large-scale three-dimensional tabular pillar layouts.

]]>**Solution of tabular mine layouts using the displacement discontinuity boundary element method**

The displacement discontinuity boundary element method assumes that the mining height (stope width) is generally much smaller than the areal extent of mining, and the excavation is approximated as a simple cut in the medium having two opposing surfaces. This approximate representation of the tabular excavation is termed a 'slit' approximation in the remainder of the paper. The relative difference in the displacement vector at each point on the excavation 'slit' surface determines the induced displacement and stress components at all points in the medium according to a wellestablished influence function theory, which allows the redistribution of stress following mining to be calculated. This approach was suggested originally for determining stress concentrations and surface subsidence deformations induced by coal mining operations in the United Kingdom^{12,13}. A later version of this procedure was introduced in South Africa by Salamon^{14}, who termed the method the 'face element principle'. The method has subsequently become known as the displacement discontinuity method^{15,16}. A variety of computer tools have been developed in a number of countries to solve large-scale tabular mine layout problems using the displacement discontinuity technique. Specific examples include MSCALC/ BESOL (Crouch^{15}), MAP3D (Wiles^{17}), MINSIM-D (Ryder and Napier^{18}), TEXAN (Napier and Malan^{19}), and special purpose tools for analysing laminated strata such as LaModel (Salamon^{20}, Heasley^{21}). In the DDM analysis the plan outline of the mined area is tessellated using square, quadrilateral, or triangular shaped elements and a specific displacement discontinuity variation is assumed within each element. The problem is solved to determine the detailed distribution of the displacement discontinuity components in all elements such that the total stress at designated points within each mined element is zero (corresponding to the boundary condition that would be expected in an open excavation) or is constrained by the limiting movement of the roof and floor when complete elastic convergence (closure) occurs. The DDM has provided a basic 'workhorse' tool to rock mechanics practitioners for routine estimation of stress distributions in pillar layouts and for assessing the potential for induced fault slip in deep-level gold and platinum mine workings. One of the most important applications of the DDM is to determine pillar stress distributions and average pillar loads for the assessment of layout stability or in the back-analysis of observed pillar failures, and the empirical determination of pillar strength criteria.

There are unfortunately a number of computational difficulties associated with the DDM technique, which are not generally appreciated and which may lead to misinterpretations of analysis results. The reasons for these difficulties and some solutions to these problems are presented in the following sections. The discussion is confined to the consideration of only the stress component normal to the plane of the seam or reef horizon.

**Application of tributary area theory to determine average pillar stress**

Let the region that is enclosed by the projection of a specific pillar outline onto the seam or reef plane be designated by Ω. Let the area of the projected pillar region in the seam or reef plane be equal to *R, *and assume that *z *defines the coordinate direction normal to this plane. Let σ* _{zz}* (

*P*) be the value of the stress component normal to the reef plane at point

*P*within the pillar region Ω. This stress component will be assumed to be compressive and, following the sign convention used in this paper, will be negative. The absolute magnitude,

*A*, of the average pillar stress component,

_{P}*, acting over the pillar is then defined to be given by the following expression:*

_{zz}In Equation [1], *dR _{p}* denotes the differential area element with respect to the points

*P*in the pillar region Ω. The evaluation of Equation [1] over a general pillar shape is, in general, very difficult and in most practical cases has to be determined numerically. However, in the particular case where each pillar is a uniformly replicated unit within an infinite, horizontal tabular layout, the average pillar stress,

*, is given by*

_{zz}*e*is the extraction ratio and is the virgin vertical stress at the seam horizon. This relationship provides a useful check for numerical procedures and is a good approximation to the average stress acting on pillars near the centre of an extensive, regular layout. In a layout comprising pillars of different sizes, Equation [2] can be interpreted as a general statement of force equilibrium before and after mining takes place, and expresses the average stress acting on the entire pillar ensemble and on the abutment regions. This average value determination is termed the simplified tributary area theory (TAT)

^{10}and is independent of the surrounding rock mass modulus. It must be emphasised that each individual pillar in a general layout will have a stress value that may be either above or below the average TAT value expressed by Equation [2], and the TAT value is therefore not necessarily a conservative estimate for the average stress over every pillar in the layout. The focus in the present paper is on the numerical computation of the average stress in individual pillars as defined by Equation [1].

Figure 4 illustrates a horizontal test layout simulated with the TEXAN code^{19} where the nominal extraction ratio is *e* = 0.75 within the pillar region and the virgin stress at the reef horizon is = -108 MPa. Using Equation [2], the tributary area estimate of the average pillar stress is *A _{P}* = 432 MPa. Figure 5 illustrates the computed average stress values for a row of pillars (from left to right) in the centre of the excavation. The layout is tessellated using 10 m square elements with constant displacement discontinuity values, and each pillar is represented by a 4 x 4 cluster of elements. The average stress in each pillar is calculated as the mean of the stress values that arise at the centres of the 4 x 4 pillar element clusters. Figure 5 illustrates that as the overall size of the layout area increases, the average pillar stress tends towards the tributary area stress asymptote predicted by Equation [2]. The average pillar stress values of the pillars adjacent to the abutments can be seen to be significantly lower than the asymptotic value. This is an illustration of how the nominal tributary area estimate of the average stress overestimates the stresses on pillars, since the effect of abutments or barrier pillars is ignored. Different results will arise if the pillar sizes are variable, or if irregular layout configurations are analysed. In these cases the numerically computed average pillar stress values in individual pillars may exceed or fall below the TAT value estimated by Equation [2]. It is also important to understand that certain numerical errors are associated with the use of the DDM, and that the choice of element size can affect the computed estimates of average pillar stress values. This is discussed in the next section.

]]>

**Effect of element size on the determination of average pillar stress values**

Errors can arise in the application of the DDM to the numerical determination of average stress values of individual pillar, depending on the choice of element size and the element shape in a geometrically complicated layout configuration. Some perspective on the nature of these errors can be gained by considering initially the simple case of a single strip pillar located between two similar parallel-sided panels in a horizontal plane. The layout configuration and characteristic dimensions are shown in Figure 6.

If the displacement discontinuity slit approximation is used, it is possible to determine the average pillar stress analytically (Salamon^{22}). The resulting expression is

where *K*(*k*) is the complete elliptic integral of the first kind defined as^{23}

^{22}

The pillar problem illustrated in Figure 6 may be solved numerically by dividing each panel into equal-sized, constant value displacement discontinuity elements. The total stress σ*zz*^{t }(*s _{i}*) at the centre

*si*of panel element

*i*is given by the expression:

where *D _{z}*(

*s*) is the value of the assumed constant elastic convergence in element

_{j}*j*and where the sum is taken over all elements,

*j*= 1 to

*N*. The sign of

*D*(

_{z}*s*) in Equation [6] is positive when the relative movement of the seam roof (reef hangingwall) is towards the seam floor (reef footwall). The influence kernel function

_{j}*F*(

_{zz}*s*) is given by

_{i}, s_{j}In Equation [7], *g _{j}* is the length of element

*j, E*is the Young's modulus of the country rock, and

*v*is the Poisson's ratio (a derivation of Equation [7] may be found in Crouch and Starfield

^{16}). It should be noted that the expression for

*F*(

_{zz}*s*) depends only on the relative distance |

_{i}, s_{j}*s*-

_{i}*s*| between the centres of elements

_{j}*i*and

*j*. In the special case where all element sizes have the same value

*g*, it is apparent that

*s*-

_{i}*s*= (

_{j}*i*-

*j*)

*g*and Equation [7] reduces to

∞ Since it is known that it follows that indicating that the influence kernels are self-equilibrating. The numerical solution to Equation [6] is

obtained by setting the total stress(*s _{i}*) at the centre of each panel element to zero and solving for the unknown values of the displacement discontinuity values

*D*(

_{z}*s*). In the specific case where there is only a single element of length

_{j}*g*at position

_{i}*si*, the unknown value of

*D*(

_{z}*s*) in Equation [6] is given by

_{i}It is of interest to note that the analytical solution for thee average elastic convergence *D _{z} *in a single, open parallelsided panel

^{22}of span

*g*is given by

which can be seen to be exactly one half of the value given by Equation [9]. Although a single, constant strength displacement discontinuity value is admittedly a very crude approximation to the elastic convergence in an isolated panel, it is apparent that the use of the kernel expressions given by Equations [7] or [8] can be expected to yield some inaccuracy in estimations of the elastic convergence that will depend on the choice of the element size. It is important also to note that since the influence kernel values are self-equilibrating, the element grid size effect will be reduced if extensive replicated layouts are considered, such as the case illustrated in Figure 4. In practice, all layouts are obviously of finite extent, and consequently errors in estimated stress values can arise when the element sizes are significant compared to the characteristic dimensions of the pillars that are of interest. A detailed analysis of this error behaviour and strategies to compensate for the grid size effect has been presented by Ryder and Napier^{18} but, unfortunately, this treatment does not address the question of how to calculate average pillar stress values. This issue is considered in more detail in the present paper.

In order to gain some perspective on this question, consider a specific example in which the pillar and panel span dimensions, *W* and *S* in Figure 6, are assigned the values *W* = 24 m and *S* = 120 m respectively. Using Equations [3], [4], and (5), it can be determined that the analytical value of the

average pillar stress is 4.55908. A numerical estimate of the average pillar stress can be obtained in two stages. In the first stage, the panel and pillar regions are divided into uniformly sized elements *g _{i}=g* (where, for convenience,

*g*is chosen to be such that

*W*and

*S*are a common integral multiple of

*g*) and Equation [6] is solved for the elastic convergence values

*D*(

_{z}*si*) in each panel. In the second stage, the average stress in the central pillar is found by again employing Equations [6] and [8] to determine the total stress values (

*sk*) at the centre

*s*of each pillar element,

_{k}*k*, using the known solution values

*D*(

_{z}*s*) and the kernel influence functions given by Equation [8]. The

_{i}numerical estimate of the average pillar stress, is calculated from the weighted average summed over the pillar elements:

Equation [11] reduces to the arithmetic average if all pillar elements have the same size. The results of these computations of the average pillar stress, when using a virgin stress value of 100 MPa, are plotted in Figure 7 as a function of the element size, *g*.

It is striking that the estimated absolute magnitude of the average pillar stress values, shown in Figure 7, follow a nearly perfect linear trend when plotted against the element grid size. The fitted trend intercept has a value of 456.05 MPa, which may be compared to the analytic solution value of 455.91 MPa (an error of about 0.03%). The linear trend displayed in Figure 7 also suggests that if two APS values of *A*_{1} and *A*_{2} are found, corresponding to grid size values of *g*_{1} and *g* _{2} respectively, then the extrapolated value (for the APS limit as the grid size tends to zero) can be estimated from the equation

In the special case where *g* 2= *g* _{1}/2, Equation [12] reduces to

where *A*_{2} corresponds to the APS estimate with the finer grid size.

In practice, it is not possible to arrange the pillar and panel element sizes to conform exactly to the edges of the layout excavations and to represent the detailed outlines of irregular pillar shapes using uniform square element tessellations. In this case, the best that can be done is to tessellate the layout with elements that are as evenly sized as possible, and to designate the element as being 'mined' if the open excavation constitutes more than one half of the element area. A finer tessellation may then be constructed either by splitting the coarse parent tessellation or by re-meshing the problem region with a finer element tessellation. Consider, for example, an annular excavation region in the *x-y* plane surrounding a central, circular pillar as shown in Figure 8. Suppose that a square region of side length 300 m is tessellated with square, constant variation displacement discontinuity elements having side lengths *g*. The centre points of the elements are assumed to be located at positions *x _{p} = x0* +

*pg*and

*y*where

_{q}= y0*p*and

*q*are integers and (

*x0,y0*) is a chosen coordinate origin. The stress component values normal to the excavation plane are given at the centres of the elements by an appropriate influence relationship analogous to Equation [6]. For the case of constant value displacement discontinuity elements, this relationship may be expressed as

where (*x _{p}, y_{q}*) is the virgin stress at point (

*x*) and

_{p}, y_{q}*D*(

_{z}*x*) is the normal displacement discontinuity component at point (

_{r}, y_{s}*x*) on the excavation plane. The summation limits

_{r}, y_{s}*M*and

*N*cover the region of interest. The stress influence function

*F*is given by the expression

_{zz}^{14,16,19}

Were *ε*_{1} = -1, *ε*_{2} = 1,* A _{i} = g(p-r*-(

*ε*

_{i}/ 2)),

*Bj = g(q-s*-(ε

_{j }/2)) and

*R*= A close examination of Equation [15] also reveals that the kernel influence function,

_{ij}*F*, depends only on the absolute magnitude of the relative distance coordinates |

_{zz}*x*| and |

_{p}- x_{r}*y*| between the source point (

_{q}-y_{s}*x*) and receiving point (

_{r}, y_{s}*x*). It may be noted as well that the influence matrix defined by Equation [15] is also self-equilibrating. The annular excavation problem depicted in Figure 8 was solved using Equations [14] and [15] with a series of different element mesh sizes, g, that are summarized in Table I. In each case, the centre of the circular pillar was assumed to be located at

_{p}, y_{q}*x*=150,

_{c}*y*= 150 and the mesh origin was chosen to be

_{c}*x*= 0,

_{0}*y*

_{0}= 0. In Table I,

*N*is the number of elements assigned to the annular excavation and

_{E}*NP*is the number of elements in the central pillar. These numbers obviously increase as the element size g decreases. The values

*S*and

_{E}*S*indicate the total area of the elements in the excavation region and in the pillar region respectively. The asymptotic values of

_{P}*S*and

_{E}*S*, as the grid size tends to zero, are 67858.4 m

_{P}^{2}and 2827.4 m

^{2}.

Figure 9 is a plot of the estimated stress values at the centres of the elements covering the central pillar region for the different mesh sizes given in Table I. Each point is plotted as a function of the radial distance of the element centre from the centre of the pillar. The magnitudes of the average pillar stress values corresponding to each element grid size are summarized in the last column of Table I and are plotted in Figure 10. It is apparent that the trend in average values is nonlinear, unlike the plane strain case shown in Figure 7, making the extrapolation procedure more uncertain. For example, if it is assumed that a quadratic extrapolation can be applied, then the average pillar stress would be given by

]]>

]]>

where *A _{1}, A*

_{2}, and

*A*

_{3}represent the average stress magnitudes computed with element grid sizes of

*g, g*/2, and

*g*/4 respectively. Applying Equation [16] to the first three entries of the last column in Table I yields an extrapolated average stress magnitude of 429.5 MPa. If Equation [16] is applied to the last three entries, the extrapolated average stress magnitude is 434.1 MPa. Alternatively, applying the linear extrapolation of Equation [13] to the last two values of the last column of Table I yields an extrapolated average stress magnitude of 434.9 MPa. An independent check on the accuracy of this estimate was made by Gordeliy

^{24}using a special axisymmetric formulation of the displacement discontinuity method

^{25}. This yielded a value of 435.0 MPa, which is in close agreement with the extrapolated value estimates. It may be noted as well that the simple extrapolation procedures described here could be extended to a more general Richardson deferred limit extrapolation

^{26}procedure, admitting a power-law variation of the average pillar stress as a function of the element grid size, but probably offering a limited improvement in accuracy in the present case.

This example demonstrates that the extrapolation technique must be treated with some caution in the case of tabular mine layouts. The main sources of numerical error arise from both the inaccuracy of constant value displacement discontinuity element convergence values and from the limitations of using a square element shape to tessellate complex layout configurations. It is important to note that for the finest grid size used in Table I (*g *= 1.25 m) the ratio of the pillar diameter to the grid size is 48. In practical layout problems, the use of such a fine grid relative to the characteristic dimension of pillars of interest may be prohibitive computationally. Unfortunately, there are no simple numerical solutions to this problem. One possible strategy is to use higher order triangular or quadrilateral shaped elements^{19}. This certainly provides a much more efficient method of representing complex layout plan forms and, combined with the extrapolation method described here, represents the best current numerical approach to the estimation of average pillar stress magnitudes. However, it is important to note that it is not possible to determine with complete accuracy the stress profiles in pillar or abutment regions, since no exact asymptotic analytic forms are known for the elastic convergence distribution near the edges of a general polygonal-shaped pillar region or mining abutment region.

A second difficulty is that the stress distribution in an actual pillar will be affected by additional factors such as the finite height of the excavation and by local face fracturing or pillar scaling. The effect of the reef and seam stiffness is discussed in the next section.

**Effect of seam or reef stiffness on average pillar stress**

The displacement discontinuity approximation to a tabular excavation takes no account of the excavation dimensions normal to the plane of the seam or reef and it simulates the pillars with 'zero height'. The pillar can therefore not deform. A popular method that is used to compensate for the seam or reef compressibility is to employ a linear spring (the so-called 'Winkler spring' approximation) to represent the assumed elastic behaviour of the seam or reef material (see, for example, Oravecz^{27}). In Equation [14], the total stress at points (*x _{p}, y_{q}*) of the seam is assumed to be given by the expression

*is a linear stiffness factor that is used to represent the seam or reef compressibility at the point (*

_{S}*x*). For example, if the mining height is

_{p}, y_{q}*H*and if it is assumed that the lateral strain parallel to the seam or reef plane is negligible, then the stiffness

*may be deduced from Hooke's law to be*

_{S}where *E _{S} *and

*v*are the Young's modulus and Poisson's ratio respectively of the seam or reef rock in the pillar or abutment regions. It is important to note that the value of

_{S}*k*depends on the ratio

_{S}*E*and will be altered if either

_{S}/H*E*or

_{S}*H*is varied. Following the adopted sign conventions used in this paper, it apparent from Equation [17] that when

*D*> 0 the seam is compressed relative to the virgin state and the total stress becomes more compressive than the virgin stress . Substituting Equation [17] into Equation [14] yields the following equality constraint relationship for the elastic convergence

_{z}*Dz*(

*x*) at point (

_{p}, y_{q}*x*) in the seam,

_{p}, y_{q}Equations [14] and [19] can be combined to solve for the elastic convergence in both the excavation and seam regions. The total stress values at all points in the seam regions can be recovered from Equation [17].

It is of some interest to compare the stress values that are computed using the linear stiffness approximation to values that are obtained by explicitly solving two adjacent rectangular section slots as illustrated in Figure 11. The explicit slot solution is obtained using an elastic boundary element program, RIFT, and the pillar stress values are evaluated along the horizontal line AB through the centre of the pillar as illustrated in Figure 11. It is assumed that the rock mass modulus is 72 000 MPa and the Poisson's ratio is 0.2. The panel spans are set to 120 m, the pillar width is equal to 24 m, and the virgin stress at the seam horizon is assumed to be = -100 MPa. Note that the pillar and panel dimensions are similar to those used to plot the graph in Figure 7, where the analytical solution for a zero height pillar was found to be equal to 455.91 MPa and the corresponding extrapolated numerical estimate was 456.05 MPa. If the seam modulus, *E _{S}*, and seam Poisson's ratio,

*v*, are equal to the rock mass values then, using Equation [18], the linear seam stiffness is

_{S}*k*= 80 000 MPa/m if the pillar height is equal to 1 m. Figure 12 shows a comparison between the pillar stress profile determined along the centre line AB of the actual pillar and the values obtained from the displacement discontinuity slit solution using the linear stiffness factor in the pillar region. It is apparent from Figure 12 that the stress profiles are very similar in this case where the pillar height:width ratio is small (1:24). If the pillar height is increased to 4 m, the equivalent linear seam stiffness is reduced to 20 000 MPa/m. Figure 13 shows a comparison between the stress profile across the central line of the 4 m pillar and the linear stiffness model. In this case, it is apparent that the stress in the actual pillar model is much lower near the pillar edges than for the linear seam stiffness model. These diagrams illustrate some important differences in the stress profiles that are estimated using an actual pillar as compared to the simple stiffness model where the pillar edge stress is highly peaked.

_{S}

]]>

The average pillar stress values, scaled to the seam horizon virgin stress, were calculated for a series of different pillar heights and are summarized in Table II. In each case the seam modulus and Poisson's ratio are assumed to be the same as the respective rock mass values. The scaled average stress values given in Table II are plotted against the scaled pillar height values (pillar height divided by the pillar width) in Figure 14. It is of particular interest to observe that for the chosen example, where the pillar elastic parameters are the same as the surrounding rock mass, the average pillar stress values are not very strongly dependent on the height of the pillar. It is also apparent that the average stress values that are calculated across the centre of the actual pillar cases are slightly larger than but very similar to, the values for the linear stiffness model. The results plotted in Figure 13 show that the linear stiffness model may yield artificially peaked stress values near the edges of the pillar.

]]>

It is of some interest to calculate the variation of the effective vertical pillar stiffness at each point across the pillar width for different pillar heights. Let the absolute difference in the vertical displacement between the roof contact and the floor contact of the pillar at a distance *x *from the pillar centre line be designated as Referring to Equation [17], the effective point pillar stiffness, *k _{eff}*, at position

*x*within the pillar is defined to be

where the virgin vertical stress component (*x*) and the total vertical stress component(*x*) are evaluated at positions *x *along the central line AB across the pillar (see Figure 11). These effective stiffness factors were computed for each pillar height value in the first column of Table II and are plotted in Figure 15. The results indicate that the effective point stiffness profile, defined according to Equation [20], is relatively uniform in the central region of each pillar but decreases more and more sharply near the pillar edges, as the pillar height becomes smaller. The central 'plateau' stiffness in each profile shown in Figure 15 is larger than the corresponding stiffness values in the second column of Table II but falls below these values near the pillar edges. It is also noteworthy that the reduced effective stiffness near the pillar edges can be expected to induce significantly different stress patterns in the pillar foundation and roof regions near the pillar edges from the stress field that is induced by using the zero height seam or reef linear stiffness approximation defined by Equation (17). This difference in the point stiffness near the pillar edges will control the propensity for the pillar to 'punch' into the surrounding rock and will play an important role in controlling possible slip on discontinuities that are parallel to the roof and floor contacts near the pillar edges. These possibilities will not be explored further here.

From the point of view of estimating average pillar stress values, it appears that the linear stiffness model will yield comparable results to cases where the actual pillar height is explicitly modelled. If the pillar material and the host rock have similar elastic properties, then the rigid seam model where the elastic convergence is zero also apparently yields a good upper bound to the average pillar normal stress. The case where the pillar material elastic properties are different from the host rock is considered in the next section.

**Reef or seam modulus different from the host rock**

Figure16 shows two examples of underground pillars where (a) the pillar rock has approximately the same modulus as the host rock and where (b) the pillar rock modulus is softer than the host rock. The pillar shown in Figure 16a is an example of a pillar in a platinum mine in the Great Dyke in Zimbabwe, where the pillar material and the rock in the hangingwall and footwall typically have the same average moduli. In contrast, the coal pillar in Figure 16b has a much lower modulus than the host rock. It is important to consider how the average pillar stress is affected by a contrast between the pillar modulus and the host rock modulus. In addition, the important question arises whether a simple rigid pillar assumption, associated with a zero width slit solution, can lead to errors when computing the average pillar stress for the two different pillar scenarios illustrated in Figure 16.

]]>

The same pillar geometry depicted in Figure 11 and discussed in the previous section was analysed using RIFT with a seam modulus *E _{S} *that was different from the host rock modulus of 72 000 MPa. The Poisson's ratio was set equal to 0.2 for both the seam material and the host rock. The results for a fixed pillar height of 4 m are shown in Table III, and the scaled average pillar stress values are plotted against the ratio of the rock modulus to the pillar modulus in Figure 17. It is again interesting to observe that the average stress associated with the actual pillar geometry is very similar to, but slightly larger than, the corresponding average stress value determined using the linear stiffness model. This appears also to hold true if the individual values of the seam modulus

*E*and the pillar height,

_{S}*H*, are varied but the ratio of

*E*is maintained at a constant value. This is illustrated in Table IV, where the linear seam stiffness is 20 000 MPa m for each row. It is notable as well that the average pillar stress of the actual pillar is remarkably insensitive to the different pillar height cases, provided the pillar modulus is changed in proportion to the height. It should be stated again that the linear stiffness model appears to yield good results in computing the elastic average pillar stress, but should be very carefully assessed if stress values are to be evaluated in the rock mass near the pillar edges.

_{S}:H

]]>

From Figure 17, it is apparent that, for the present example, the average pillar stress will be overestimated if the seam modulus is assumed to be very stiff. (As the rock modulus: seam modulus tends to zero, the pillar effectively becomes rigid). For cases where the seam modulus is much lower than the host rock modulus, it is then recommended that the linear stiffness model should be used with the pillar stiffness estimated using Equation [18].

**Conclusions**

The computation of average pillar stress values is important when attempting to establish criteria for pillar design and in the analysis of the stability of tabular pillar layouts. One of the default 'classic' numerical methods that are used to determine pillar stresses is the displacement discontinuity method. In many instances it is not clearly understood that this approach does suffer from some limitations, particularly in relation to the fact that in coarse element simulations, the results can depend on the chosen element mesh size. The nature of this error is highlighted in the paper and some strategies are suggested to bound the magnitude of these errors. It is demonstrated as well that the popular linear stiffness approximation to pillar or seam compressibility does appear to allow reasonably accurate estimates to be made of the average pillar stress when either the pillar height is varied or when the seam modulus differs from the host rock modulus. At the same time, it is noted that although this approximation has some utility in allowing the determination of average pillar stress values, it may lead to significant errors when analysing the stress distribution off the reef plane or along the seam horizon in positions that are close to the edges of pillars or abutments. It is important to explore these effects in relation to pillar foundation failure mechanisms in future studies.

A practical implication of this study is that if the seam modulus is noticeably lower than that of the host rock, such as for coal seams, it is advisable to use a linear stiffness model for the pillars instead of a 'rigid' pillar assumption. This added complexity seems unnecessary, however, when simulating hard rock pillars in mines where the seam modulus is very similar to that of the surrounding rock.

]]>**Acknowledgements**

The authors wish to acknowledge the assistance of Dr Elizaveta Gordeliy in performing an independent check on the estimated value of the average stress in the circular pillar example used in the paper. We also acknowledge the helpful comments relating to tributary area theory made by a referee. In addition we wish to acknowledge our appreciation to Professor Anthony Peirce for his detailed review and crosschecking of the paper as well as comments on the Richardson extrapolation technique.

**References**

1. HEDLEY, D.G. F. and GRANT, F. Stope-and-pillar design for Elliot Lake Uranium Mines. *CIM Bulletin*, vol. 65, 1972. pp. 37-44. [ Links ]

2. OZBAY, M.U., RYDER, J.A., and JAGER, A.J. The design of pillar systems as practiced in shallow hard-rock tabular mines in South Africa, *Journal of the South African Institute of Mining and Metallurgy*, Jan/Feb 1995. pp. 7-18. [ Links ]

3. MARTIN, C.D. and MAYBEE, W.G. The strength of hard rock pillars, *International Journal of Rock Mechanics and Mining Science*, vol. 37, 2000. pp. 1239-1246. [ Links ]

4. MALAN, D.F. and NAPIER, J.A.L. The design of stable pillars in hard rock mines: A problem solved? *Journal of the Southern African Institute of Mining and Metallurgy*, 2011. [ Links ]

5. SALAMON, M.D.G. and ORAVECZ, K.I. Rock Mechanics in Coal Mining. Coal Mining Research Controlling Council, Chamber of Mines, Johannesburg, 1976. [ Links ]

6. OZBAY, M.U. The stability and design of yield pillars located at shallow and moderate depths. *Journal of the South African Institute of Mining and Metallurgy*, vol. 89. no. 3, 1989. pp. 73-79. [ Links ]

7. SALAMON, M.D.G. and MUNRO, A.H. A study of the strength of coal pillars, *Journal of the South African Institute of Mining and Metallurgy*, vol. 68, 1967. pp. 56-67. [ Links ]

8. VAN DER MERWE, J.N. New pillar strength formula for South African coal. *Journal of the South African Institute of Mining and Metallurgy*, 2003. pp. 281-292. [ Links ]

9. WATSON, B.P., RYDER, J.A., KATAKA, M.O., KUIJPERS, J.S., and LETEANE, F.P. Merensky pillar strength formulae based on back analysis of pillar failures at Impala Platinum. *Journal of the Southern African Institute of Mining and Metallurgy*, vol. 108, 2008. pp. 449-461. [ Links ]

10. RYDER, J.A. and JAGER, A.J. A textbook on rock mechanics for tabular hard rock mines. *Safety in Mines Research Advisory Committee*, Johannesburg, 2002. [ Links ]

11. SPENCER, D. A case study of a pillar system failure at shallow depth in a chrome mine. Hagan, T.O. (ed.). *Proceedings SARES99*, 2nd Southern African Rock Engineering Symposium, 1999. pp. 53-59. [ Links ]

12. HACKETT, P. An elastic analysis of rock movements caused by mining. *Transactions of the Institution of Mining Engineers*, 118, 1959. pp. 421-435. [ Links ]

13. BERRY, D.S. and SALES, T.W. An elastic treatment of ground movement due to mining-Part III. Three dimensional problem, transversely isotropic ground. *Journal of the Mechanics and Physics of Solids*, vol. 10, 1962. pp. 73-83. [ Links ]

14. SALAMON, M.D.G. Elastic analysis of displacements and stresses induced by the mining of seam or reef deposits-Part I: Fundamental principles and basic solutions as derived from idealised models. *Journal of the South African Institute of Mining and Metallurgy*, vol. 63, 1963. pp. 128-149. Part II: Practical methods of determining displacement, strain and stress components from a given mining geometry. *Journal of the South African Institute of Mining and Metallurgy*, vol. 64, 1964. pp. 197-218. [ Links ]

15. CROUCH, S.L. Computer simulation of mining in faulted ground. *Journal of the South African Institute of Mining and Metallurgy*, 1979. pp.159-173. [ Links ]

16. CROUCH, S.L. and STARFIELD, A.M. Boundary Element Methods in Solid Mechanics. George Allen & Unwin, London, 1983. [ Links ]

17. WILES, T.D. Reliability of numerical modeling predictions. *International International Journal of Rock Mechanics and Mining Science*, vol. 43, 2006. pp. 454-472. [ Links ]

18. RYDER, J.A. and NAPIER, J.A.L. Error analysis and design of a large-scale tabular mining stress analyser. *5th International Conference on Numerical Methods in Geomechanics*, Nagoya, Japan, 1985. pp. 1549-1555. [ Links ]

19. NAPIER, J.A.L. and MALAn, D.F. The computational analysis of shallow depth tabular mining problems, *Journal of the Southern African Institute of Mining and Metallurgy*, vol. 107, 2007. pp. 725-742. [ Links ]

20. SALAMON, M.D.G. Deformation of stratified rock masses: A laminated model, *Journal of the South African Institute of Mining and Metallurgy*, vol. 91, 1991, pp. 9-25. [ Links ]

21. HEASLEY, K.A. Numerical Modeling of Coal Mines with a Laminated Displacement-Discontinuity Code. Ph.D. Thesis. Colorado School of Mines, Golden, Colorado, 1998. [ Links ]

22. SALAMON, M.D.G. Two-dimensional treatment of problems arising from mining tabular deposits in isotropic or transversely isotropic ground. *Int. International Journal of Rock Mechanics and Mining Science*, vol. 5, 1968. pp. 159-185. [ Links ]

23. ANDREWS, G.E., ASKEY, R., and ROY, R. Special Functions. Cambridge University Press, 1999. [ Links ]

24. GORDELIY, E. Personal communication. March, 2011. [ Links ]

25. GORDELIY, E. and DETOURNAY, E. Displacement discontinuity method for modeling axisymmetric cracks in an elastic half-space. *International Journal of Solids and Structures*, 2011. [ Links ]

26. KOPAL, Z. Numerical Analysis. John Wiley & Sons, New York, Second Edition, 1961. [ Links ]

27. ORAVECZ, K.I. Loading of coal pillars in bord and pillar workings, Chamber of Mines of South Africa, Research Report 40/73, 1973. [ Links ]

Paper received Feb. 2011; revised paper received Jul. 2011

*© The Southern African Institute of Mining and Metallurgy, 2011. SA ISSN 0038-223X/3.00 + 0.00.*