NON-LINEAR FINITE-ELEMENT MODELLING OF ROOM AND PILLAR MINE WORKINGS INCLUDING THE STRAIN-SOFTENING BEHAVIOUR OF THE ROCK MASS SALIM BENSEHAMDI and ABDELBAKI SERIANI About the authors Salim Bensehamdi University Badji Mokhtar, Faculty of Earth Science, Mining Engineering Department Bp 12 (23000) Annaba, Algeria E-mail: seriani7@yahoo.fr Abdelbaki Seriani University Badji Mokhtar, Faculty of Earth Science, Mining Engineering Department Bp 12 (23000) Annaba, Algeria E-mail: seriani7@yahoo.fr Abstract A two-dimensional model adopting post-failure criteria was used to simulate the behaviour of the rock mass and the development of yield zones around room and pillar mine workings. The model conformed to the strain-softening behaviour of the rock mass and accounted for its post-failure residual strength. The structural-stability-analysis approach accounted for the main features of the mine structures' yield produced during loading through changes in the rock material's stiffness and the subsequent evolution of the stresses. A comprehensive parametric analysis was performed and the inevitable effect of the interaction of the roof, pillar and floor on the overall stability limit of the mine was investigated. The numerical results clearly showed that the finite-element linear models could not realistically represent the true behaviour of the mine structure. However, they clearly demonstrated the limitations of the finite-element linear solutions in representing the true behaviour of the mine structure, particularly when the rock-mass structure is relatively weak, and that a nonlinear approach was justified. Keywords non-linear FE analysis, yielding, plastic zones, room and pillar mining, residual strenght, stability assessment 1 INTRODUCTION One of the most difficult design problems in practical rock mechanics arises in conditions of the rock mass's complex non-linear constitutive behaviour, including structural discontinuities, and the non-homogeneity of the medium [7]. Field measurements and laboratory tests have shown the presence of the strain-weakening behaviour of the rock mass and have indicated that in many cases the assumption that the rock is linearly elastic leads to calculated stresses and displacements that disagree significantly with the measured values. In particular, as the rock mass around an excavation may exist in the post-yield state [2], [3], a realistic approach should incorporate the effect of the post-yield behaviour in the analysis [5], [10]. In the current finite-element analysis this has usually been achieved by adopting an elastic, perfectly plastic model and gradually changing the material coefficients of the yielding materials using a quasi-elastic finite-element solution [4]. In this paper an advanced elastic-plastic LUSAS finite-element program [6] was used to predict the distribution of stresses during the plastic and elastic strain states, and to simulate the possible mechanisms of the yield of the rock mass around a mine opening. This accounted for the residual strength after the failure. 2 DESCRIPTION OF THE MATERIAL MODEL The solutions of existing non-linear material models that are used to simulate rock mechanics and other geotech-nical problems are somewhat limited in their ability to properly reflect the behaviour of the complete stressstrain curve of most rock materials. [5] Both laboratory and field measurements have indicated that immediately after the peak strength level is reached, the stress-strain curve drops with a negative gradient and then flattens out at a residual strength. The classical solution adopting the constant-load incremental technique has been shown to fail when the solution reaches the limit points ACTA GeOTeCHNICfi SLOVeNICfi, 2008/2 5. S. BENSEHAMDI & A. SERIANI: NON-LINEAR FINITE-ELEMENT MODELLING OF ROOM AND PILLAR MINE WORKINGS INCLUDING on the material stress-strain curve, as shown in Fig. 1 and Fig. 2, [4], where either load- or displacement reversal occurs, resulting in a singular stiffness matrix, which automatically leads to the failure of the solution to converge. The LUSAS finite-element programs [6] used in this analysis contain a model that allows a simulation of the material's behaviour well into the strain-softening part of the stress-strain curve. Displacement, LI Figure 1. LUSAS modified arc length incrementation solution. A general method that may follow the solution path through the well-known post-peak portion of the material's stress-strain curve is the modified arc-length technique [1]. The salient characteristic of the arc-length technique is that the load level does not remain constant at each load increment, where during each iteration the load level is modified so that convergence near the limit points A and B can be achieved, Fig. 1. A typical, modelled, material stress-strain characteristic is illustrated in Fig. 2, where the complete stress-strain curve of the tested rock material is matched by a series of elastic and plastic lines, K0 , K1, and Kn , corresponding to the material's performance in the elastic and plastic states [5]. The initial portion of the model, as matched by the plastic line k1 , represents the work-hardening behaviour of the material, which is characterized by increasing stress with plastic deformation, while the second portion of the model shows the behaviour of the strain-softening state of the material, having a residual strength that decreases with plastic deformation, and finally an ideally plastic state of the material, where the deformation increases at constant stress. Figure 2. Modelled material stress-strain curve. 3 CONSTITUTIVE LAW In the elastic range, the way in which the stress and the strain are related for a material under load is described qualitatively by its constitutive behaviour. For an isotropic body undergoing plane-strain deformation, the stress-strain relation follows Hooke's law Eq. (1), that is: Me = [E]e Me Me = £x II £y ^z (1) [ El =■ (1-v) (1 + v)(1-2v) 1-v 1-v 1 0 1-2v 2(1-v) For non-linear material behaviour the plastic state is specified by: 1. A yield function to specify the onset of plastic deformation, Eq. (2): F({a} ,| K|)=0 (2) v 1 0 v 10. ACTA GEOTECHNICA SLOVENICA, 2008/2 S. BENSEHAMDI & A. SERIANI: NON-LINEAR FINITE-ELEMENT MODELLING OF ROOM AND PILLAR MINE WORKINGS INCLUDING where F: Yield function {o}: Stress vector |K|: Hardening, softening parameter In classical plasticity, stress states that provide a positive value of the yield function cannot exist. However, in numerical models, positive values of the yield function indicate that yielding should occur and the stress state is modified by accumulating plastic strains until the yield criterion is reduced to zero. This process is known as the plastic-corrector phase or return mapping. 2. A flow rule to define the plastic straining is given by Eq. (3): dF S (£>P = A (3) dM where 5{e|p: Increment vector of the plastic strain dF : Direction of the plastic strain 3{a} X: Lagrangian plastic multiplier defining the magnitude of the plastic strain 3. A hardening and softening rule to define the evolution of the yield surface with plastic strain. This is defined by describing the evolution of the yield function in relation to the effective plastic strain, eegp, using a series of straight-line segments, Fig.2. Kefip = o y +| Kn|.e, effp (4) where Keffp: Effective hardening, softening slope aY : Initial yield stress |Kn|: Hardening, softening slopes During an increment of stress, 6(a), changes in the strain are assumed to be the sum of an elastic and a plastic component. Thus, the concept of total strain will be invoked, Eq. (5): S{£} = S {£}e+S {£}p (5) The elastic strain components, 5(e)e , are related to the stress components by a matrix of constant [E], known as the stiffness matrix Eq. (6): S (e}e =[E]-1.S {£} (6) Taking account of the elastic and plastic components produces Eq. (7): S {e}e =[Ep.S {e} + A dF (7) The elastic-plastic stress and strain increments may be related by the following equation. Eq. (8): where S {£}e =[ E L .S {£} (8) [E] is the elasto-plastic stiffness matrix During the elastic-plastic analysis, the material's stiffness matrix is updated by the new, computed elastic-plastic stiffness matrix, [E] , at each increment of the finitely element solution [2], [3]. 4 THE NON-LINEAR FINITE-ELEMENT TECHNIQUE The technique for performing a non-linear finite- element analysis is illustrated in Fig. 3 (next page). This technique is summarized as follows: 1) First, load increments are applied to the mine structure and the strains and hence the stresses are found at the Gaussian points in the elements. For each increment of load, an initial material stiffness is used and the elastic solution is obtained. 2) During the solution the courses of all the finite elements are checked for yield. If the stresses at the Gaussian points lie within the previously prescribed yield surface, the stress update has been completed. Otherwise, stress lying outside the yield surface must be returned to the yield surface by plastic straining. During each iteration cycle, computed stresses and strains are added to the total already accumulated and a new material-stiffness matrix is reformulated for the next load increment. Within each load increment the system of Eq. (9): S M,=[E1 .S{£}, (9) must be solved for the strain increment ¿{e} where S{a|i : Increment of stress during iteration, i [E]iep : Updated elasto-plastic stiffness matrix ¿{e} : Strain increment during iteration 5 {£}, =(<5 {£}e +< {£} p ) 3) Steps (1) and (2) are repeated for all the increments of load that constitute the total load applied to the structure. 10. ACTA GEOTECHNICA SLOVENICA, 2008/2 S. BENSEHAMDI & A. SERIANI: NON-LINEAR FINITE-ELEMENT MODELLING OF ROOM AND PILLAR MINE WORKINGS INCLUDING (Start ] Read Data Compute Elastic Stiffness Matrix Read in Applied Load Increment I Has Failure Criterion Been Exceeded ? No Solve Elastic Problem Yes - Compute Elastic Strain Increments - Compute Stress Increments and Add to Total Stresses left from the Last Load Increment < Has Convergence Been Achieved ? Yes 1 Update Stiffness Matrix r 1JÎ4E+Î> (b) (Pa) ])!)L>9Ï ■c ;:îîl.; I411E+4? éUÎWJ SflJltmP -«. JîtTE'ÎF -0 1 BB^I > D !■ V. » il iliiEOr1 j.ï.jnE'ir & H H 3BE H ¿T E Hîlt+tÎ Figure 8. Shear-stress contours in a yielding plastic FE solution (a), and unyielded elastic FE solution (b), of a room and pillar mine structure at the final equilibrium steady state. 10. ACTA GEOTECHNICA SLOVENICA, 2008/2 S. BENSEHAMDI & A. SERIANI: NON-LINEAR FINITE-ELEMENT MODELLING OF ROOM AND PILLAR MINE WORKINGS INCLUDING A relationship exists between the yield deformations of the pillar, where three distinct stress zones develop as the pillar is yielding: - A zone of local yielding in a stress-relieved area, - A zone of transition from a yielding to the solid state in the pillar, where the stress -concentration area was developed, - A zone of pillar core that has a uniform confined stress. The overall pillar stability depends upon the geometrical development of these principal zones. For example, an increase in the yielding pillar area results in a decrease of the solid core, and hence a decrease in its bearing performance, which automatically leads to the instability of the system. Fig. 9 shows the yielding state of the mine structure where the growth of the plastic zone for each load increment is given by the yield Gaussian points. At 50% of the maximum applied load, Fig. 9A, the first plastic zones are shown to occur at the immediate roof-pillar intersection, over the cross-section of the pillar and, to a limited extent, at the floor-pillar intersection. In this case the solid pillar area is over 50% of the total Û « ÏÎ Y / \t i Yielded Gauss Points 1 < < _ - - - A- 1st Load increment, ov= 50% Figure 9. Spread of the yield zone at A A m 1ft Y ill! w n. V 'i 0 p 4i P \V 9 . J "l — B- 2nd Load Increment, av= 70% \ Y \ ft? PIP inr m IRr if If? r ii IMI im V - / ? 5 K3 / j • ft Mi * . rt ■n n n S S ifN --- C- 4th Load Increment, ov= 100% applied load levels; the extraction ratio r = 60%. 10. ACTA GEOTECHNICA SLOVENICA, 2008/2 S. BENSEHAMDI & A. SERIANI: NON-LINEAR FINITE-ELEMENT MODELLING OF ROOM AND PILLAR MINE WORKINGS INCLUDING area, and the pillar is assumed to be stable. As the load is increased, the plastic zone spreads upwards into the immediate roof and laterally into the pillar, Fig. 9B. The effective pillar core is reduced to less than half of the pillar area, but can still provide support. At the final stage of loading, the yielded zone reaches its maximum as the mine structure passes to its final equilibrium steady state, as shown in Fig. 9C, where the plastic zone is seen to increase considerably in both the roof and the pillar, while the floor remains almost unaffected. The effective pillar area is reduced to more than half of the total pillar area and overall pillar instability is likely. At low extraction ratios of 0.25 and 0.40, as shown in Fig. 10 (A, B), the developed yield zone is confined to an area around the rock-mass opening interface. In both cases the area of effective support is shown to be larger than the outer yield zone. Under these conditions, the pillar is almost at its maximum capability and the mine structure is stable. The non-linear finite-element analysis was further extended to analyse the influence of the roof, pillar and floor interactions on the overall pillar stability. In order A- Kxt. 25% C- F.vt. |iat)0= ft] ft Figure 10. Spread of the yield zone at the full applied load level and various extraction ratios, r. 10. ACTA GEOTECHNICA SLOVENICA, 2008/2 S. BENSEHAMDI & A. SERIANI: NON-LINEAR FINITE-ELEMENT MODELLING OF ROOM AND PILLAR MINE WORKINGS INCLUDING to find out the effect of this interaction as well as the effect of the change of their properties on the overall pillar stability, three models were analysed using different properties for the roof, pillar and floor. Fig. 11 shows the predicted vertical and horizontal stresses of the pillar for model A, having a stiff roof, a weak pillar and a weak floor; and model B, having a weak roof, a weak pillar and a stiff floor; and model C, where the roof and the floor are taken to be stiff and the pillar is kept weak. It is clear from Fig. 12 ( next page) that for the same pillar material in the three suggested models and the different combinations of material properties for the roof and floor, model C shows complete yielding of the pillar, while in model B the area of effective support is considerably reduced, but to a lesser extent than with model C. Finally, in model A, a moderate effective area is indicated in the pillar, which yields to a certain degree. It can be concluded from these findings that the interaction of the roof, pillar and floor has a significant effect on the pillar's overall stability. Therefore, when assessing the stability of a given room and the pillar mine workings, the roof, pillar and floor should be considered as an integrated structure. This approach generally yields more accurate results. Figure 11. Vertical and horizontal pillar stresses, determined by the elastic-plastic finite-element solution for different models. 10. ACTA GEOTECHNICA SLOVENICA, 2008/2 S. BENSEHAMDI & A. SERIANI: NON-LINEAR FINITE-ELEMENT MODELLING OF ROOM AND PILLAR MINE WORKINGS INCLUDING IR)' ■y / / t ; / — — á | St T ft r r f Model A N X f) ft ft? m ft ft / ri q \ J n ffi ft -X, t i \ J LI — . I Oi? [í m n mmmsm TT •W' -Í- as — ■if m ifife 1& (SÍ m i? i r t n Model B * / Model C Figure 12. Spread of the yield zone in the mine structure for the three models. 7 CONCLUSION Based on the modelling results presented above, the following conclusions may be drawn: - The strain-weakening behaviour of rock layers around a mine opening was simulated using an advanced non-linear elastic-plastic finite-element model, which accounts for the residual strength after failure. This simulation was achieved using a step-by-step iterative computational procedure in which the rock material's stiffness was updated after each run according to the initially applied load, until the solution reached the final equilibrium steady state of the mine structure. The numerical results clearly demonstrate the limitations of the linear model when it comes to realistically representing the overall structure behaviour Fig. 5, particularly when the structure of the rock mass is relatively weak, and that a non-linear approach was justified. - As progressive mining was simulated in the model, the maximum load in the pillar increased and trans- 10. ACTA GEOTECHNICA SLOVENICA, 2008/2 S. BENSEHAMDI & A. SERIANI: NON-LINEAR FINITE-ELEMENT MODELLING OF ROOM AND PILLAR MINE WORKINGS INCLUDING ferred laterally towards the centre of the pillar as yielding occurred, this trend of lateral stress transfer and load build up in the core of the pillar is consistent with that observed in situ. - Based on the simulation results, the upper and lower bounds of the overall room and pillar stability were obtained. Parameters such as the roof, the pillar width, wp , the extraction height, hp , the variability of the pillar, roof and floor strength have been shown to significantly affect the mechanism of roof, pillar and floor yield. references [1] Riks, S.F. (1979). An incremental approach to the solution of snapping and buckling problems. International Journal of Solid Structure 15, 29-551. [2] Nayak, G. C. and Zienkiewicz, O. C. (1972). Elasto-plastic stress analysis. A Generalization for various constitutive relations including strain-softening. Int. J. Num Meth. Eng., 5, 113-135. [3] Zienkiewicz, O. C. Valliappan, S. and King, (1969). Elasto-plastic solutions of engineering problems 'Initial stress' finite element approach. Int. J. Num. Meth. Eng.. 1. 75-100. [4] Crisfield, M. A. (1981). A fast incremental/iterative solution procedure that handles snap-through. Computers and Structures, 13, 55-62 [5] Cheng, G. and Karmis, M. (1988). Computer modelling of yield pillar behaviour using post-failure criteria. Proceedings of '88 7th Int. Conf. Ground Control in Mining. West Virginia University, Morgantown, USA. [6] LUSAS, (1989). Advanced nonlinear finite element analysis.FEA.Ltd, Forge House High street, Kingston upon Thames, Surrey, KTN 1HN. [7] Park, D. W. and Ash, N. F. (1985). Stability analysis of entries in deep coal mines using finite element method. Mining Science and Technologie, 3, 11-20. [8] Salamon , M.D.G., Badr, S., Mandoza, R., Ozbay, M. U. (2003). Pillar failure in deep coal seems: Numerical simulation. ISRM, 2003 - technology road map for rock mechanics. South Africa. [9] Colwell, M. G., (1998). Chain pillar design (calibration of ALPS), ACARP research project C6063 Australian coal research Ltd. 66p. [10] Mark, C. and Barton T. (1996). The uniaxial compressive strength of coal: should it be used to design pillar? Proceedings of the 14th conference on ground control in mining, West Virginia University, August 1-3, 63-71. 10. ACTA GEOTECHNICA SLOVENICA, 2008/2