# A Fractal and Numerical Simulation Coupled Study of Fracture Network during Coal Mining Excavation

Authors: Yanan Gao; Feng Gao; Man-chu Ronald Yeung
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

This paper features a numerical study that is carried out by using discontinuous deformation method (DDA) and fractal geometry. The configurations of rock strata calculated by DDA were imported into a code that is written by using VC++ called “Fractal” to calculate the fractal dimension of the rock strata. As illustrated, a long wall mining case in China was presented. The relationship of the fractal dimension, excavation length, stress, and movement of strata were discussed. The evolution of fractal dimension can be considered as an index of instability or failure. The method proposed in this paper can be employed to predict the period weighting in long wall mining engineering.


---

## Body

## 1. Introduction

Numerical methods which are rapidly developing and have been widely employed in various engineering fields like geotechnical engineering, mining engineering, civil engineering, and so forth from the past decades [1, 2] can be used to study different geometry and time scales of engineering cases or lab tests and the results obtained by these methods can be repeated with ease. Enormous amount of information extracted from a numerical analysis includes almost all the physical parameters such as displacement, stress state, and energy. The advantages mentioned above rarely exit in theoretical or experimental methods.The important information obtained from a numerical analysis is displacements, stresses, and failure area. Besides, the evolution of fracture network can also be simulated by using discontinuous method, such as UDEC [3] software and DDA (discontinuous deformation method) [4, 5].Study of fracture network is of great significance to underground engineering, such as deep mining, since stability of surrounding rock, the movement of strata, gas, and water flow are all associated with it. Wang et al. [6] set up a physical model to investigate the fracture network evolution of overburden rock mass. They pointed out in their study that there exists a strong relation between the ground pressure and evolution and configuration of the fracture network. Fractal dimension was employed as an index to study the evolution fracture network.The fractal theory was widely used to study the geomaterials, especially rocks, after the pioneer work of Xie and Chen [7]. However, compared to lab tests, which are used to study material properties, the relative larger scale problem such as overburden rock movement is rarely studied except the work of Wang et al. [6]. Inspired by the previous study, we attempt to explore more information based on numerical results and fractal dimension.


## 2. Methodology and Basic Theory

Two tools are used in this study; one is the numerical method—DDA—and the other is the fractal theory and fractal dimension.


### 2.1. Basic Theory of DDA

The DDA is a discrete element method following the definition by Cundall and Hart [3]. The method was developed by Shi [4, 5]; since then it has received considerable attention by the geoengineering community and has been under continuous development and application. The DDA is based on inverse analysis method inspired by an in situ rock block model experiment. It was developed to back-calculate the deformed geometry of the block model from experimentally measured displacements and deformations using a least-square formulation.The displacement function of DDA can be written as(1)
u
=
a
0
+
a
1
(
x
-
x
0
)
+
a
2
(
y
-
y
0
)
,
v
=
b
0
+
b
1
(
x
-
x
0
)
+
b
2
(
y
-
y
0
)
,
where u and v are the horizontal and vertical displacement of the rock block, respectively, and a
i and b
i (i
=
0,1
,
2) are unknowns.In DDA, the large deformation and displacement are calculated by accumulating the small deformation and displacement with a time matching scheme. Based on the definition of Cauchy strain, (1) can be written as
(2)
u
=
u
0
+
ε
x
x
(
x
-
x
0
)
+
(
1
2
γ
x
y
-
r
0
)
(
y
-
y
0
)
,
v
=
v
0
+
(
1
2
γ
x
y
+
r
0
)
(
x
-
x
0
)
+
ε
y
y
(
y
-
y
0
)
,
where, the ε
x
x, ε
y
y, and γ
x
y are horizontal strain, vertical strain, and shear strain, respectively.Here, the deformation of theith block is assumed as [D
i]:
(3)
[
D
i
]
=
[
u
0
,
v
0
,
γ
0
,
ε
x
,
ε
y
,
γ
x
y
]
T
.We thus have the displacement as the following form:(4)
[
u
v
]
=
[
T
i
]
[
D
i
]
=
[
∑
j
=
1
6
t
1
j
d
j
∑
j
=
1
6
t
2
j
d
j
]
,
where
(5)
[
T
i
]
=
[
1
0
-
(
y
-
y
0
)
(
x
-
x
0
)
0
(
y
-
y
0
)
2
0
1
(
x
-
x
0
)
0
(
y
-
y
0
)
(
x
-
x
0
)
2
]
(
u
0
,
v
0
) is the rigid body displacement of centroid (
x
0
,
y
0
); γ
0 is rotation angle.Minimization of the potential energy of the system of blocks, following the FEM convention, results in the following equation:(6)
[
K
11
K
12
K
13
⋯
K
1
n
K
21
K
22
K
23
⋯
K
2
n
K
31
K
32
K
33
⋯
K
3
n
⋮
⋮
⋮
⋱
⋮
K
n
1
K
n
2
K
n
3
⋯
K
n
n
]
[
D
1
D
2
D
3
⋮
D
n
]
=
[
F
1
F
2
F
3
⋮
F
n
]
.
D
j is made of 6 × 1 submatrices that contain the six unknowns of each block; K
i
j is composed of 6 × 6 stiffness submatrices associated with the corresponding degrees of freedom of block j; and F
i is a set of 6 × 1 force submatrices of ith block. K
i
i depends on the material properties of ith block and K
i
j
(
i
≠
j
) on the contacts between blocks.


### 2.2. General Illustration of Fractal Dimension

The fractal is a word from Lingua, which means fracture. Fractal geometry, developed by Mandelbrot et al. [8, 9], allows describing irregular forms which are more complex than Euclidean shapes. Fractal has become a powerful tool for natural structures modeling such as fracture surfaces, cracks, and fault traces [7, 10–14]. Both micro- and macrostructure have fractal property and can be characterized by the fractal dimension. A typical fractal example called Koch curve is shown in Figure 1. The Koch curve has a fractal dimension of 1.262.Figure 1
Koch curve.In this paper we employed the software called “Fractal” to calculate the fractal dimension of fracture network and to study the configurations and evolution of fracture network of overburden rocks during the mining excavation. The flow chart is shown in Figure2.Figure 2
Flow chart of the present study.


## 2.1. Basic Theory of DDA

The DDA is a discrete element method following the definition by Cundall and Hart [3]. The method was developed by Shi [4, 5]; since then it has received considerable attention by the geoengineering community and has been under continuous development and application. The DDA is based on inverse analysis method inspired by an in situ rock block model experiment. It was developed to back-calculate the deformed geometry of the block model from experimentally measured displacements and deformations using a least-square formulation.The displacement function of DDA can be written as(1)
u
=
a
0
+
a
1
(
x
-
x
0
)
+
a
2
(
y
-
y
0
)
,
v
=
b
0
+
b
1
(
x
-
x
0
)
+
b
2
(
y
-
y
0
)
,
where u and v are the horizontal and vertical displacement of the rock block, respectively, and a
i and b
i (i
=
0,1
,
2) are unknowns.In DDA, the large deformation and displacement are calculated by accumulating the small deformation and displacement with a time matching scheme. Based on the definition of Cauchy strain, (1) can be written as
(2)
u
=
u
0
+
ε
x
x
(
x
-
x
0
)
+
(
1
2
γ
x
y
-
r
0
)
(
y
-
y
0
)
,
v
=
v
0
+
(
1
2
γ
x
y
+
r
0
)
(
x
-
x
0
)
+
ε
y
y
(
y
-
y
0
)
,
where, the ε
x
x, ε
y
y, and γ
x
y are horizontal strain, vertical strain, and shear strain, respectively.Here, the deformation of theith block is assumed as [D
i]:
(3)
[
D
i
]
=
[
u
0
,
v
0
,
γ
0
,
ε
x
,
ε
y
,
γ
x
y
]
T
.We thus have the displacement as the following form:(4)
[
u
v
]
=
[
T
i
]
[
D
i
]
=
[
∑
j
=
1
6
t
1
j
d
j
∑
j
=
1
6
t
2
j
d
j
]
,
where
(5)
[
T
i
]
=
[
1
0
-
(
y
-
y
0
)
(
x
-
x
0
)
0
(
y
-
y
0
)
2
0
1
(
x
-
x
0
)
0
(
y
-
y
0
)
(
x
-
x
0
)
2
]
(
u
0
,
v
0
) is the rigid body displacement of centroid (
x
0
,
y
0
); γ
0 is rotation angle.Minimization of the potential energy of the system of blocks, following the FEM convention, results in the following equation:(6)
[
K
11
K
12
K
13
⋯
K
1
n
K
21
K
22
K
23
⋯
K
2
n
K
31
K
32
K
33
⋯
K
3
n
⋮
⋮
⋮
⋱
⋮
K
n
1
K
n
2
K
n
3
⋯
K
n
n
]
[
D
1
D
2
D
3
⋮
D
n
]
=
[
F
1
F
2
F
3
⋮
F
n
]
.
D
j is made of 6 × 1 submatrices that contain the six unknowns of each block; K
i
j is composed of 6 × 6 stiffness submatrices associated with the corresponding degrees of freedom of block j; and F
i is a set of 6 × 1 force submatrices of ith block. K
i
i depends on the material properties of ith block and K
i
j
(
i
≠
j
) on the contacts between blocks.


## 2.2. General Illustration of Fractal Dimension

The fractal is a word from Lingua, which means fracture. Fractal geometry, developed by Mandelbrot et al. [8, 9], allows describing irregular forms which are more complex than Euclidean shapes. Fractal has become a powerful tool for natural structures modeling such as fracture surfaces, cracks, and fault traces [7, 10–14]. Both micro- and macrostructure have fractal property and can be characterized by the fractal dimension. A typical fractal example called Koch curve is shown in Figure 1. The Koch curve has a fractal dimension of 1.262.Figure 1
Koch curve.In this paper we employed the software called “Fractal” to calculate the fractal dimension of fracture network and to study the configurations and evolution of fracture network of overburden rocks during the mining excavation. The flow chart is shown in Figure2.Figure 2
Flow chart of the present study.


## 3. Case Study

A long wall mining case of Jining number 2 mine (China) is selected for study. The length and height of the model as shown in Figure3 are 450 m and 110 m, respectively. The key stratum and the inferior key stratum are marked in Figure 3. The load of overburden rock is assumed as 12 MPa to model the in situ stress of 600 m underground. The excavation process is divided into 6 steps. The length of every step is 30 m. The first step begins at the place that is 150 m away from the left boundary. The key stratum and inferior key stratum, which are often thick and hard strata and play major roles in ground control, are marked in Figure 3 [15–17]. The mechanical properties of the model are listed in Table 1.Table 1
Mechanical properties of rocks.
Rock
Property
E/GPa
μ
Thickness/m
Note
Sand and soil
15
0.29
50
5 m of each layer
Siltstone
36
0.2
10
key stratum
Mudstone
16
0.22
12
2 m of each layer
Fine sandstone
19
0.17
5/7
Inferior key stratum and floor
Coal
8.8
0.23
3
Excavation objects
Coarse sandstone
16
0.20
13
Roof, 3.25 m of each layerFigure 3
Numerical model.


### 3.1. Movement and Stress Distribution

As the excavation may result in the change of fracture network, the configurations of rock blocks after every excavation step were recorded as shown in Figure4.Configurations after excavation.
(a)
(b)
(c)
(d)
(e)
(f)After 30 m excavation, the direct roof collapses (Figure4(a)); after 60 m excavation, the inferior key stratum is fractured due to the excavation disturbance (Figure 4(b)); after 90 m excavation, failure of the key stratum happens for the first time (Figure 4(c)); after 120 m excavation (Figure 4(d)), the fracture of inferior key stratum and the collapse of direct roof occurred consistently, while the fracture of key stratum does not extend; after 150 m excavation (Figure 4(e)), the key stratum was fractured again; during that time and afterwards, the fracture of inferior key stratum and the collapse of direct roof also happened consistently.Figure5 plots the evolution of maximum vertical displacement and the first principal stress of key stratum with respect to excavation step. The maximum vertical displacement increased during the excavation. It increases sharply after the excavation of 60 m. This sharp increase may be due to the fracture of the key stratum (90 m excavation). The increase became moderate after 120 m excavation and almost reached an equilibrium value after 150 m excavation. While the first principal stress had a jump after 90 m excavation, it dropped after 150 m excavation.Figure 5
Plot of first principal stress and displacement of key stratum with excavation length.As mentioned above and illustrated in Figure4, the key stratum is fractured after 90 m and 150 m of excavation. It can be concluded that the sharp increases of both displacement and stress are due to first time fracture of the key stratum. Meanwhile, the fractured blocks had relative movement to the adjacent ones (Figure 3(c)), which can also cause an increase of stress. Afterwards, the subsidence of the key stratum increased sharply and reached the maximum value when excavation was about 150 m. This means that the compaction degree of the rock mass under the key stratum increases before 150 m excavation. Afterwards, the stresses in key stratum released.In general, the first time fracture of the key stratum increased the displacement and the stress; the second time fracture of the key stratum caused stress release. To understand the reason of this phenomenon, we carried out the study on the stress around the work face, which is the key factor during the excavation.Based on Figure6, we can see that the maximum stress at the work face increased consistently before 120 m excavation and had dropped after that. The stress accumulation around the work face was due to the fact that the compaction degree increased before 120 m excavation. The stress accumulation got released after 150 m excavation. This release of stress around the work face caused the stress redistribution of the whole stope. This can be the reason for key stress releases observed in the key stratum in the later stage (after 150 m excavation).Figure 6
Plot of maximum stress around work face and fractal dimension with excavation length.


### 3.2. Fractal Study

In the previous section, the stress and displacement were analyzed. In this section, the fractal dimension of the fracture network was calculated and the relationship between the stress and fractal dimension was studied. To calculate the fracture dimension, the configurations of the fracture network (Figure3) were imported into “Fractal.” The fractal dimension of each image was calculated by using box-counting method [7].


#### 3.2.1. Fractal Dimension and Movement and Stress of Key Stratum

Obviously, the fracture network of the overburden rock changes with the configuration of the rock stratum. The displacement and fractal dimension were thus firstly studied. Figure6 shows that fractal dimension increases before 150 m excavation; afterwards it decreased. According to Figure 5, the displacement of key stratum increased with excavation length and almost reached the maximum value after 150 m of excavation. This implies that the fractal dimension may increase with the subsidence of the overburden rock. Meanwhile, after 150 m of excavation, the compaction degree of the rock mass increased, while the fractal dimension suffered a slight decrease. However, after 180 m of excavation, the displacement of the key stratum does not increase, while the fractal dimension increased again. The increase in fractal dimension may be due to the reconfiguration of the fracture network caused by the movement of key stratum since the key stratum was fractured again after 150 m of excavation. In general, the fractal dimension changes with the subsidence of rock strata and decreases with the compaction degree. The fracture of key stratum can also cause an increase in fractal dimension.For the stress in key stratum, before 150 m of excavation, the stress of key stratum increased with the excavation length. When the fractal dimension decreased, the stress reached its maximum value.


#### 3.2.2. Fractal Dimension and Stress around Work Face

From Figure6, we can see that the stress around the work face and fractal dimension increased before 120 m of excavation, while fractal dimension slightly decreased after 150 m of excavation. The maximum stress around the work face also decreased after 150 m of excavation. This implies that the decrease in fractal dimension may be a forewarning of the release in stress at the work face (periodic weighting). It can be obtained that the maximum principal stress as well as fractal increased before the periodic weighting; the fractal dimension increases when the stress increases again, that is, the beginning of the next weighting period.


## 3.1. Movement and Stress Distribution

As the excavation may result in the change of fracture network, the configurations of rock blocks after every excavation step were recorded as shown in Figure4.Configurations after excavation.
(a)
(b)
(c)
(d)
(e)
(f)After 30 m excavation, the direct roof collapses (Figure4(a)); after 60 m excavation, the inferior key stratum is fractured due to the excavation disturbance (Figure 4(b)); after 90 m excavation, failure of the key stratum happens for the first time (Figure 4(c)); after 120 m excavation (Figure 4(d)), the fracture of inferior key stratum and the collapse of direct roof occurred consistently, while the fracture of key stratum does not extend; after 150 m excavation (Figure 4(e)), the key stratum was fractured again; during that time and afterwards, the fracture of inferior key stratum and the collapse of direct roof also happened consistently.Figure5 plots the evolution of maximum vertical displacement and the first principal stress of key stratum with respect to excavation step. The maximum vertical displacement increased during the excavation. It increases sharply after the excavation of 60 m. This sharp increase may be due to the fracture of the key stratum (90 m excavation). The increase became moderate after 120 m excavation and almost reached an equilibrium value after 150 m excavation. While the first principal stress had a jump after 90 m excavation, it dropped after 150 m excavation.Figure 5
Plot of first principal stress and displacement of key stratum with excavation length.As mentioned above and illustrated in Figure4, the key stratum is fractured after 90 m and 150 m of excavation. It can be concluded that the sharp increases of both displacement and stress are due to first time fracture of the key stratum. Meanwhile, the fractured blocks had relative movement to the adjacent ones (Figure 3(c)), which can also cause an increase of stress. Afterwards, the subsidence of the key stratum increased sharply and reached the maximum value when excavation was about 150 m. This means that the compaction degree of the rock mass under the key stratum increases before 150 m excavation. Afterwards, the stresses in key stratum released.In general, the first time fracture of the key stratum increased the displacement and the stress; the second time fracture of the key stratum caused stress release. To understand the reason of this phenomenon, we carried out the study on the stress around the work face, which is the key factor during the excavation.Based on Figure6, we can see that the maximum stress at the work face increased consistently before 120 m excavation and had dropped after that. The stress accumulation around the work face was due to the fact that the compaction degree increased before 120 m excavation. The stress accumulation got released after 150 m excavation. This release of stress around the work face caused the stress redistribution of the whole stope. This can be the reason for key stress releases observed in the key stratum in the later stage (after 150 m excavation).Figure 6
Plot of maximum stress around work face and fractal dimension with excavation length.


## 3.2. Fractal Study

In the previous section, the stress and displacement were analyzed. In this section, the fractal dimension of the fracture network was calculated and the relationship between the stress and fractal dimension was studied. To calculate the fracture dimension, the configurations of the fracture network (Figure3) were imported into “Fractal.” The fractal dimension of each image was calculated by using box-counting method [7].


### 3.2.1. Fractal Dimension and Movement and Stress of Key Stratum

Obviously, the fracture network of the overburden rock changes with the configuration of the rock stratum. The displacement and fractal dimension were thus firstly studied. Figure6 shows that fractal dimension increases before 150 m excavation; afterwards it decreased. According to Figure 5, the displacement of key stratum increased with excavation length and almost reached the maximum value after 150 m of excavation. This implies that the fractal dimension may increase with the subsidence of the overburden rock. Meanwhile, after 150 m of excavation, the compaction degree of the rock mass increased, while the fractal dimension suffered a slight decrease. However, after 180 m of excavation, the displacement of the key stratum does not increase, while the fractal dimension increased again. The increase in fractal dimension may be due to the reconfiguration of the fracture network caused by the movement of key stratum since the key stratum was fractured again after 150 m of excavation. In general, the fractal dimension changes with the subsidence of rock strata and decreases with the compaction degree. The fracture of key stratum can also cause an increase in fractal dimension.For the stress in key stratum, before 150 m of excavation, the stress of key stratum increased with the excavation length. When the fractal dimension decreased, the stress reached its maximum value.


### 3.2.2. Fractal Dimension and Stress around Work Face

From Figure6, we can see that the stress around the work face and fractal dimension increased before 120 m of excavation, while fractal dimension slightly decreased after 150 m of excavation. The maximum stress around the work face also decreased after 150 m of excavation. This implies that the decrease in fractal dimension may be a forewarning of the release in stress at the work face (periodic weighting). It can be obtained that the maximum principal stress as well as fractal increased before the periodic weighting; the fractal dimension increases when the stress increases again, that is, the beginning of the next weighting period.


## 3.2.1. Fractal Dimension and Movement and Stress of Key Stratum

Obviously, the fracture network of the overburden rock changes with the configuration of the rock stratum. The displacement and fractal dimension were thus firstly studied. Figure6 shows that fractal dimension increases before 150 m excavation; afterwards it decreased. According to Figure 5, the displacement of key stratum increased with excavation length and almost reached the maximum value after 150 m of excavation. This implies that the fractal dimension may increase with the subsidence of the overburden rock. Meanwhile, after 150 m of excavation, the compaction degree of the rock mass increased, while the fractal dimension suffered a slight decrease. However, after 180 m of excavation, the displacement of the key stratum does not increase, while the fractal dimension increased again. The increase in fractal dimension may be due to the reconfiguration of the fracture network caused by the movement of key stratum since the key stratum was fractured again after 150 m of excavation. In general, the fractal dimension changes with the subsidence of rock strata and decreases with the compaction degree. The fracture of key stratum can also cause an increase in fractal dimension.For the stress in key stratum, before 150 m of excavation, the stress of key stratum increased with the excavation length. When the fractal dimension decreased, the stress reached its maximum value.


## 3.2.2. Fractal Dimension and Stress around Work Face

From Figure6, we can see that the stress around the work face and fractal dimension increased before 120 m of excavation, while fractal dimension slightly decreased after 150 m of excavation. The maximum stress around the work face also decreased after 150 m of excavation. This implies that the decrease in fractal dimension may be a forewarning of the release in stress at the work face (periodic weighting). It can be obtained that the maximum principal stress as well as fractal increased before the periodic weighting; the fractal dimension increases when the stress increases again, that is, the beginning of the next weighting period.


## 4. Conclusions

In this paper, the numerical method, discontinuous deformation analysis, was employed to calculate a long wall mining case of Jining number 2 coal mine, and for further study, the code “Fractal” was used to calculate the fractal dimension of the fracture network. The evolution of the movement and stress of the overburden rock as well as the fractal dimension were studied.(1)
The release of stress in key stratum can be due to the stress release in the work face.
(2)
The fractal dimension changes with the subsidence of rock strata and decreases with the compaction degree. The fracture of key stratum can also cause an increase of fractal dimension.
(3)
The decrease of fractal dimension indicates that the stress might have reached the maximum value.
(4)
The decrease of fractal dimension may serve as a forewarning of periodic weighting.The method in this paper extends the results of numerical study and provides a new way to predict the strata movement and ground pressure in mining engineering.



---

*Source: 158194-2014-02-23.xml*