Modelling of Longwall Mining-Induced Strata Permeability Change

Abstract

The field measurement of permeability within the strata affected by mining is a challenging and expensive task, thus such tests may not be carried out in large numbers to cover all the overburden strata and coal seams being affected by mining. However, numerical modelling in conjunction with a limited number of targeted field measurements can be used efficiently in assessing the impact of mining on a regional scale. This paper presents the results of underground packer testing undertaken at a mine site in New South Wales in Australia and numerical simulations conducted to assess the mining-induced strata permeability change. The underground packer test results indicated that the drivage of main headings (roadways) had induced a significant change in permeability into the solid coal barrier. Permeability increased by more than 50 times at a distance of 11.2–11.5 m from the roadway rib into the solid coal barrier. The tests conducted in the roof strata above the longwall goaf indicated more than 1,000-fold increase in permeability. The measured permeability values varied widely and strangely on a number of occasions; for example the test conducted from the main headings at the 8.2–8.5 m test section in the solid coal barrier showed a decline in permeability value as compared to that at the 11.2–11.5 m section contrary to the expectations. It is envisaged that a number of factors during the tests might have had affected the measured values of permeability: (a) swelling and smearing of the borehole, possibly lowering the permeability values; (b) packer bypass by larger fractures; (c) test section lying in small but intact (without fractures) rock segment, possibly resulting in lower permeability values; and (d) test section lying right at the extensive fractures, possibly measuring higher permeability values. Once the anomalous measurement data were discarded, the numerical model results could be seen to match the remaining field permeability measurement data reasonably well.

Appendix: Layered Cosserat Model

For simplicity and clarity of presentation only a two-dimensional plane strain model formulation will be described in this report although the formulation is fully implemented in three dimensions. Using the Cartesian coordinates (x ₁, x ₂), the material point displacement can be defined by a translational vector (u ₁, u ₂) and by a rotation Ω₃. Here, axis x ₃ is aligned in the out of plane direction and axis x ₂ is perpendicular to the layers.

The two-dimensional Cosserat model has four non-symmetric stress components σ ₁₁, σ ₂₂, σ ₂₁, σ ₁₂ and two couple stresses m ₃₁, m ₃₂. When the rock layers are aligned in the x ₁-coordinate direction, the moment stress term m ₃₂ vanishes. The four stresses are conjugate to four deformation measures γ ₁₁, γ ₂₂, γ ₂₁, γ ₁₂ defined by:

$$\gamma_{ij} = \frac{{\partial u_{j} }}{{\partial x_{i} }} - \varepsilon_{3ij} \Omega_{3}$$

(A1)

where ɛ _3ij is the permutation tensor, and the couple stress m ₃₁ is conjugate to the respective curvature κ ₁ defined by:

$$\kappa_{1} = \frac{{\partial \Omega_{3} }}{{\partial x_{1} }}$$

(A2)

The elastic stress strain relationships are described by:

$$\varvec{\sigma} = \left[ {D_{e} } \right] \, \varvec{e}_{e}$$

(A3)

where $\varvec{\sigma} = \left\{ {\sigma_{11} ,\sigma_{22} ,\sigma_{21} ,\sigma_{12} ,{\text{m}}_{31} } \right\}$, $\varvec{e} = \left\{ {\gamma_{11} ,\gamma_{22} ,\gamma_{21} ,\gamma_{12} ,\kappa_{ 1} } \right\}$ and (A4).

$$D \, = \, \left[ {\begin{array}{*{20}c} {A_{11} } & {A_{12} } & 0 & 0 & 0 \\ {} & {A_{22} } & 0 & 0 & 0 \\ {} & {} & {G_{11} } & {G_{12} } & 0 \\ {} & {\text{symm}} & {} & {G_{22} } & 0 \\ {} & {} & {} & {} & {B_{1} } \\ \end{array} } \right]$$

(A5)

here,

$$A_{11} = \frac{E}{{1 - \nu^{2} - \frac{{\nu^{2} \left( {1 + \nu^{2} } \right)}}{{1 - \nu^{2} + \frac{E}{{hk_{\text{n}} }}}}}}$$

(A6)

$${\text{A}}_{22} = \frac{1}{{\frac{{1 - \nu - \nu^{2} }}{{{\it{E}}\left( {1 - \nu } \right)}} + \frac{1}{{{\it{hk}}_{\text{n}} }}}}$$

(A7)

$$A_{12} = A_{21} = \frac{\nu }{1 - \nu }A_{11}$$

(A8)

$$\frac{1}{{G^{\prime}}} = \frac{1}{G} + \frac{1}{{bk_{\text{s}} }}$$

(A9)

$$B = \frac{{D_{f} }}{b}\left( {1 - \frac{{G^{\prime}}}{G}} \right)$$

(A10)

where E is the Young’s modulus of the intact layer, ν is the Poisson’s ratio, h is the layer thickness, G is the shear modulus of the intact layer, k _n and k _s are the joint normal and shear stiffnesses.

The layer interfaces can exhibit three different modes of behaviour: (a) elastically connected with the interface normal and shear stiffness, (b) plastic with frictional sliding and (c) disconnected with tensile opening. Similarly the rock layer may either deform elastically or may sustain some plastic deformation as well. With this in mind, the rate of the deformation tensor is decomposed into elastic and plastic parts (Adhikary and Guo 2002):

$$\dot{e} = \dot{e}_{\text{e}} + \dot{e}_{\text{p}}$$

(A11)

In a manner similar to the conventional plasticity theory the rate of plastic deformation is assumed to be equal to:

$$\dot{e}_{\text{p}} = \dot{\lambda }\frac{\partial g}{{\partial \dot{\sigma }}}$$

(A12)

where $\dot{\lambda }$ is the so-called plastic multiplier and g is the plastic potential function. Then the incremental elasto-plastic relationships in the general form can be expressed as usual:

$$\dot{\sigma } = \left[ {D_{\text{ep}} } \right] \, \dot{e}$$

(A13)

where, $\dot{\sigma }$ and $\dot{e}$ are the incremental stress and strain, and

$$D_{\text{ep}} = D_{\text{e}} - \alpha \frac{{D_{\text{e}} \left\{ {\frac{\partial g}{\partial \sigma }} \right\}\left\{ {\frac{\partial f}{\partial \sigma }} \right\}^{\text{T}} D_{\text{e}} }}{{\left\{ {\frac{\partial f}{\partial \sigma }} \right\}^{\text{T}} D_{\text{e}} \left\{ {\frac{\partial g}{\partial \sigma }} \right\}}}$$

(A14)

Here f ≤ 0 is the yield function, g is the plastic potential and α is defined as:

$$\alpha = 1{\text{ if }}f = 0 \quad {\text{ and }}\quad \dot{\lambda } \, ({\text{plastic multiplier}}) > 0; \, 0{\text{ if }}f < 0 \, {\text{ and}}\, /{\text{or }}\dot{\lambda } \le 0$$

(A15)

The course of derivation leading to Eq. A14 is exactly the same as in standard continua so it will not be discussed in detail here. The yield and plastic potential functions adopted in this study will be introduced.

The yield criteria for interface sliding, $f_{\text{s}}^{\text{joint}}$, and the corresponding plastic potential function, $g_{\text{s}}^{\text{joint}}$, for a joint parallel to the 1-axis are defined as (here tensile stresses are assumed to be positive):

$$f_{\text{s}}^{\text{joint}} = \left| {\sigma_{21} } \right| + \sigma_{22} \tan \varphi^{\text{joint}} - c^{\text{joint}} = 0, \, g_{\text{s}}^{\text{joint}} = \left| {\sigma_{21} } \right| + \sigma_{22} \tan \psi^{\text{joint}}$$

(A16)

where Φ^joint, Ψ^joint and c ^joint designate the angle of friction, dilation angle and the cohesion of the joints respectively. Similarly, the yield criterion for the tensile opening and the corresponding plastic potential function are written as:

$$f_{\text{t}}^{\text{joint}} = \sigma_{22} - \sigma_{\text{ten}}^{\text{joint}} = 0, \, g_{\text{t}}^{\text{joint}} = \sigma_{22}$$

(A17)

where $\sigma_{\text{ten}}^{\text{joint}}$ is the joint tensile strength.

The Cosserat theory incorporates asymmetric stresses and couple stresses in the rock layer. For simplicity, the rock yield function is formulated on the basis of symmetric part of the stress tensor with the incorporation of the moment stresses. In this study, the couple stress is assumed to introduce additional axial stress in the rock layer. In a beam (rock layer) subjected to a bending moment, axial stress and strain vary linearly across the depth of the section. As the bending moment is increased the yield stress is attained first at the outer fibres. Once yield stress is attained the rock layer could be considered broken and subsequently a zero tensile strength could be assigned. The magnitude of this bending moment in a beam can be calculated in terms of ultimate yield stress $\sigma_{\text{ten}}^{\text{rock}}$ in the following manner:

$$M^{\text{yield}} = \sigma_{\text{ten}}^{\text{rock}} bh^{2} /6$$

(A18)

Here, b is the breadth of the beam in the out of plane direction x ₃ and h is the beam thickness.

The plastic couple stress can be expressed as the plastic moment per unit area as:

$$m_{31}^{\text{yield}} = M^{\text{yield}} /bh$$

(A19)

Equations (3.16) and (3.17) yields:

$$\sigma_{\text{ten}}^{\text{rock}} = 6m_{31}^{\text{yield}} /h$$

(A20)

Similar to the joint failure modes, the rock matrix is assumed to fail either in tension or shear. Tensile strength of rocks is often observed to be 10–40 times less than the uniaxial compressive strength. Thus, when there is relatively large moment stress in the rock layer due to bending, it is most likely that the rock layers will fail in tension. When the moment stress in the rock layer is small (in the absence of layer bending) both shear and tensile failure remain a possibility and will largely be determined by the conventional stresses. In this case asymmetry in stresses will be small. Thus, in order to capture the tensile failure of the rock layer subjected to bending, the effective normal stress in the rock layer is defined by:

$$\sigma_{11}^{N} = \frac{{6\left| {m_{31} } \right|}}{h} + \sigma_{11}$$

(A21)

where, σ ₁₁ is the conventional normal stress acting along the direction of layering.

Then the rock layer yield criterion can be defined in terms of the symmetric part of the stress tensors as:

$$f_{\text{s}}^{\text{rock}} = (\sigma_{11}^{N} - \sigma_{22} )^{2} + (\sigma_{12} + \sigma_{21} )^{2} - \left( {\left( {\sigma_{11}^{N} + \sigma_{22} } \right)\sin \varphi_{r} - 2C_{r} \cos \varphi_{r} } \right)^{2} = 0$$

(A22)

and

$$f_{t}^{\text{rock}} = \sigma_{\hbox{max} } - \sigma_{ten}^{rock} = 0, \, \left( {{\text{here tension is }} + {\text{ve}}} \right)$$

(A23)

Similar to the joint plastic potential functions, the rock plastic potential functions are obtained simply by replacing rock friction angle by rock dilation angle. While formulating D _ep, B ₁ is either made equal to zero if rock layer is yielding or made equal to Eh ²/12(1 − ν ²) if joint alone is yielding. The finite element formulation of 2D Cosserat model is fully described in Mühlhaus (1993) and Adhikary and Dyskin (1997). Of course, when the bending moment vanishes (i.e. h = 0) the Classical continuum model is recovered.