A visco-elastic theoretical model for analysis of dynamic ground subsidence due to deep underground mining

Applied Mathematical Modelling xxx (2015) xxx–xxx

Contents lists available at ScienceDirect

Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm

A visco-elastic theoretical model for analysis of dynamic ground subsidence due to deep underground mining Wen-Xiu Li a,b,⇑, Chong-Yang Gao a, Xia Yin b, Ji-Fei Li c, Dong- Liang Qi a, Jing-Chao Ren a a

College of Civil Engineering and Architecture, Hebei University, Baoding 071002, PR China Hainan Technology and Business College, Haikou 570203, PR China c Research Institute of Geotechnical Engineering, Hebei University, Baoding 071002, PR China b

a r t i c l e

i n f o

Article history: Received 2 March 2013 Received in revised form 25 December 2014 Accepted 5 January 2015 Available online xxxx Keywords: Rock rheology theory Underground mining Dynamic ground subsidence Visco-elastic model Pillarless sublevel caving method

a b s t r a c t The dynamic ground subsidence due to mining is a complicated time-dependent process. Based on the theory of rock rheology, a theoretical model for the prediction and analysis of the dynamic subsidence due to deep underground mining is developed. The model is used for Chengjiao mine, an underground iron ore mine using pillarless sublevel caving method, to predict and analyse the subsidence in West and North mining areas. The theoretical results obtained were compared with actual subsidence data from West and North mining areas in Chengjiao Iron mine, eastern China. The agreement of the theoretical results with the filed measurements shows that the model is satisfactory and the formulae obtained are valid and thus can be effectively used for predicting the ground movement due to deep underground mining by pillarless sublevel caving method. Ó 2015 Elsevier Inc. All rights reserved.

1. Introduction Since the 1990s, there has been a growing demand for iron ore and steel products worldwide. With the development of world economy, the demand for mineral keeps increasing. However, shallow mineral resources are dwindling and iron mines in world have to get used to deep mining. Both observed and theoretical analysis show that the surface subsidence is continual and slow with a long duration of movement due to deep underground mining of iron deposits by the pillarless sublevel caving method. As the result of the engineering development, the depth at which iron deposits are excavated has increased substantially reaching up to 1.2 km in China and 3 km in Canada. The underground mining at such great depths could induce geostress redistribution and rotation to the area surrounding mine workings. The change in geostress regimes (especially the high horizontal tectonic stress), in turn, leads to the large deformation for the occurrence of the rockmass displacement (surrounding rock failure), such as the stability of underground roadway, ground surface subsidence. Many of mining practice proved that the range of overlying strata and ground surface movements (subsidence and deformation) have creep characteristics. It is the main factor of dynamic ground subsidence [1–10]. In addition, exploitation of minerals using caving methods such as the pillarless sublevel caving method will result in surface subsidence. For example, Chengjiao iron mine, in which the pillarless sublevel caving was used, there has been widespread deformation of land surface

⇑ Corresponding author at: College of Civil Engineering and Architecture, Hebei University, Baoding 071002, PR China. Tel.: +86 312 5079493; fax: +86 312 5079375. E-mail address: [email protected] (W.-X. Li). http://dx.doi.org/10.1016/j.apm.2015.01.003 0307-904X/Ó 2015 Elsevier Inc. All rights reserved.

Please cite this article in press as: W.-X. Li et al., A visco-elastic theoretical model for analysis of dynamic ground subsidence due to deep underground mining, Appl. Math. Modell. (2015), http://dx.doi.org/10.1016/j.apm.2015.01.003

2

W.-X. Li et al. / Applied Mathematical Modelling xxx (2015) xxx–xxx

Nomenclature in situ stress (MPa, kN/m2) c unit weight (kN/m3) H mining depth (i.e. overburden thickness, m) x ground surface coordinate in horizontal direction (dependent on the coordinate system xoz) t time (year) Wb (x, t) subsidence of visco-elastic beam at point x at the moment t (m) Wp (x, t) subsidence of visco-elastic foundation at the point x at the moment t (m) W (x, t) dynamic ground subsidence due to underground deep mining (m) W (0, t) the roof deflection (subsidence) before filling in mining boundary (m) U (x, t) dynamic ground horizontal displacement due to deep mining (m) r stress (MPa) e(x, t) strain on the line parallel to the neutral axis (horizontal deformation of top layer of the beam (mm/m) e_ ¼ @@te rate of change of e z distance from point x to the neutral axis (m) q(x, t) radius of curvature of the neutral axis (m) Eb elasticity modulus of the bending beam (GPa) Ep elastic modulus of the visco-elastic foundation (iron seam) (GPa) elastic modulus of filling body (GPa) Ek gb viscosity coefficient of the bending beam (Pa s) gp viscosity coefficient of visco-elastic foundation (Pa s) gk viscosity coefficient of filling body (Pa s) I moment of inertia of cross-section (m4) A the cross-section area of the beam (m2) M(x, t) moment of cross section of the beam (kN m) q(x,t) external load of the rock beam (MPa) u(x,t) foundation reaction force (MPa) w(x,t) reaction force of filling body for the visco-elastic beam (MPa) m mining thickness (m) S mining width of the sub-critical mining (m) L ‘‘half-wavelength’’ of stress wave and is a constant to be determined (m) l the width of gob (m) Wmax maximum surface subsidence (m) n subsidence coefficient and is a constant to be determined (non-dimensional parameter) L0 mining width of the critical mining (m) R radius of influence of ground subsidence (m) d the influence angle of ground subsidence (°) m rate of face advance (i.e. mining rate; m/year) Pz

far away from the mining area. This phenomenon can cause environmental problems and damage to surface and subsurface structures. In order to protect the environment and structures from these damages, precise ground movement prediction is essential. On the other hand, ground subsidence is a time dependent dynamic process. The development of ground movement, from the immediate roof to the surface, has a dynamic character and it is related to the progress of underground mining and time. In other words, ground surface subsidence due to deep mining is a dynamic process (time-dependent and rate-dependent) which obeys rock rheological principles. Ground movement prediction methods which can pre-calculate the final and intermediate stages of this process are important in mine design [1,3,4]. Deep rock masses show particular rheological properties, because of their complex geophysical environment and stress fields. Ground subsidence and deformation due to the underground deep mining, are complex processes affected by the mining method and velocity, timing of the mining process, rheological properties of overlying strata [6]. Therefore, it is necessary to study the dynamic ground subsidence due to deep mining. In fact, ground movement is a major problem associated with deep underground mining of iron ore body, causing damage to subsurface structures and surface buildings. To avoid adverse impacts of ground subsidence, a reliable prediction is essential. The prediction and analysis of ground movement due to underground mining have been studied by many scholars in this field and valuable results have been obtained [1–36]. Many scientists from the UK, USA, China, Australia and other countries attempted to predict the ground subsidence of mines using numerous methods, such as the profile function methods [12– 14,36], the influence function methods [15–20,36] and the void diffusion method [21]. Also, many theoretical studies are carried out using stochastic [22], elastic [23], and visco-elastic methods [6,8–10,24–26]. Recent efforts include the use of finite-element method [27,28], fuzzy system method [29], boundary element method [30], distinct element method [31],

Please cite this article in press as: W.-X. Li et al., A visco-elastic theoretical model for analysis of dynamic ground subsidence due to deep underground mining, Appl. Math. Modell. (2015), http://dx.doi.org/10.1016/j.apm.2015.01.003

W.-X. Li et al. / Applied Mathematical Modelling xxx (2015) xxx–xxx

3

FLAC method [32], and artificial neural network method [33] in the prediction of subsidence. In China, numerical modelling for the prediction of subsidence caused by underground mining has been carried out [34,35]. However, most prediction methods developed to analyse ground surface movement caused by underground mining have, in general, concentrated on the final state of displacements. Valuable investigation was taken up by Li, Guo and Hou [6]. In their paper, the viscoelasticity was applied to describe ground subsidence due to underground mining by pillarless sublevel caving method. However, this method appeared to be not suitable for the deep mining areas in China [37]. And marked progress has not been achieved in describing the dynamic ground subsidence due to deep mining. In the present work, the authors attempt to predict dynamic ground subsidence using a visco-elastic model (VEM). Detailed analyses were carried out to predict subsidence caused by underground mining of iron deposits. The validation of the model has been performed by predicting the subsidence profiles for the Chengjiao iron mine in China. We have found a reasonably good match between the observed and predicted dynamic ground subsidence curves. 2. The visco-elastic theoretical models 2.1. Basic assumptions With the theory of rheology, the nature of the process of dynamic ground subsidence due to deep mining can be revealed more effectively in time and space domains. This paper is mainly concerned with the application of viscoelasticity to dynamic surface subsidence. In this paper, to study the creep properties of argillaceous siltstone, creep tests were carried out on the intact rock samples for argillaceous siltstone from three drill holes of Chengjiao mine [37]. Creep is the time-dependent strain or deformation under constant axial stress. Kelvin visco-elastic creep model is selected to represent the creep behaviour based on the creep test results of weak rock mass [6,24,37]. According to the research results of Liu [24], Li et al. [6] and Hebei University [37], shows that the deformation characteristics of soft rock can be described by the Kelven model. In fact, the Kelvin rheological model has been widely used to describe viscoelastic deformation caused by underground excavation. In order to predict the dynamic surface subsidence due to deep mining, we will use Kelven model to describe dynamic ground movement due to deep mining by the pillarless sublevel caving method in soft rock strata. The analysis will be based on the following assumptions: (1) The laminated overburden and iron seam are horizontal (Fig. 1). (2) The roof and seam are consistent with Kelvin rheological model (Fig. 2). The constitutive equation is as follows [6,24,37]:

r ¼ E  e þ g  e_ ¼ E  e þ g 

@e ; @t

ð1Þ

where, r, e = stress and strain, respectively, e_ ¼ @@te = rate of change of e, and E, g = rheological parameters of the medium. (3) The roof and the iron seam are viscoelastic media, put forward the idea that the principle of ground movement is similar to that of roofs. (4) The in-situ stress at a given depth is

Pz ¼ cH;

ð2Þ

where, Pz is the vertical stress in MPa and H is the depth in metres, c is the unit weight. (5) The vertical strain of the foundation at any point is proportional to the deflection of the beam at that point. The properties of a rock mass are commonly represented by physical models consisting of elastic and viscous elements. A two-parameter viscoelastic rheological model of the Kelvin type is used in the modelling of the system as follows. Fig. 2 shows the two element model which is suitable for determining properties of a rock mass. It is now possible to calculate the dynamic ground subsidence due to deep mining by pillarless sublevel caving method.

P

z

X O

roof seam

Z Fig. 1. The horizontal overburden and iron seam.

Please cite this article in press as: W.-X. Li et al., A visco-elastic theoretical model for analysis of dynamic ground subsidence due to deep underground mining, Appl. Math. Modell. (2015), http://dx.doi.org/10.1016/j.apm.2015.01.003

4

W.-X. Li et al. / Applied Mathematical Modelling xxx (2015) xxx–xxx

Fig. 2. Kelvin viscoelastic model.

2.2. Basic theoretical models The strength of a rockmass changes with time, therefore, dynamic ground movement in an iron mine is also time dependent. According to the existing research results, we can build a new mathematical model for prediction of dynamic surface subsidence due to deep mining by the pillarless sublevel caving method. Based on the engineering geological conditions of Chengjiao mine in Shandong Province, the iron seam and laminated overburden are horizontal (Fig. 3). According to the creep experimental data and results of theoretical analysis [6,24], the roof and iron seam are consistent with Kelvin rheological model, which belongs to the visco-elastic rockmass. The seam is considered as a visco-elastic foundation and the roof as visco-elastic cylindrical bending plate on the foundation. Thus the unit width beam (roof) can be used to characterise the cylindrical bending plate (roof) [24]. Assuming that there is no gap between beam and foundation, i.e. the deflection of the beam is equal to the foundation subsidence at any point and any moment of time. Then we have:

W b ðx; tÞ ¼ W p ðx; tÞ ¼ Wðx; tÞ;

ð3Þ

where, Wb(x,t) is subsidence of visco-elastic beam at point x at the moment t, Wp(x,t) is subsidence of visco-elastic foundation at the point x at the moment t. The horizontal deformation of top layer of the beam meet the following condition [6,24,38]:

eðx; tÞ ¼ z=qðx; tÞ:

ð4Þ

According to the constitutive Eq. (1), we can get:

rðx; tÞ ¼ Eb 

  z @ 1 ; þ gb  z  @t qðx; tÞ qðx; tÞ

ð5Þ

where, q(x,t) is the curvature radius of the neutral axis, Eb is elasticity modulus of the cylindrical bending beam, gb is viscosity coefficient of the beam, z is the distance from point x to the neutral axis. The moment of cross section of the beam can be expressed as the following equation.

Mðx; tÞ ¼

Z

rðx; tÞ  z  dA ¼ I  A



  Eb @ 1 ; þ gb  @t qðx; tÞ qðx; tÞ

Pz

ð6Þ

x

o Overlying strata 1

Sublevel 2

Surrounding rock

n

z Fig. 3. The pillarless sublevel caving method.

Please cite this article in press as: W.-X. Li et al., A visco-elastic theoretical model for analysis of dynamic ground subsidence due to deep underground mining, Appl. Math. Modell. (2015), http://dx.doi.org/10.1016/j.apm.2015.01.003

5

W.-X. Li et al. / Applied Mathematical Modelling xxx (2015) xxx–xxx

R where, I ¼ A z2  dA is the moment of inertia of cross-section, A is the cross-section area of the beam. According to the principle of differential geometry, the curvature of the beam can be approximately represented by the following equation

1 ¼ qðx; tÞ 

@ 2 Wðx;tÞ @x2

1þ



@Wðx;tÞ @x

2 32



@ 2 Wðx; tÞ : @x2

ð7Þ

Thus, substituting (7) into (6), we can derive the following equation

Mðx; tÞ ¼ Igb 

@ 3 Wðx; tÞ @ 2 Wðx; tÞ þ IE  : b @x2 @t @x2

ð8Þ

According to the theory of mechanics of materials, the load q(x,t) of the beam can be written in the following form

qðx; tÞ ¼

@ 2 Mðx; tÞ @ 5 Wðx; tÞ @ 4 Wðx; tÞ ¼ Igb  þ IEb  : 2 4 @x @x @t @x4

ð9Þ

From Fig. 4, we know that the external load of the beam is composed of two parts: (1) The upper surface of the beam by effect of the overlying strata weight and the load is uniformly distributed: Pz = cH; (2) The lower surface of beam will bear the foundation reaction force u(x,t), the external load of the beam can be written in the following form:

qðx; tÞ ¼ cH  uðx; tÞ;

ð10Þ

where, q(x,t) is the external load of the rock beam in MPa and u(x,t) is the foundation reaction force in MPa. According to the assumption (5), the reaction forces of the foundation are proportional at every point to the deflection of the beam at that point. From Eq. (1) and reference [24] , the foundation reaction force can be written as:

uðx; tÞ ¼ Ep ep þ gp e_ p ¼

gp @ Ep Wðx; tÞ þ  Wðx; tÞ; m m @t

ð11Þ

where, Ep is elastic modulus of the visco-elastic foundation (iron seam) in GPa; gp is the viscosity coefficient of visco-elastic foundation in Pas; m is mining thickness in m. According to the research results of the reference [24], u(x, t) is compressive stress and can also be a tensile stress. Substituting (11) into (10), can obtain the load of beam is

 qðx; tÞ ¼ Pz 

 gp @ Ep Wðx; tÞ : Wðx; tÞ þ m m @t

ð12Þ

Substituting (12) into (9), we can get the partial differential equation of visco-elastic beam deflection curve:

Igb 

@ 5 Wðx; tÞ @ 4 Wðx; tÞ gp @Wðx; tÞ Ep þ IE  þ þ Wðx; tÞ ¼ Pz : b @x4 @t @x4 @t m m

ð13Þ

For the cantilever part of the visco-elastic foundation beam, the foundation reaction force is zero. Therefore, the partial differential equation of visco-elastic beam deflection curve can be written as the following form.

Igb 

@ 5 Wðx; tÞ @ 4 Wðx; tÞ þ IEb  ¼ Pz : 4 @x @t @x4

ð14Þ

According to the initial and boundary conditions, we can solve equations (13) and (14). Based on the research results of Liu [24] and Li et al. [6], let

P z

X O l

φ(x,t) Z Fig. 4. The foundation reaction force u(x,t).

Please cite this article in press as: W.-X. Li et al., A visco-elastic theoretical model for analysis of dynamic ground subsidence due to deep underground mining, Appl. Math. Modell. (2015), http://dx.doi.org/10.1016/j.apm.2015.01.003

6

W.-X. Li et al. / Applied Mathematical Modelling xxx (2015) xxx–xxx

 Eb mPz Wðx; tÞ ¼ Uðx; tÞ exp  t þ ; gb Ep

ð15Þ

where, U(x, t) is rockmass horizontal displacement due to mining [24]. Substituting (15) into (13), we have

Igb 

 Eb gp @ 5 Wðx; tÞ gp @Wðx; tÞ 1 Uðx; tÞ ¼ 0: Ep  þ þ 4 @x @t @t m m gb

ð16Þ

We assume that: Eb/gb = Ep/gp = k, then the Eq. (16) can be simplified into the form below.

" # @ @ 4 Uðx; tÞ Ep þ Uðx; tÞ ¼ 0: @t @x4 mIEb

ð17Þ

From Eq. (17), we can obtain

@ 4 Uðx; tÞ Ep þ Uðx; tÞ ¼ f ðxÞ; @x4 mIEb

ð18Þ

where, f(x) is an arbitrary function depends on x. The general solution of Eq. (18) is

Uðx; tÞ ¼ expðaxÞ½a1 ðtÞ sinðaxÞ þ a2 ðtÞ cosðaxÞ þ expðaxÞ½a3 ðtÞ sinðaxÞ þ a4 ðtÞ cosðaxÞ þ FðxÞ;

ð19Þ

where, a = [Ep/(4mIEb)]1/4; Eb = E0/(1l2); E0 is the elastic modulus of the material; l is the Poisson’s ratio of the material; F(x) is an arbitrary function depends on f(x). Substituting Eq. (19) into (15), we can get the deflection curve equation of visco-elastic foundation beam.

Wðx; tÞ ¼ expðktÞfexpðaxÞ½a1 ðtÞsinðaxÞ þ a2 ðtÞcosðaxÞ þ expðaxÞ½a3 ðtÞsinðaxÞ þ a4 ðtÞ cosðaxÞ þ FðxÞg þ mPz =Ep : ð20Þ According to the boundary conditions, we can get the following results. (1) When x ? 1, the deflection equation of the beam (roof) is

lim Wðx; tÞ ¼ mPz =Ep :

ð21Þ

x!1

From Eqs. (20) and (21), we have

lim FðxÞ ¼ 0;

a1 ðtÞ ¼ a2 ðtÞ ¼ 0:

x!1

ð22Þ

(2) When t ? 0, the deflection equation of the beam is

limWðx; tÞ ¼ mPz =Ep :

ð23Þ

t!0

From Eqs. (20) and (23), we have

FðxÞ ¼  expðaxÞ½a3 ð0Þ sinðaxÞ þ a4 ð0Þ cosðaxÞ:

ð24Þ

When x ? 1, F(x) = 0. Substituting Eqs. (22 and 24) into (20), we can get

Wðx; tÞ ¼ expðktÞ expðaxÞ½A3 ðtÞ sinðaxÞ þ A4 ðtÞ cosðaxÞ þ mPz =Ep ;

ð25Þ

where,

A3 ðtÞ ¼ a3 ðtÞ  a3 ð0Þ; A4 ðtÞ ¼ a4 ðtÞ  a4 ð0Þ:

ð25aÞ

From Eq. (25a) and the reference [24], we can obtain 2

2

A3 ðtÞ ¼ 

Pz l Pz l expðktÞ þ ; 4IEb a2 4IEb a2

A4 ðtÞ ¼ 

Pz l ð2 þ alÞ Pz l ð2 þ alÞ expðktÞ þ ; 4IEb a2 4IEb a3

2

ð26Þ 2

ð27Þ

where, l is the width of gob (m); a = [Ep/(4mIEb)]1/4. Substituting Eqs. (26 and 27) into (25), we can get the deflection equation of the beam.

Please cite this article in press as: W.-X. Li et al., A visco-elastic theoretical model for analysis of dynamic ground subsidence due to deep underground mining, Appl. Math. Modell. (2015), http://dx.doi.org/10.1016/j.apm.2015.01.003

7

W.-X. Li et al. / Applied Mathematical Modelling xxx (2015) xxx–xxx

" Wðx; tÞ ¼ ½1  expðktÞ expðaxÞ

# 2 Pz l Pz lð2 þ alÞ MPz þ sinð a xÞ cosðaxÞ þ : 4IEb a3 4IEb a2 Ep

ð28Þ

Eq. (28) is the subsidence of roof (rock beam). The last item (MPz/Ep) in Eqs. (5)–(25), (25a), (26)–(34)has no practical engineering significance, because it is happened in prior to mining, and should therefore be deleted it. So we get the following equation.

 p  P L2 l2  p  2L p   z cos Wðx; tÞ ¼ ½1  expðktÞ exp  x l sin þ l x x ; þ L 4IEb p2 L p L

ð29Þ

where, L is the ‘‘half wavelength’’ of stress wave and is an undetermined parameter. Eq. (29) is the subsidence of roof and the final subsidence equation can be written as follows:

 p  P lL2  p  2L p  z cos Wðx; tÞ ¼ ½1  expðktÞ exp  x l sin þ l þ x x : L 4IEb p2 L p L

ð30Þ

According to the assumed conditions, we can choose the coordinate system shown in Fig. 5. Fig. 5 shows the following five points: (1) When x > 0, the foundation of visco-elastic beam is unmined seam; (2) When x < 0, after mining, the void space (mined-out area) generally is filled with waste stone (beam foundation is the filling body). (3) There is a uniform load Pz in the upper surface of the beam. (4) When x > 0, the foundation reaction force is u(x,t). (5) When x < 0, the foundation reaction force is w(x,t). When x > 0, from Eq. (25), we can write the curve equation of viscoelastic beam deflection.

Wðx; tÞ ¼ expðktÞ expðpx=LÞ½A3 ðtÞ sinðpx=LÞ þ A4 ðtÞ cosðpx=LÞ;

ð31Þ

When x < 0, from Fig. 5 and Eq. (1), the reaction force of foundation can be described by the following equation [24].

wðx; tÞ ¼ Ek ek þ gk e_ k ¼

Ek g @ ½Wðx; tÞ  Wð0; tÞ þ k ½Wðx; tÞ  Wð0; tÞ; m m @t

ð32Þ

where, Ek is the elastic modulus of filling body; gk is the viscosity coefficient of filling body; m is the mining thickness (Here, the mining thickness is approximately equal to the thickness of the filling body); W (0, t) is the roof deflection (i.e. subsidence) before filling in mining boundary. Substituting (32) into (10), we can get the load q(x,t) acting on the rock beam.

qðx; tÞ ¼ Pz  wðx; tÞ ¼ Pz 

  Ek g @ ½Wðx; tÞ  Wð0; tÞ þ k ½Wðx; tÞ  Wð0; tÞ : m m @t

ð33Þ

Substituting (33) into (9), we can obtain the partial differential equation of beam deflection curve.

Igb

@ 5 Wðx; tÞ @ 4 Wðx; tÞ gk @ Ek þ IEb þ ½Wðx; tÞ  Wð0; tÞ þ ½Wðx; tÞ  Wð0; tÞ ¼ Pz : 4 @x @t @x4 m @t m

ð34Þ

Considering the boundary conditions and let Eb/gb = Ep/gp = k, from Eq. (34), the roof deflection (subsidence) equation is obtained:

Wðx; tÞ ¼ expðktÞfexpðbxÞ½B1 ðtÞ sinðbxÞ þ B2 ðtÞ cosðbxÞ  mPz =Ek g þ Wð0; tÞ þ mP z =Ek ;

ð35Þ

1/4

where, b = [Ek/(4IEbm)] , Ek is the elastic modulus of filling body (GPa), Eb is the elasticity modulus of the bending beam (GPa), I is the moment of inertia of cross-section (m4), m is the mining thickness (m).

P

Z

O

X

φ(x,t)

Ψ(x,t) Z

Fig. 5. Calculation model for analysis of the beam bending.

Please cite this article in press as: W.-X. Li et al., A visco-elastic theoretical model for analysis of dynamic ground subsidence due to deep underground mining, Appl. Math. Modell. (2015), http://dx.doi.org/10.1016/j.apm.2015.01.003

8

W.-X. Li et al. / Applied Mathematical Modelling xxx (2015) xxx–xxx

2.3. Theoretical models of ‘‘sub-critical mining’’ In general, when the mining depth is much larger than the mining width, the roof subsidence cannot be fully propagate to the ground surface, this mining condition is called the ‘‘sub-critical mining’’ (Fig. 6). Then, the ground maximum subsidence will continue increase with the increment of gob size of underground. The surface subsidence basin is called ‘‘incomplete basin’’. This paper mainly studies the dynamic ground subsidence under the condition of finite mining (sub-critical mining). In fact, one finite mining subsidence is the superposition of two semi-infinite mining subsidences. i.e.

W f ¼ W s1 þ W s2 ; where, Wf is the ground subsidence due to finite mining in the integral interval [S 0]; Ws1 and Ws2 are the subsidences due to semi-infinite mining in the integral intervals (1, 0] and [0, 1), respectively. In interval I (Fig. 6): 1 6 x 6 S, the roof deflection (subsidence) curve can be described by Eq. (13). Let Eb/gb = Ep/gp = k, then we can obtain the solution of differential equation (13) [24]. Considering the boundary conditions, then can obtain the expression of ground subsidence in this interval. When x ? 1, the roof deflection equation is: mPz/Ep, namely

lim Wðx; tÞ ¼ mP z =Ep :

ð36Þ

x!1

When t ? 0, the roof deflection equation is:

limWðx; tÞ ¼ mPz =Ep :

ð37Þ

t!0

Then, we can get the following dynamic subsidence equation in interval 1 6 x 6 S.

Wðx; tÞ ¼ expðLx  ktÞ½A1 ðtÞ  sinðLxÞþA2 ðtÞ  cosðLxÞ þ mPz =Ep ;

ð38Þ

where, A1, A2 are function parameters to be determined. In interval II: S 6 x 6 0, the roof deflection (subsidence) can be described by Eq. (34). When t ? 0, the roof deflection is W (0, 0). Then, we can obtain the expression of ground subsidence in this interval:

Wðx; tÞ ¼ expðktÞfexpðLxÞ½B1 ðtÞ sinðLxÞ þ B2 ðtÞ cosðLxÞ þ expðLxÞ½B3 ðtÞ sinðLxÞ þ ½B4 ðtÞ cosðLxÞg þ Wð0; tÞ þ

mPz ½1  expðktÞ; Ek

ð39Þ

where, B1, B2, B3, B4 are function parameters to be determined. In interval III: 0 6 x 6 1, the roof deflection is Eq. (11). Setting Eb/gb = Ep/gp = k, and combined with the boundary conditions: (1) when x ! 1, the roof deflection is: mPz/Ep, namely, satisfy the following equation:

lim Wðx; tÞ ¼ mP z =Ep :

ð40Þ

x!1

(2) when t ! 0, we have

limWðx; tÞ ¼ mPz =Ep :

ð41Þ

t!0

Thus, we can obtain the expression of ground subsidence in this interval:

Wðx; tÞ ¼ expðLx  ktÞ½C 1 ðtÞ  sinðLxÞþC 2 ðtÞ  cosðLxÞ þ mP z =Ep ;

ð42Þ

where, C1, C2 are function parameters to be determined.

O

X

W Pz

S Z Fig. 6. The ground subsidence for the sub-critical mining.

Please cite this article in press as: W.-X. Li et al., A visco-elastic theoretical model for analysis of dynamic ground subsidence due to deep underground mining, Appl. Math. Modell. (2015), http://dx.doi.org/10.1016/j.apm.2015.01.003

9

W.-X. Li et al. / Applied Mathematical Modelling xxx (2015) xxx–xxx

According to the boundary conditions, we have [24,37]:

Wð0; tÞ ¼ WðS; tÞ;

@Wð0; tÞ @WðS; tÞ @ 2 Wð0; tÞ @W 2 ðS; tÞ ¼ ; ¼ : @x @x @x2 @x2

ð42aÞ

Therefore, we can solve the coefficients (A1, A2, B1, B2, B3, B4, C1, C2) in the equations (38, 39) and (42). Substituting these coefficients into the original equations, and then we can get the final expression of surface subsidence W(x, t).

Wðx; tÞ ¼ W 1 ðx; tÞ þ W 2 ðx; tÞ þ W 3 ðx; tÞ;

ð43Þ

W 1 ðx; tÞ ¼ 0:5W max TðtÞfexp½a1 ðx þ SÞ cos½a1 ðx þ SÞ  expða1 xÞ cosða1 xÞ; 1 < x 6 S;

ð43aÞ

W 2 ðx; tÞ ¼ 0:5W max TðtÞf2  expða1 xÞ cosða1 xÞ  exp½a1 ðx þ SÞ cos½a1 ðx þ SÞg; S < x 6 0;

ð43bÞ

W 3 ðx; tÞ ¼ 0:5W max TðtÞfexpða1 xÞ cosða1 xÞ  exp½a1 ðx þ SÞ cos½a1 ðx þ SÞg; 0 < x < 1;

ð43cÞ

where, TðtÞ ¼ ½1  expðktÞ=½1 þ expðbÞ cos b; a1 = p/L, b = pS/L, L is called ‘‘half-wavelength’’ and is a constant to be determined; Wmax = mn, m is mining thickness, n is subsidence factor. According to the theoretical formulae derived above, we can predict and analyse the dynamic ground subsidence due to underground deep mining. 3. Application of visco-elastic theoretical models to engineering examples In order to demonstrate the application of the visco-elastic theoretical models some examples are given of the practical application of the above theoretical results. The visco-elastic model is verified with the observed dynamic subsidence data from the Chengjiao iron mine, which is located in the eastern part of China. This iron mine is in one of the most important ore regions in China. Mining started in the mine in 1971. The strata are horizontal in the mining area and the lithology of overlying strata is mainly argillaceous siltstone (soft rock strata). The mine consists of two mining areas (West and North), using the pillarless sublevel caving method in deep orebody and the largest mining depth will reach 1100 m. The maximum thickness of overlying rock is 520 m, and the final mining width will exceed 2700 m; the mining width will exceed five times the thickness of the overlying strata. According to the geological conditions in the mine, the roof (overlying strata) is considered as a viscoelastic plate. In the two-dimensional case, the roof is considered as a visco-elastic beam. According to the geological data and test results, the following formula was established [37]:

L ¼ 2R;

ð44Þ

where, L is called half-wavelength in metres; R = H/tand is radius of influence of ground subsidence in metres; H is the overburden thickness (i.e. mining depth) in metres; d is the influence angle of ground subsidence in degrees; tand is a parameter depending on the mining method and the rock properties.According to the analysis results of the geology data in the mining area, the West-1 and North-2 mining areas are belong to the ‘‘sub-critical mining’’ in Chengjiao mine. Mining depth is 512 m and mining thickness is 40 m in West mining area. And mining depth is 520 m and thickness is 50 m in North mining area. Based on the field data in the West-1 and North-2 mining areas, the modelling parameters are given in Table 1. All other engineering parameters are directly taken from geomechanics test results.In general, the finite mining-induced subsidence can be called ‘‘finite subsidence’’ and semi-infinite mining-induced subsidence can be called ‘‘semi-infinite subsidence’’. In the process of calculation, therefore, the finite subsidence is considered as the superposition of two semi-infinite subsidences. Thus, we can draw the following theoretical formulae of the finite mining-induced subsidence (horizontal strata visco-elastic models): 2

Wðx; tÞ ¼ 0:5W max ðk=½ðk þ pm=LÞ þ ðpm=LÞ2 Þ  fexpð2p þ ðpm=LÞtÞ  ½ðk þ pm=LÞ  /ð2p ðpm=LÞtÞ  ðpm=LÞ  sinð2L  v tÞ  ðk þ pm=LÞ  expðkL=m  p  ktÞg; t 2 ½L=m; 2L=m;

ð45Þ

h i 2 Wðx; tÞ ¼ 0:5W max  ðkðk þ gÞ= ðk þ gÞ þ g 2 Þ  ½1  expðk  pÞ  expð2k  ktÞ h i 2  0:5W max  k= ðk  kÞ þ g 2  fexpð2p  gtÞ  ½ðk  gÞ  /ð2p  gtÞ  g  sinð2p  gtÞ  ðk  gÞ  expð2k  ktÞg þ W max  ½1  expð2k  ktÞ; t 2 ½2L=m; 1Þ;

ð46Þ

Table 1 The values of engineering parameters. Miningarea

Mining thickness m (m)

H (m)

L0 (m)

d (°)

Eb (GPa)

gb (GPa s)

n

West-1 North-2

40 50

512 520

560 530

68 66

2.0 2.2

2.60 2.56

0.31 0.28

Please cite this article in press as: W.-X. Li et al., A visco-elastic theoretical model for analysis of dynamic ground subsidence due to deep underground mining, Appl. Math. Modell. (2015), http://dx.doi.org/10.1016/j.apm.2015.01.003

10

W.-X. Li et al. / Applied Mathematical Modelling xxx (2015) xxx–xxx

0

0.2

W/Wmax

0.4

0.6

0.8

1 0

predicted measured 2

4

6

8

10

12

14

time (years) Fig. 7. The data points and the theoretical curve for dynamic subsidence due to deep mining (West 21).

0

0.2

W/Wmax

0.4

0.6

0.8

predicted measured

1 0

2

4

6

8

10

12

14

time (years) Fig. 8. The data points and the theoretical curve for dynamic subsidence due to deep mining (North 56).

where, g = pm/L, k = kL/m; m is the rate of face advance (m/year); k = Eb/gb = Ep/gp; L is the ‘‘half-wavelength’’ of stress wave and is a constant to be determined (m). Based on Eqs. (45) and (46), we can predict the dynamic ground subsidence due to deep mining by pillarless sublevel caving method. The theoretical results are compared with the measured data (observation stations West 21 and North 56; see Figs. 7 and 8). From Figs. 7 and 8, we can draw the following conclusions: (1) From the results measured of the surface observation stations West 21 and North 56, the dynamic ground subsidences have experienced more than 13 years and 14 years, respectively. (2) The engineering actual situation is consistent with the theoretical results. When time tends to infinity, that is considered the final state of surface subsidence. Here, we analyse the ground subsidence of the West mining area in the iron mine. The prediction results of surface subsidence due to underground mining by using the horizontal strata visco-elastic model (VEM) and the stochastic medium model (SMM) are shown in Fig. 9. According to theoretical results and Fig. 9, we can obtain the following conclusions: (1) From Fig. 9 shows, surface uplift phenomenon is appeared after underground mining in Chengjiao ming area. (2) The agreement of the theoretical results with the filed measurements shows that the model is satisfactory. (3) The results of this investigation show that the horizontal strata visco-elastic model (VEM) is more apt to propagate vertical movements due to deep mining by pillarless sublevel caving method, and better able to fit observed data than the stochastic medium model (SMM). Please cite this article in press as: W.-X. Li et al., A visco-elastic theoretical model for analysis of dynamic ground subsidence due to deep underground mining, Appl. Math. Modell. (2015), http://dx.doi.org/10.1016/j.apm.2015.01.003

W.-X. Li et al. / Applied Mathematical Modelling xxx (2015) xxx–xxx

11

0 2

W /m

4 6 8 10

Measured VEM SMM

12 -600 -500 -400 -300 -200 -100

0

100 200 300 400 500

x/m Fig. 9. The data points and the theoretical curves (western mining area).

4. Conclusions According to the results of theoretical analysis, we can get the following conclusions. (1) The predicted results accord with practical engineering situation. The agreement of the theoretical results with the filed measurements shows that the model is satisfactory and the formulae obtained are valid and thus can be effectively used for predicting the ground movement due to deep underground mining by pillarless sublevel caving method. (2) From the comparison between the measured data and the theoretical curves, the proposed visco-elastic model (VEM) to forecast the ground subsidence is in line with the actual situation of mining engineering. (3) The predicted results show that the visco-elastic model (VEM) is better than the classical stochastic medium model (SMM). (4) The parameter values determined (b = 66°–68°, n = 0.28–0.31, gb = 2.56–2.60) in accordance with engineering practice in Chengjiao mine. And the parameter values can provide important scientific basis for later engineering prediction.

Acknowledgements Comments from two anonymous reviewers greatly improved the manuscript. The Project Supported by Natural Science Foundation – Steel and Iron Foundation of Hebei Province is gratefully acknowledged (No. E2011201114). The efforts of all members of the Chengjiao iron mine involved in this project, and their respective support staff in collecting and processing the extensive quantity of data from the monitoring program, is acknowledged. References [1] A. Jarosz, M. Karmis, A. Sroka, Subsidence development with time – experiences from longwall operations in the Appalachian coalfield, Int. J. Min. Geol. Eng. 8 (3) (1990) 261–273. [2] Z. Zeng, X. Kou, Application of viscoelasticity to study the time-dependent surface subsidence caused by underground mining, Eng. Geol. 32 (1992) 279–284. [3] X. Cui, J. Wang, Y. Liu, Prediction of progressive surface subsidence above longwall coal mining using a time function, Int. J. Rock Mech. Min. Sci. 38 (2001) 1057–1063. [4] X.G. Lian, J. Andrew, S.R. Jose, H.Y. Dai, Extending dynamic models of mining subsidence, Trans. Nonferrous Met. Soc. China 21 (2011) 536–542. [5] X. Yu, T. Dang, H. Pan, J. Wu, Theological characteristics of surface dynamic subsidence by mining, J. Xian Univ. Sci. Technol. 23 (2) (2003) 131–134. [6] W.X. Li, Y.G. Guo, X.B. Hou, A visco-elastic model for analysis of ground subsidence due to underground mining by pillarless sublevel caving method, Eng. Mech. 26 (7) (2009) 227–231. [7] T. Wang, X. Yan, X. Yang, H. Yang, Dynamic subsidence prediction of ground surface above salt cavern gas storage considering the creep of rock salt, Sci. China Technol. Sci. 12 (53) (2010) 3197–3202. [8] R.P. Singh, R.N. Yadav, Prediction of subsidence due to coal mining in Raniganj coalfield, West Bengal, India, Eng. Geol. 39 (1995) 103–111. [9] A.J. Astin, A visco-elastic analysis of ground movement due to an advancing coal face, J. Eng. Math. 2 (1968) 9–22. [10] D. Grgic, F. Homand, D. Hoxha, A short- and long-term rheological model to understand the collapses of iron mines in Lorraine, France, Comput. Geotech. 30 (2003) 557–570. [11] H.J. King, J.T. Whetton, Mechanics of mine subsidence, in: Proc. Eur. Congr. Ground Movement, Leeds, 1957, pp. 27–36. [12] H. Hoffman, The effects of direction of working and rate of advance on the scale-deformation of a self-loaded stratified model of a large body of ground, in: Proc. Int. Conf. Strata Control, New York, N.Y., 1964, pp. 397–411. [13] B. Kumar, N.C. Saxena, B. Singh, A new hypothesis for subsidence prediction, J. Mines Met. Fuels 31 (1983) 459–465. [14] W.X. Li, Applications of Fuzzy Mathematics in Mining and Geotechnical Engineering, The Press of Metallurgy Industry, Beijing, 1998. [15] S. Knothe, Observations of surface movements under influence of mining and their theoretical interpretation. in: Proc. Eur. Congr. Ground Movement, Leeds 1957, pp. 210-218.

Please cite this article in press as: W.-X. Li et al., A visco-elastic theoretical model for analysis of dynamic ground subsidence due to deep underground mining, Appl. Math. Modell. (2015), http://dx.doi.org/10.1016/j.apm.2015.01.003

12

W.-X. Li et al. / Applied Mathematical Modelling xxx (2015) xxx–xxx

[16] M. Karmis, A. Jarosz, P. Schilizzi (Monitoring and prediction ground movements above underground mines in the Eastern United States), in: Proceedings of 6th International Conference on Ground Control in Mining, West Virginia University, 1987, pp. 184–194. [17] M. Karmis, A. Jarosz, Z. Agioutantis, Predicting subsidence with a computer, Coal 26 (12) (1989) 54–61. [18] M. Karmis, Z. Agioutantis, A. Jarosz, Subsidence prediction techniques in the United States: a state-of-the-art review, Miner. Res. Eng. 3 (3) (1990) 197– 210. [19] G. Ren, D.J. Reddish, B.N. Whittaker, Mining subsidence and displacement prediction using influence function methods, Min. Sci. Technol. 5 (1987) 89– 104. [20] G. Ren, B.N. Whittaker, D.J. Reddish, Mining subsidence and displacement prediction using influence function methods for steep seams, Min. Sci. Technol. 8 (1989) 235–252. [21] Q.W. Hao, W.M. Ma, Void diffusion models for analysis and prediction of mine subsidence. Rock Mechanics Contributions and Challenges. Proc. U.S. Syrup. Rock Mechanics 31 (1990) 203-210. [22] J. Litwiniszyn, On Certain Linear and Non Linear Strata Theoretical Models. in: Proceedings of the International Conference on Strata Control, New York, N.Y. 1964, pp. 83–394. [23] M.D.G. Salamon (The role of linear models in the estimation of surface ground movements induced by mining tabular deposits), in: J.D. Geddes (Ed.), Proceedings of the Conference on Large Ground Movement and Structures, Pentech Press, Cardiff, 1977, pp. 187–208. [24] B. Liu, Introduction of Mining Rock Mechanics, Hunan Sciences & technology Press, Changsha, 1982. 9. [25] D.S. Berry (Progress in the analysis of ground movement due to mining), in: J.D. Geddes (Ed.), Proceedings of the Conference on Large Ground Movement and Structures, Pentech Press, Cardiff, 1977, pp. 781–810. [26] D.J. Reddish, Study of ground strain in relation to mining subsidence (Ph.D. thesis), University of Nottingham, U.K. 1984. [27] D.Y. Fitzpatrick, Modelling of mining subsidence mechanisms and prediction of ground movements (Ph.D. thesis), University of Nottingham, UK. 1987. [28] T.J. Jones, K.K. Kohli (Subsidence over a room and pillar mine in the Appalachian coal province and the use of subsidence predictive methods-a comparative analysis), in: Proc. U.S. Symp. Rock Mech. 26 (1985) 179–187. [29] W.X. Li, H.N. Li, Fuzzy system models (FSMs) for analysis of rock mass displacement caused by underground mining in soft rock strata, Expert Syst. Appl. 36 (3) (2009) 4637–4645. [30] H.J. Siriwardane, J. Amanat (Numerical modeling of the behavior of overburden rock mass associated with longwall mining), in: Proc. U.S. Symp. Rock Mech. 26 (1983) 171–177. [31] K.B. McNabb, Three dimensional numerical modeling of surface subsidence induced by underground mining. Div. Geomech. CS1RO, Mt. Waverly, Australia, Tech. Rep. 146, 1987, 20 pp. [32] L.R. Alejano et al, FDM predictive methodology for subsidence due to flat and inclined coal seam mining, Int. J. Rock Mech. Min. Sci. 36 (1999) 475–491. [33] A. Tomaz, T. Goran, Prediction of subsidence due to underground mining by artificial neural networks, Comput. Geosci. 29 (5) (2003) 627–637. [34] M. Ma, B. Liu, Z. Fang, The Analyzing of ground surface movement due to mining using 3-D Finite Element method on microcomputers, China Min. Mag. 2 (2) (1993) 70–74. [35] F. Shi, J. He, Numerical simulation of surface movement in Insufficient Mining Region, Jiangxi Nonferrous Met. 24 (4) (2007) 16–17. [36] S.S. Peng, Surface Subsidence Engineering, SME, New York, 1992. [37] Hebei University, Experimental Study on the Mechanical Properties of Soft Rock of Deep Mine (Research Report), 2009, 221 pp. [38] A.D. Kerr, Elastic and viscoelastic foundation models, J. Appl. Mech. 31 (1964) 491–499.

Please cite this article in press as: W.-X. Li et al., A visco-elastic theoretical model for analysis of dynamic ground subsidence due to deep underground mining, Appl. Math. Modell. (2015), http://dx.doi.org/10.1016/j.apm.2015.01.003