You do not have permission to edit this page, for the following reason:

You are not allowed to execute the action you have requested.


You can view and copy the source of this page.

x
 
1
<!-- metadata commented in wiki content
2
3
4
'''A BP neural network-based micro particle parameters calibration and an energy criterion for applying strength reduction method in MatDEM to evaluate 3D slope stability'''
5
6
Wei Jiang*<sup>1, 2</sup>, Yihong Tan<sup>1</sup>, Jinzhou Yan<sup>3</sup>, Ye Ouyang<sup>1</sup>, Zhaoyu Fu<sup>1</sup>, Qiang Feng<sup>1</sup>
7
8
1.Hubei Key Laboratory of Disaster Prevention and Mitigation, China Three Gorges University, Yichang, Hubei, 443002, China
9
10
2. Key Laboratory of Geological Hazards on Three Gorges Reservoir Area, Ministry of Education, Yichang, Hubei, 443002, China
11
12
3. Yichang Construction Investment & Development Co., Ltd., Yichang, Hubei, 443002, China
13
14
Correspondence: Wei Jiang ([mailto:jiangweilion@ctgu.edu.cn jiangweilion@ctgu.edu.cn])
15
-->
16
==Abstract==
17
18
To enhance the applicability of discrete element method in 3D slope stability analysis, a BP neural network-based micro parameter calibration method and an energy criterion are proposed by taking MatDEM as an example. Firstly, the relationship between the micro particle parameters and the shear strengths of particle aggregate are represented by using the BP neural network. And then the micro particle parameters are obtained for the given shear strengths by using a correction calibration. Next, the energy conversions are investigated for the stable and instable slope models in MatDEM. From a view of practical application, the abrupt in variation tendency and magnitude of the kinetic energy is selected for indicating the emergence of the limit equilibrium state of a slope. Finally, the effectiveness of the proposed improvements is testified by taking Baijiabao landslide as an example. Results verify that the calibration method established in this study is applicable to provide the micro particle parameters when the shear strength is constantly reduced, and the factor of safety determined by the kinetic energy criterion reflects the landslide stability at the global level.
19
20
'''Keywords''': 3D slope stability analysis, MatDEM, strength reduction method, micro parameters, failure criterion
21
22
==1. Introduction==
23
24
Slope stability analysis is a classic problem in geotechnical engineering. Three dimensional (3D) analyses has some significant advantages over two dimensional (2D) analyses in the slope stability analysis that the geometry structure, boundary conditions and spatial variations of a slope is fully took into account. The limit equilibrium method (LEM) mainly developed in the last century is still the most popular methodology among engineers. A considerable amount of research has developed 3D LEM from different point of views such as the method of columns, rigorous solutions [1] or utilizing lower bound theorems [2]. The progress in the development of 3D LEM theory till 2021 has been sumarized by Firincioglu [3]. To achieve an accurate result by 3D LEM, slope and slip surface geometry including geological structures must be precisely determined in the model, which is not always possible in natural instabilities with complex geology. 
25
26
With the advances of computers, numerical analysis methods have become affordable to engineers, such as the finite element method (FEM) [4-5], discrete element method (DEM) [6-7] and others [8-15]. These numerical methods fall into two categories, continuum-based and discontinuum-based numerical methods. The latter kind of methods, e.g. DEM, overcome the limitation of the assumption of macro-continuity in the former kind of methods, and are capable of simulating the generation and development process of the slip surface. Due to these advantages, DEM is becoming increasing popular for slope stability analysis in situation where the continuum-based numerical model is inapplicable [16-18].
27
28
When analyzing the slope stability using numerical methods, the strength reduction method (SRM) [19] is a commonly used strategy for determining the factor of safety (FOS) of a slope. The SRM was initially proposed for FEM analysis of slope stability [20]. SRM continuously reduces the shear strength of soil and rock until a slope reaches a limit equilibrium, and then the reduction factor is considered as the FOS of a slope. Theoretically, SRM can be apllied in DEM similarity to the application of SRM in FEM. But, two critical issues are expected to be addressed when analyzing the slope stability by applying SRM in DEM.
29
30
The first critical issue is how to reduce the shear strength of soil and rock in the DEM modeling. DEM simulates the mechanical properties of soil and rock from the micro point of view, and the relationship between the macro parameters of soil and rock and the micro parameters of particles is complicated. Theoretical and numerical approaches were used to reveal this relationship for soil and rock simulated by different models, e.g. the bonded particle model [21-22], the flat jointed model [23], the close-packed lattice model [24], and the coarse grained particle model [25]. But, a rigid theory is still not available to ensure the magnitude of the relationship between macro parameters and micro parameters [16]. The trial-and-error calibration remains the usual approach to obtain the micro parameters of particles corresponding to the macro parameters [26], even if it exhibits high computing costs. In slope stability analysis by using SRM in DEM, the shear strength is constantly reduced, which means that the micro parameters are necessary to be adjusted correspondingly. To avoid the cumbersome calibration on the micro parameters, the shear strength reduction is usually simplified as the reduction of several micro parameters [27-28] based on the regression analysis of the effect of micro parameter on macro property. The simplification may results in the FOS deviating from its original definition, in view of the complicated relationship between the shear strength and the micro parameters. 
31
32
With the rapid development of neural network technology, recently many researchers attempted to establish an intelligent method to represent the relationship between macro and micro parameters for DEM simulation. Albeit some sophisticated algorithms have been employed for accomplishing this object, e.g. the improved simulated annealing algorithm [29], a sequential quasi-Monte Carlo filter [30] and the non-dominated sorting genetic algorithm [31], the backpropagation (BP) neural network is still the most commonly used algorithm in the previous attempts [32-35] due to its simplicity and widely available software. However, BP neural network is possible to obtain a result with a low precision when the output layer has more neurons than the input layer, even if a considerable amount of samples are collected [35]. Therefore, BP neural network has potential to be used in the reduction of the shear strength in the DEM modeling, while a reliable way is to be developed.
33
34
The second critical issue is how to determine the FOS when applying SRM in DEM to analyze the slope stability. The target is to establish an appropriate criterion to indicate that a slope is in the limit equilibrium state. When analyzing the slope stability by using SRM in FEM, whether the slope is in the limit equilibrium state is usually evaluated by using three criteria: (1) penetration of the plastic zone in the slope, (2) convergence of the unbalance force, and (3) an abrupt change in the displacements of selected characteristic points. These criteria cannot be simply duplicated into the slope stability analysis by using SRM in DEM. In many DEM modelings, the particles are treated as rigid bodies, and thus there will be no plastic zone, which induces that it is impossible to employ the first criterion for the DEM analysis of the slope stability. In case some particles fall or a small region of the slope is deformed, the unbalance forces in the DEM modeling will not convergence, while actually the slope is still stable overall. Thus, it is unreasonable to judge whether the slope is unstable based on the second criterion [36]. Similarly, in views of the fact that the displacement characteristic of some selected points cannot reflect the overall behavior of the entire slope, it is improper to determine whether the entire slope is unstable based on the third criterion [37]. To avoid the limitations of the traditional criteria, the abrupt change in the vector graphic of global displacements [38], the variation coefficient of global displacements [36], and the slope failure extent [16] have been used as the failure criterion in the slope stability analysis by using SRM in DEM. These criteria are capable to represent the slope stability at the global level, but they are inapplicable in some cases. For example, the vector graphic of global displacements may be quite disordered for a complex slope, and it is difficult to evaluate significant changes. Thus, a reasonable criterion is expected to determine the FOS when analyzing the slope stability by using SRM in DEM.
35
36
The fast GPU matrix computing discrete element method (MatDEM) [39] is a newly developed code based on the DEM theory and MATLAB, which has been verified as an effective numerical simulations by many examples [40-43]. Since matrix operations and high-performance GPU calculations are used in the code, its computational efficiency is significantly improved. By taking MatDEM as an example, two improvements are proposed in this study to address the aforementioned two issues. Firstly, to precisely represent the reduction of the shear strengths of the geological material in the MatDEM modeling, the relationship between the macro shear strengths and the micro particle parameters is represented by an overdetermined BP neural network, and a correction calibration is developed to obtain the micro particle parameters corresponding to the original and reduced shear strengths. Then, to reasonably determine the FOS of a slope at the global level, the energy conversion and heat generation during the failure of a slope are fully investigated, and thus an energy criterion is developed for the application of SRM in MatDEM to evaluate the slope stability. Finally, the validity of the proposed strategies is testified by taking Baijiabao landslide situated in Zigui County, Hubei Province of China as a case study.
37
38
==2. Fundamental principles of the MatDEM==
39
40
===2.1 Contact model of particles===
41
42
The basic model of the MatDEM code is a series of 3D close-packed elastic particles as shown in [[#img-1|Figure 1]]a, and these particles interact through the spring forces ([[#img-1|Figure 1]]b).
43
44
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;width:auto;" 
45
|-style="background:white;"
46
|
47
{|
48
|+
49
|-
50
|style="text-align: center;padding:10px;"| [[Image:Review_852894090255-image1.png|250x250px]]
51
|-
52
|style="text-align: center;font-size: 75%;padding-bottom:10px;"|(a) 3D close-packed elastic particles
53
|-
54
|style="text-align: center;padding:10px;"| [[Image:Review_852894090255-image2.png|474px]]
55
|-
56
|style="text-align: center;font-size: 75%;padding-bottom:10px;"|(b) Normal and tangential springs between two particles
57
|}
58
|-
59
| style="background:#efefef;text-align:left;padding:10px;font-size: 75%;"| '''Figure 1'''. The contact model of particles in MatDEM
60
|}
61
62
63
The normal force <math>F_n</math> and the normal relative displacement <math>X_n</math> between two adjacent particles are simulated by a normal spring:
64
65
{| class="formulaSCP" style="width: 100%; text-align: center;" 
66
|-
67
| 
68
{| style="text-align: center; margin:auto;" 
69
|-
70
| <math>F_n=\left\{\begin{array}{c}
71
K_nX_n,\\
72
0,\\
73
K_nX_n,
74
\end{array}\mbox{ }\begin{array}{c}
75
X_n\leq X_b,\mbox{ intact }\\
76
X_n\geq 0,\mbox{ broken}\\
77
X_n<0,\mbox{ broken}
78
\end{array}\mbox{ }\begin{array}{c}
79
\qquad \mbox{(a)}\\
80
\qquad\mbox{(b)}\\
81
\qquad\mbox{(c)}
82
\end{array}\right.</math>
83
|}
84
| style="width: 5px;text-align: right;white-space: nowrap;" | (1)
85
|}
86
87
where <math>K_n</math> is the stiffness of the normal spring, and <math>X_b</math> is the breaking displacement. Initially, the particles are interconnected with their adjacent particles and subjected to tensile or compressive spring forces as Eq. (1a). In case of <math>X_n</math> between the particle pair exceeding <math>X_b</math>, the spring breaks and the tensile force ceases to exist between them (Eq. (1b)). However, the compressive force may act between them when they return to a compressive status (Eq. (1c)). It is noticed that the tensile force is positive, and the compressive force is negative in Eq. (1).
88
89
The shear force <math>F_s</math> and the tangential relative displacement <math>X_s</math> between two particles are simulated by a tangential spring:
90
91
{| class="formulaSCP" style="width: 100%; text-align: center;" 
92
|-
93
| 
94
{| style="text-align: center; margin:auto;" 
95
|-
96
| <math>F_\mbox{s}=K_sX_s</math>
97
|}
98
| style="width: 5px;text-align: right;white-space: nowrap;" | (2)
99
|}
100
101
where <math>K_s</math> is the stiffness of the tangential spring. The tangential spring also has a fracture criterion based on the Mohr-Coulomb yield criterion:
102
103
{| class="formulaSCP" style="width: 100%; text-align: center;" 
104
|-
105
| 
106
{| style="text-align: center; margin:auto;" 
107
|-
108
| <math>F_{smax}=\left\{ \begin{array}{c}
109
F_{s0}-\mu_pF_n,\\
110
-\mu_pF_n,
111
\end{array} \mbox{ }\begin{array}{c}
112
\mbox{intact }\\
113
\mbox{ broken}
114
\end{array}\mbox{ }\begin{array}{c}
115
\qquad \mbox{(a)}\\
116
\qquad\mbox{(b)}
117
\end{array}\right.</math>
118
|}
119
| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
120
|}
121
122
where <math>F_{smax}</math> is the inter-particle shear resistance, <math>F_{s0}</math> is the initial inter-particle shear resistance when no normal force exists, and <math>\mu_p</math> is the inter-particle coefficient of friction. When the tangential force <math>F_s</math> exceeds <math>F_{smax}</math> in Eq. (3a), the tangential spring will break, whereupon the shear force <math>F_s</math> is limited to be less than or equal to the inter-particle shear resistance <math>F_{smax}</math> of the broken spring. If the tangential spring is broken and the magnitude of external shear force exceeds the limit <math>F_{smax}</math> in Eq. (3b), two particles begin slipping, and the slipping friction between them is <math>F_{smax}</math> in Eq. (3b).
123
124
===2.2 Energy system of the MatDEM modeling===
125
126
The MatDEM modeling follows the law of energy conversation and Newton’s law of motion. The energies of a model fall into two categories, mechanical energy and heat. A brief introduction of their compositions is given as follows.
127
128
(1) The mechanical energy in the system consists of elastic potential energy, gravitational potential energy and kinetic energy. The elastic potential energy <math>E_e</math> is the sum of the strain energy of normal and tangential springs between particles:
129
130
{| class="formulaSCP" style="width: 100%; text-align: center;" 
131
|-
132
| 
133
{| style="text-align: center; margin:auto;" 
134
|-
135
| <math>E_e=\frac{1}{2}K_nX_n^2+\frac{1}{2}K_sX_s^2</math>
136
|}
137
| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
138
|}
139
140
141
The gravitational potential energy <math>E_g</math> of a particle is:
142
143
{| class="formulaSCP" style="width: 100%; text-align: center;" 
144
|-
145
| 
146
{| style="text-align: center; margin:auto;" 
147
|-
148
| <math>E_g=mgh</math>
149
|}
150
| style="width: 5px;text-align: right;white-space: nowrap;" | (5)
151
|}
152
153
where <math>m</math> is the particle mass, <math>g</math> denotes the gravity acceleration, and <math>h</math> represents the height above the reference level. The kinetic energy <math>E_k</math> of a particle is computed as:
154
155
{| class="formulaSCP" style="width: 100%; text-align: center;" 
156
|-
157
| 
158
{| style="text-align: center; margin:auto;" 
159
|-
160
| <math>E_k=\frac{1}{2}mv^2</math>
161
|}
162
| style="width: 5px;text-align: right;white-space: nowrap;" | (6)
163
|}
164
165
in which <math>v</math> is the scalar of particle velocity.
166
167
(2) The heat in the system comes from the viscous damping, the spring breaking and the friction. Damping is used in the MatDEM to weaken the elastic wave energy of the model and dissipate the kinetic energy in the system. The damping force <math>\boldsymbol{F}_d</math> is given by
168
169
{| class="formulaSCP" style="width: 100%; text-align: center;" 
170
|-
171
| 
172
{| style="text-align: center; margin:auto;" 
173
|-
174
| <math>\boldsymbol{F}_d=-\eta \boldsymbol{v}</math>
175
|}
176
| style="width: 5px;text-align: right;white-space: nowrap;" | (7)
177
|}
178
179
where <math>\eta</math> is the damping coefficient and <math>\boldsymbol{v}</math> is the vector of particle velocity. Because the time step of the simulation is very small, the particle velocity is assumed to be constant in a step. Viscous heat <math>Q_d</math> generated by viscous damping is calculated by
180
181
{| class="formulaSCP" style="width: 100%; text-align: center;" 
182
|-
183
| 
184
{| style="text-align: center; margin:auto;" 
185
|-
186
| <math>Q_d=-\boldsymbol{F}_d\cdot d\boldsymbol{x}</math>
187
|}
188
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
189
|}
190
191
in which <math>d\boldsymbol{x}</math> denotes the particle displacement in the current time step.
192
193
When an intact spring breaks, the spring force reduces, and the elastic potential energy of the inter-particle normal spring and/or tangential spring will dissipate into heat. Thus, breaking heat is equal to the reduction of the elastic potential energy <math>E_e</math>. The calculation of breaking heat is based on the status of the inter-particle normal force. If the inter-particle normal force is tensile, both the normal and tangential spring forces will reduce to zero, and thus the breaking heat <math>Q_b</math> is the sum of elastic potential energy of both normal and tangential springs:
194
195
{| class="formulaSCP" style="width: 100%; text-align: center;" 
196
|-
197
| 
198
{| style="text-align: center; margin:auto;" 
199
|-
200
| <math>Q_b=\frac{1}{2}K_nX_n^2+\frac{1}{2}K_sX_s^2</math>
201
|}
202
| style="width: 5px;text-align: right;white-space: nowrap;" | (9)
203
|}
204
205
206
If the normal force is compressive, the normal spring force will not be changed, and the tangential spring force reduced from <math>F_{smax}</math> in Eq. (3a) to <math>F_{smax}</math> in Eq. (3b). The breaking heat <math>Q_b</math> is the reduction of the elastic potential energy of the tangential spring, which can be calculated by
207
208
{| class="formulaSCP" style="width: 100%; text-align: center;" 
209
|-
210
| 
211
{| style="text-align: center; margin:auto;" 
212
|-
213
| <math>Q_b=\frac{\left({\left(F_{s0}-{\mu }_pF_n\right)}^2-{\left({\mu }_pF_n\right)}^2\right)}{2K_s}=</math><math>\frac{F_{s0}^2-2{\mu }_pF_nF_{s0}}{2K_s}</math>
214
|}
215
| style="width: 5px;text-align: right;white-space: nowrap;" | (10)
216
|}
217
218
219
Friction heat <math>Q_f</math> generated during the sliding process is defined as the product of the average sliding friction and the effective sliding distance:
220
221
{| class="formulaSCP" style="width: 100%; text-align: center;" 
222
|-
223
| 
224
{| style="text-align: center; margin:auto;" 
225
|-
226
| <math>Q_f=\vert 0.5\left(\boldsymbol{F}_{s1}+\boldsymbol{F}_{s2}\right)\cdot d\boldsymbol{S}\vert </math>
227
|}
228
| style="width: 5px;text-align: right;white-space: nowrap;" | (11)
229
|}
230
231
in which <math>\boldsymbol{F}_{s1}</math> and <math>\boldsymbol{F}_{s2}</math> are the sliding friction forces at respectively the beginning and the end of the current time step, and <math>d\boldsymbol{S}</math> denotes the effective sliding distance. More details about Eq. (11) can be found in [39].
232
233
The total energy of a MatDEM model is the sum of all mechanical energy and heat:
234
235
{| class="formulaSCP" style="width: 100%; text-align: center;" 
236
|-
237
| 
238
{| style="text-align: center; margin:auto;" 
239
|-
240
| <math>E_{total}=E_e+E_g+E_k+Q</math>
241
|}
242
| style="width: 5px;text-align: right;white-space: nowrap;" | (12)
243
|}
244
245
where <math>Q</math> represents the sum of viscous heat <math>Q_d</math>, breaking heat <math>Q_b</math> and friction heat <math>Q_f</math>:
246
247
{| class="formulaSCP" style="width: 100%; text-align: center;" 
248
|-
249
| 
250
{| style="text-align: center; margin:auto;" 
251
|-
252
| <math>Q=Q_d+Q_b+Q_f</math>
253
|}
254
| style="width: 5px;text-align: right;white-space: nowrap;" | (13)
255
|}
256
257
258
Based on the law of energy conversation, the total energy of an isolated system is constant. If the system deforms under external force, the increment of the total energy must be equal to the work done by the external force.
259
260
==3. A BP neural network-based micro particle parameters calibration==
261
262
===3.1 The structure of the BP neural network to predict the shear strength===
263
264
In this study, the BP neural network is used to establish the relationship between the shear strength and micro parameters for MatDEM modeling. BP neural network is a multilayer feed-forward neural network which consists of an input layer, implicit layer and output layer as shown in [[#img-2|Figure 2]] .
265
266
<div id='img-2'></div>
267
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;width:auto;" 
268
|-style="background:white;"
269
|
270
{|
271
|+
272
|-
273
|style="text-align: center;padding:10px;"| [[Image:Review_852894090255-image16.png|354px]]
274
|}
275
|-
276
| style="background:#efefef;text-align:left;padding:10px;font-size: 75%;"| '''Figure 2'''. Typical structure of BP neural network
277
|}
278
279
280
When the value of neurons in the input layer is given, neurons in the implicit layer and the output layer are calculated according to the following forward propagation:
281
282
{| class="formulaSCP" style="width: 100%; text-align: center;" 
283
|-
284
| 
285
{| style="text-align: center; margin:auto;" 
286
|-
287
| <math>u_i^l={\sigma }_i^l\left(\sum_{j=1}^{k_{l-1}}w_{ji}^lu_j^{l-1}+ b_i^l\right)</math>
288
|}
289
| style="width: 5px;text-align: right;white-space: nowrap;" | (14)
290
|}
291
292
where <math>l</math> is the number of layers in the implicit layer, <math>u_i^l</math> denotes the <math>i</math>-th neuron in the <math>l</math>-th layer, <math>{\sigma }_i^l</math> is the activation function, <math>k_l</math> is the number of neurons in the <math>l</math>-th layer, <math>w_{ji}^l</math> is the weight and <math>b_i^l</math>  is the deviation. The training of BP neural network indicates constantly adjusting the weights and deviations to minimize the differences between the output value and the true value.
293
294
The structure of the BP neural network prediction model for the shear strength is plotted in [[#img-3|Figure 3]]. The micro parameters are placed in the input layer and the shear strengths are placed in the output layer. It is noticed that only the micro parameters in terms of the particle contact model is involved in the BP neural network. Although some other micro parameters, e.g. the particle radius and the particle density, and the particle packing pattern have been verified to influence the macro mechanical properties of geotechnical materials in the numerical model [31], their effect is ignored in this study. The reason is that the particle radius and the particle density will be constant when the numerical model of the slope has been established, and all particles are prescribed to be close-packed in MatDEM.
295
296
<div id='img-3'></div>
297
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;width:auto;" 
298
|-style="background:white;"
299
|
300
{|
301
|+
302
|-
303
|style="text-align: center;padding:10px;"| [[Image:Review_852894090255-image18.png|396px]]
304
|}
305
|-
306
| style="background:#efefef;text-align:left;padding:10px;font-size: 75%;"| '''Figure 3'''. BP neural network prediction model for the shear strength
307
|}
308
309
310
Considering that an ultimate aim of this study is to acquire the micro parameters corresponding to the original and reduced shear strength, the established BP neural network seems unfavorable. The reasons are explained as follows. If placing the micro parameters in the output layer and the shear strengths in the input layer, the output layer have more neurons than the input layer, which forms an underdetermined mathematical problem. In this situation, the BP neural network is possible to result in a result of low precision. For example, to calibrate micro particle parameters for Particle Flow Code (PFC) modeling, a BP neural network was built by placing three macro geotechnical mechanical parameters in the input layer and four micro parameters in the output layer [35]. Although the BP neural network database consists of four hundred samples, results of some testing samples were still unsatisfactory. Thus, the micro parameters and the shear strengths are placed in the input layer and the output layer respectively, to promise the validity of the BP neural network. Base on the BP neural network predict model of the shear strength, a calibration of the micro particle parameters will be developed similar to the work in [34].
311
312
===3.2 Virtual 3D direct shear tests to measure the shear strength===
313
314
The direct shear test and the triaxial compression test are two commonly used tests to measure the shear strength of geotechnical materials in laboratory, and their virtual counterparts have been established in DEM simulation by researchers [30,44]. Under the same micro parameter conditions, the two virtual tests are possible to result in different values of the shear strength [44]. Considering that the shear strength of geotechnical materials in the latter case study is obtained by using direct shear test, virtual 3D direct shear test is used to measure the shear strength in this study.
315
316
As shows in [[#img-4|Figure 4]], the virtual 3D direct shear test model consists of two parts, the specimen and the shear box. The specimen is composed of 7636 spherical particles, whose radiuses vary from 0.83 mm to 1.20 mm. The radius of the specimen is 30 mm and the height is 20 mm. The shear box consists of 7374 walls. To simulate the shearing process, the imposed boundary conditions are the horizontal velocity of the lower box and the normal pressure on the top wall. By prescribing the horizontal displacement of the lower box in each time step to be 0.05 mm, the shear stress was recorded until the specimen was finally damaged.
317
318
<div id='img-4'></div>
319
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;width:auto;" 
320
|-style="background:white;"
321
|
322
{|
323
|+
324
|-
325
|style="text-align: center;padding:10px;"| [[Image:Review_852894090255-image19.png|360px]]
326
|}
327
|-
328
| style="background:#efefef;text-align:left;padding:10px;font-size: 75%;"| '''Figure 4'''. Virtual 3D direct shear test model for measuring the shear strength of particle aggregate
329
|}
330
331
332
333
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;width:auto;" 
334
|-style="background:white;text-align:center;padding:10px;"
335
|
336
{|
337
|-
338
| [[Image:Review_852894090255-image19.png|360px]]
339
|}
340
|-
341
| style="background:#efefef;text-align:left;padding:10px;font-size: 75%;"| '''Figure 4'''. Virtual 3D direct shear test model for measuring the shear strength of particle aggregate
342
|}
343
344
345
Given the micro particle parameters as <math>K_n =5.0 \times 10^6</math> N/mm, <math>K_s =1.0\times 10^6</math> N/mm, <math>X_b =1.0\times 10^{-5}</math> m, <math>F_{s0}=1.0\times 10^6</math> Pa and <math>\mu_p=0.20</math>, four normal pressures, 100 kPa, 200kPa, 300kPa and 400 kPa, are applied on the top wall to perform the virtual 3D direct shear test. The shear stress−shear displacement curves under the four normal pressures are plotted in [[#img-5|Figure 5]]a. The trend of the four curves is consistent with each other. The increasing of the imposed normal pressures induces a higher peak value of the shear stress, and the shear displacement where the shear stress reaches its peak value increases slightly, which is similar to the result of the previous virtual direct shear test [44]. The Coulomb formula
346
347
{| class="formulaSCP" style="width: 100%; text-align: center;" 
348
|-
349
| 
350
{| style="text-align: center; margin:auto;" 
351
|-
352
| <math>\tau =C + \sigma \tan\varphi </math>
353
|}
354
| style="width: 5px;text-align: right;white-space: nowrap;" | (15)
355
|}
356
357
is employed to determine the shear strength of the particle aggregate as shown in [[#img-5|Figure 5]]b. Finally, the friction angle <math>\varphi</math> is 32.4°, and the cohesion <math>C</math> is 76 kPa.
358
359
<div id='img-5'></div>
360
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;width:auto;" 
361
|-style="background:white;"
362
|
363
{|
364
|-
365
|style="padding:10px;"| [[Image:Review_852894090255-image21.png|270px]]
366
|-
367
|style="text-align:center;font-size: 75%;padding-bottom:10px;"| (a) Shear stress−shear displacement curves
368
|-
369
|  style="text-align:center;padding:10px;"| [[Image:Review_852894090255-image22.png|246px]]
370
|-
371
|  style="text-align:center;font-size: 75%;padding-bottom:10px;"|(b) Determination of the shear strengths
372
|}
373
|-
374
| style="background:#efefef;text-align:left;padding:10px;font-size: 75%;"| '''Figure 5'''. Virtual 3D direct shear test model for the shear strength measurement
375
|}
376
377
===3.3 A correction calibration of micro particle parameters based on the BP neural network===
378
379
Based on the previous studies on the influence of the micro particle parameters in MatDEM on the mechanical properties of the particle aggregate [39], four levels are considered for <math>F_{s0}</math> and <math>\mu_p</math>, and three levels are considered for <math>K_n</math>, <math>K_s</math> and <math>X_b</math> when establishing the BP neural network database. [[#tab-1|Table 1]] lists the level values of micro particle parameters. By using the full orthogonal combination, 432 combinations of micro particle parameters are generated. And then, the virtual 3D direct shear test is executed to obtain their corresponding shear strengths. Finally, the BP neural network database is composed by 432 groups of micro particle parameters and the corresponding shear strengths.
380
381
<div class="center" style="font-size: 75%;">'''Table 1'''. The level values of micro particle parameters to establish the neural network database</div>
382
383
<div id='tab-1'></div>
384
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;font-size:85%;width:auto;" 
385
|-style="text-align:center"
386
! <math>K_s</math> (N/m) !! <math>K_n</math> (N/m) !! <math>F_{s0}</math> (Pa) !! <math>\mu_p</math> !! <math>X_b</math> (m) 
387
|- style="text-align:center"
388
| <math>1.0\times 10^5</math>, <math>1.0\times 10^6</math>, <math>1.0\times 10^7</math>
389
| <math>5.0\times 10^5</math>, <math>5.0\times 10^6</math>, <math>5.0\times 10^7</math>
390
| <math>1.0\times 10^5</math>, <math>5.0\times 10^5</math>, <math>1.0\times 10^6</math>, <math>5.0\times 10^6</math>
391
|<math>0.10, 0.20, 0.30, 0.40</math>
392
|<math>1.0\times 10^{-6}</math>, <math>1.0\times 10^{-5}</math>, <math>1.0\times 10^{-4}</math>
393
|}
394
395
396
The training of the BP neural network actually successfully establishes the mapping relationship between the micro particle parameters in the input layer and the shear strengths in the output layer, which can be denoted as the following formula.
397
398
{| class="formulaSCP" style="width: 100%; text-align: center;" 
399
|-
400
| 
401
{| style="text-align: center; margin:auto;" 
402
|-
403
| <math>\left(C,\tan\varphi \right)=f\left(K_n,K_s,X_b,F_{s0},{\mu }_p\right)</math>
404
|}
405
| style="width: 5px;text-align: right;white-space: nowrap;" | (16)
406
|}
407
408
409
Assuming that the shear strengths (outputs) are continuous in the micro particle parameter space, the micro particle parameters (inputs) corresponding to the target shear strengths could be resolved similarly to the resolution of the nonlinear equations, because the gradient at any given input could be computed in the BP neural network. In the first place, an objective function is defined as
410
411
{| class="formulaSCP" style="width: 100%; text-align: center;" 
412
|-
413
| 
414
{| style="text-align: center; margin:auto;" 
415
|-
416
| <math>Error={\left(\frac{C-C^\ast }{C^\ast }\right)}^2+{\left(\frac{\tan\varphi -\tan{\varphi }^\ast }{\tan{\varphi }^\ast }\right)}^2</math>
417
|}
418
| style="width: 5px;text-align: right;white-space: nowrap;" | (17)
419
|}
420
421
where <math>C^\ast</math> and <math>{\varphi }^\ast</math> are the target cohesion and the target friction angle, respectively. And then, the gradient descent method is used to adjust the inputs until the objective function have a value smaller than the prescribed threshold value. For a simplicity expression, the input (<math>K_n</math>, <math>K_s</math>, <math>X_b</math>, <math>F_{s0}</math>, <math>\mu_p</math>) is denoted by <math>\boldsymbol{X}</math>, the output (<math>C</math>, <math>\tan\varphi</math>) is denoted by <math>\boldsymbol{Y}</math>, and thus Eq. (16) can be written as <math>\boldsymbol{Y}=f(\boldsymbol{X})</math>. The resolution process is explained as follows:
422
423
(1) Determine the initial exploration input <math>\boldsymbol{X}_0</math> by using the K-means clustering analysis [45]. K-means clustering is a vector quantization method originally from signal processing and can partition <math>n</math> observations into <math>k</math> clusters where each observation belongs to the cluster with the nearest mean. The samples in the neural network database are partitioned into <math>k</math> clusters by using K-means clustering analysis. K is takens as 15 in this study. Distances between the cluster centers and the target are computed by Eq. (17), and the nearest cluster center is selected as <math>\boldsymbol{X}_0</math>.
424
425
(2) Calculate the response outputs <math>\boldsymbol{Y}_i</math> for the exploration point of <math>\boldsymbol{X}_i</math> at step <math>i</math>, and then invoke the objective function. When the objective function has a value smaller than the threshold value, the resolution is terminated and <math>\boldsymbol{X}</math> takes the value of <math>\boldsymbol{X}_i</math>. If not, continue to perform step (3).
426
427
(3) Calculate the gradient of the objective function at <math>\boldsymbol{X}_i</math>, and update <math>\boldsymbol{X}</math> along the descent direction of the objective function. The formulation to update <math>\boldsymbol{X}</math> is
428
429
{| class="formulaSCP" style="width: 100%; text-align: center;" 
430
|-
431
| 
432
{| style="text-align: center; margin:auto;" 
433
|-
434
| <math>\Delta \boldsymbol{X}=-\eta \frac{\partial \left(Error\right)}{\partial \boldsymbol{X}}\vert _{\boldsymbol{X}=\boldsymbol{X}_i},\boldsymbol{X}_{i+1}=</math><math>\boldsymbol{X}_i+\Delta \boldsymbol{X}</math>
435
|}
436
| style="width: 5px;text-align: right;white-space: nowrap;" | (18)
437
|}
438
439
in which <math>\Delta \boldsymbol{X}</math> is the increment used to update <math>\boldsymbol{X}_i</math> and <math>\eta</math> is the step length. <math>\Delta \boldsymbol{X}</math> is obtained by multiplying the negative gradient of the objective function by <math>\eta</math>. On account of the nonlinearity of forward propagation rule, <math>\eta</math> is introduced for performing a step search approximation along the error gradient descent direction.
440
441
(4) Let <math>i=i+1</math> and return to step (2).
442
443
Forward propagation of the established BP neural network can be regarded as a substitution of the virtual 3D direct shear test. The samples in the neural network database are always finite, which indicates that the substitution is certain to be approximate. Therefore, the inputs obtained from the mathematical resolution may have a low accuracy. An additional step called as “verifying test” is used to determine whether the inputs are acceptable. The verifying test performs virtual 3D direct shear test by using the inputs obtained from the mathematical resolution, and calculates the error between the resulted shear strength and the target by using Eq. (17). If the objective function has a value smaller than the threshold value, the inputs are acceptable. Otherwise, a correction should be introduced for obtain a more reasonable input.
444
445
The failure of the verifying test indicates that the current BP neural network is still not accurate enough to obtain a fine mathematical resolution. A feasible strategy to improve the precision of BP neural network is adding some new samples to the database. Although the current BP neural network is not accurate enough, it provides a rough search direction for the resolution and the inputs resulted by the last round of the mathematical resolution must be a bit closer to the optimal solution than the initial inputs. Considering that, a new sample consisting of the resulted inputs in the last round and the corresponding shear strength is added to the database, and then the neural network is retrained. After the retraining of the BP neural network, the mapping relationship in Eq. (16) has been replaced by a new relationship. Then, the mathematical resolution is executed again by selecting the resulted inputs in the last round as the initial exploration input, and the verifying test is performed again. The “mathematical resolution−verifying test−neural network correction” flow is repeated till the objective function has a value smaller than the threshold value in verifying test. The entire flow of the correction method to calibrate the micro particle parameters is plotted in [[#img-6|Figure 6]] .
446
447
<div id='img-6'></div>
448
449
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;width:auto;" 
450
|-style="background:white;"
451
|
452
{|
453
|+
454
|-
455
|style="text-align: center;padding-bottom:10px;"| [[Image:Review_852894090255-image26.png|600px]]
456
|}
457
|-
458
| style="background:#efefef;text-align:left;padding:10px;font-size: 75%;"| '''Figure 6'''. The entire flow of the correction method to calibrate the micro particle parameters
459
|}
460
461
===3.5 Validity of the proposed correction calibration===
462
463
To verify the validity of the proposed calibration method, four different target shear strengths are given, and the optimal micro particle parameters are searched by prescribing the threshold value as <math>2.0\times 10^{-4}</math> for Eq. (17). After three or four rounds of mathematical resolution, the resulted micro particle parameters passed the verifying test. [[#tab-2|Table 2]] lists the micro parameters provided by the mathematical resolution and the shear strengths provided by the verifying test in each round.
464
465
<div class="center" style="font-size: 75%;">'''Table 2'''.  Micro particle parameters and shear strengths in each round</div>
466
467
<div id='tab-2'></div>
468
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;font-size:85%;width:auto;" 
469
|-style="text-align:center"
470
! Target shear strength !! Round !! <math>K_n</math> (N/m) !! <math>K_s</math> (N/m) !! <math>X_b</math> (m) !! <math>F_{s0}</math> (Pa) !! <math>\mu_p</math> !! <math>C</math> (kPa) !! <math>\varphi (^o)</math>
471
|-style="text-align:center"
472
|  rowspan='4' | <math>C^\ast=30.0</math> kPa <br> <math>{\varphi }^\ast =19.6^o</math>
473
|  1
474
|  <math>0.27\times 10^7</math>
475
|  <math>5.27\times 10^6</math>
476
|  <math>1.95\times 10^{-5}</math>
477
|  <math>1.25\times 10^{5}</math>
478
|  <math>0.320</math>
479
|  <math>21.7</math>
480
|  <math>30.2</math>
481
|-style="text-align:center"
482
|  <math>2</math>
483
|  <math>0.47\times 10^7</math>
484
|  <math>1.84\times 10^6</math>
485
|  <math>2.72\times 10^{-5}</math>
486
|  <math>3.67\times 10^{5}</math>
487
| <math> 0.220</math>
488
|  <math>28.4</math>
489
|  <math>22.8</math>
490
|-style="text-align:center"
491
|  <math>3</math>
492
|  <math>2.03\times 10^7</math>
493
|  <math>2.32\times 10^6</math>
494
|  <math>1.84\times 10^{-5}</math>
495
|  <math>6.42\times 10^{5}</math>
496
|  <math>0.143</math>
497
|  <math>29.5</math>
498
|  <math>18.9</math>
499
|-style="text-align:center"
500
|  <math>4</math>
501
|  <math>1.27\times 10^7</math>
502
|  <math>2.09\times 10^6</math>
503
|  <math>2.04\times 10^{-5}</math>
504
|  <math>5.95\times 10^{5}</math>
505
|  <math>0.141</math>
506
|  <math>30.2</math>
507
|  <math>19.6</math>
508
|-style="text-align:center"
509
| rowspan='3' |<math>C^\ast=28.5</math> kPa <br> <math>{\varphi }^\ast = 18.6^o</math>
510
|<math> 1 </math>
511
|<math> 2.15\times 10^7 </math>
512
|<math> 4.92\times 10^6 </math>
513
|<math>2.36\times 10^{-5}</math>
514
|<math>3.29\times 10^{5}</math>
515
|<math>0.114</math>
516
|<math>26.3</math>
517
|<math>16.7</math>
518
|-style="text-align:center"
519
|<math>2</math>
520
|<math>1.63\times 10^7 </math>
521
|<math>4.40\times 10^6 </math>
522
|<math>2.02\times 10^{-5}</math>
523
|<math>3.17\times 10^{5}</math>
524
|<math>0.127</math>
525
|<math>29.2</math>
526
|<math>17.9</math>
527
|-style="text-align:center"
528
|<math>3</math>
529
|<math>2.38\times 10^7 </math>
530
|<math>4.75\times 10^6 </math>
531
|<math>1.61\times 10^{-5}</math>
532
|<math>4.52\times 10^{5}</math>
533
|<math>0.129</math>
534
|<math>28.6</math>
535
|<math>18.5</math>
536
|-style="text-align:center"
537
| rowspan='4' |<math>C^\ast=45.0</math> kPa <br> <math>{\varphi }^\ast =28.0^o</math>
538
|<math>1</math>
539
|<math>4.03\times 10^7 </math>
540
|<math>3.82\times 10^6 </math>
541
|<math>5.25\times 10^{-5}</math>
542
|<math>2.07\times 10^{5}</math>
543
|<math>0.324</math>
544
|<math>48.9</math>
545
|<math>31.5</math>
546
|-style="text-align:center"
547
|<math>2</math>
548
|<math>2.94\times 10^7 </math>
549
|<math>4.43\times 10^6 </math>
550
|<math>3.85\times 10^{-5}</math>
551
|<math>2.18\times 10^{5}</math>
552
|<math>0.296</math>
553
|<math>45.7</math>
554
|<math>29.7</math>
555
|-style="text-align:center"
556
|<math>3</math>
557
|<math>3.16\times 10^7 </math>
558
|<math>5.26\times 10^6 </math>
559
|<math>4.19\times 10^{-5}</math>
560
|<math>1.95\times 10^{5}</math>
561
|<math>0.268</math>
562
|<math>45.3</math>
563
|<math>28.9</math>
564
|-style="text-align:center"
565
|<math>4</math>
566
|<math>3.25\times 10^7 </math>
567
|<math>4.47\times 10^6 </math>
568
|<math>4.37\times 10^{-5}</math>
569
|<math>2.23\times 10^{5}</math>
570
|<math>0.272</math>
571
|<math>44.8</math>
572
|<math>27.8</math>
573
|-style="text-align:center"
574
| rowspan='4' |<math>C^\ast =12.0 </math> kPa <br> <math>{\varphi }^\ast =33.0^o</math>
575
|<math>1</math>
576
|<math>3.27\times 10^7 </math>
577
|<math>4.13\times 10^6 </math>
578
|<math>1.74\times 10^{-5}</math>
579
|<math>2.95\times 10^{5}</math>
580
|<math>0.330</math>
581
|<math>18.5</math>
582
|<math>30.5</math>
583
|-style="text-align:center"
584
|<math>2</math>
585
|<math>3.65\times 10^7 </math>
586
|<math>3.29\times 10^6 </math>
587
|<math>1.86\times 10^{-5}</math>
588
|<math>3.20\times 10^{5}</math>
589
|<math>0.306</math>
590
|<math>16.1</math>
591
|<math>28.5</math>
592
|-style="text-align:center"
593
|<math>3</math>
594
|<math>2.45\times 10^7 </math>
595
|<math>2.94\times 10^6 </math>
596
|<math>1.66\times 10^{-5}</math>
597
|<math>2.41\times 10^{5}</math>
598
|<math>0.358</math>
599
|<math>14.3</math>
600
|<math>32.0</math>
601
|-style="text-align:center"
602
|<math>4</math>
603
|<math>2.36\times 10^7 </math>
604
|<math>3.15\times 10^6 </math>
605
|<math>1.57\times 10^{-5}</math>
606
|<math>2.46\times 10^{5}</math>
607
|<math>0.379</math>
608
|<math>11.9</math>
609
|<math>33.2</math>
610
|}
611
612
613
Results in [[#tab-2|Table 2]] shows that <math>Error</math> in Eq. (17) decreases whith the mathematical resolution round increasing, which implies that the target shear strengths are approximated gradually. Taking the case <math>C^\ast</math>=30.0 kPa, <math>{\varphi }^\ast=19.6^o</math> as an example, the relative errors of the resulted shear strengths are respectively 27.7% and 63.5% for the cohesion and the friction coefficient in the first round, by using the target shear strengths as the reference. In the second round, they decrease to 5.3% and 18.1% for the cohesion and the friction coefficient, respectively. Ultimately, they decrease to 0.7% and 0.2% in the fourth round for the cohesion and the friction coefficient, respectively. Thus, the proposed method has the ability to obtain the micro particle parameters that well reflects the shear strengths of particle assembly.
614
615
When performing the slope stability analysis by applying SRM in MatDEM, the shear strengths of the geological material in the MatDEM modeling are constantly adjusted. Once a new value of the shear strengths is given, a new combination of micro particle parameters is to be searched. At the moment, the established neural network database will provide a new initial exploration input <math>\boldsymbol{X}_0</math> by using the K-means clustering analysis, and then the proposed method is performed to obtain the micro particle parameters.
616
617
==4. An energy criterion to evaluate the slope stability==
618
619
As mentioned in Section 2.2, the MatDEM modeling follows the law of energy conversation and Newton’s law of motion. When slope instability happens, the displacements of soil and rock will induce the conversion between various kinds of energies. Consequently, the energy conversion of an instable slope in MatDEM must be different with that of a stable slope. In this section, a simple 2D slope is taken as an example to investigate the energy conversion difference between the stable and instable slope. And then, an energy criterion is established to evaluate the slope stability.
620
621
===4.1 Investigation on the energy conversion difference between the stable and instable slope===
622
623
As  shows in [[#img-7|Figure 7]], the slope is 7 m height and the ratio of slope is 7: 3. And, the slope is assumed to be homogeneous and composed by one kind of geotechnical material for simplicity. The MatDEM model of the slope contains 6305 particles. The radius of particles varies from 0.064 m to 0.096 m. The density is taken as 2.04&#x00d7;10<sup>3</sup> kg/m<sup>3</sup>, and the gravity acceleration <math>g= 9.8</math> m/s<sup>2</sup>. All particles at the left, right and bottom boundaries are fixed as the boundary conditions of the model.
624
625
<div id='img-7'></div>
626
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;width:auto;" 
627
|-style="background:white;"
628
|
629
{|
630
|+
631
|-
632
|style="text-align: center;padding:10px;"| [[Image:Review_852894090255-image27.png|300px]]
633
|}
634
|-
635
| style="background:#efefef;text-align:left;padding:10px;font-size: 75%;"| '''Figure 7'''. A simple homogeneous slope 2D model
636
|}
637
638
639
In the numerical simulation of the slope model, three groups of micro particle parameters listed in [[#tab-3|Table 3]] are selected, which are denoted as Case 1, 2 and 3 respectively. Considering that the chief purpose here is to reveal the energy conversion difference between the stable and instable slope, the shear strengths are not measured for the three groups of micro particle parameters by using the virtual numerical test.
640
641
<div class="center" style="font-size: 75%;">'''Table 3'''. Micro particle parameters to be used in the investigation of energy conversion</div>
642
643
<div id='tab-3'></div>
644
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;font-size:85%;width:auto;" 
645
|-style="text-align:center"
646
! Parameter !! <math>K_n</math> (N/m) !! <math>K_s</math> (N/m) !! <math>X_b</math> (m) !! <math>F_{s0}</math> (Pa) !! <math>\mu_p</math>
647
|-style="text-align:center"
648
|  Case 1
649
|  <math>7.80\times 10^6</math>
650
|  <math>3.10\times 10^6</math>
651
|  <math>1.10\times 10^{-4}</math>
652
|  <math>5.46\times 10^5</math>
653
|  <math>0.020</math>
654
|-style="text-align:center"
655
|  Case 2
656
|  <math>2.80\times 10^8</math>
657
|  <math>1.50\times 10^8</math>
658
|  <math>2.60\times 10^{-5}</math>
659
|  <math>4.90\times 10^5</math>
660
|  <math>0.120</math>
661
|-style="text-align:center"
662
|  Case 3
663
|  <math>9.60\times 10^7</math>
664
|  <math>2.80\times 10^7</math>
665
|  <math>9.30\times 10^{-6}</math>
666
|  <math>3.90\times 10^5</math>
667
|  <math>0.309</math>
668
|}
669
670
671
After 10 time steps of the computation, the slope is verified to be instable in Case 1 but stable in Case 2 and 3. It is noticed that one time step in MatDEM consists of a number of sub-time steps, and the value of a sub-time step is self-determined by the software. In the computation of the slope model, one time step is composed of 63 sub-time steps, and a sub-time step takes 0.0001<math>s</math>, 0.00015<math>s</math> and 0.0002<math>s</math> for Case 1, Case 2 and Case 3, respectively. The displacement cloud charts are plotted in [[#img-8|Figure 8]]a and [[#img-9|Figure 9]]a for Case 1 and Case 2, respectively. The displacement cloud chart of Case 3 is not provided because the displacement difference between Case 2 and Case 3 is not noticeable. [[#img-8|Figure 8]]b and [[#img-9|Figure 9]]b illustrate the variations of the total energy, elastic potential energy, kinetic energy, gravitational potential energy and heat with the time steps for Case 1 and Case 2, respectively. Some interesting phenomena are observed in the energy-time steps curves, and a further investigation is valuable for the establishment of an energy criterion to evaluate the slope stability.
672
673
<div id='img-8'></div>
674
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;width:auto;" 
675
|-style="background:white;"
676
|
677
{|
678
|+
679
|-
680
|style="text-align: center;padding:10px;"| [[Image:Review_852894090255-image28-c.png|324px]]
681
|-
682
|style="text-align: center;font-size: 75%;padding-bottom:10px;"|(a) Displacement cloud chart
683
|-
684
|style="text-align: center;padding:10px;"| [[Image:Review_852894090255-image29.png|312px]]
685
|-
686
|style="text-align: center;font-size: 75%;padding-bottom:10px;"|(b) Energy variations with time step
687
|}
688
|-
689
| style="background:#efefef;text-align:left;padding:10px;font-size: 75%;"| '''Figure 8'''. Displacements and the energy conversion for the slope in Case 1
690
|}
691
692
693
<div id='img-9'></div>
694
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;width:auto;" 
695
|-style="background:white;"
696
|
697
{|
698
|+
699
|-
700
|style="text-align: center;padding:10px;"| [[Image:Review_852894090255-image30-c.png|330px]]
701
|-
702
|style="text-align: center;font-size: 75%;padding-bottom:10px;"|(a) Displacement cloud chart
703
|-
704
|style="text-align: center;padding:10px;"| [[Image:Review_852894090255-image31.png|306px]]
705
|-
706
|style="text-align: center;font-size: 75%;padding-bottom:10px;"|(b) Energy variations with time step
707
|}
708
|-
709
| style="background:#efefef;text-align:left;padding:10px;font-size: 75%;"| '''Figure 9'''. Displacements and the energy conversion for the slope in Case 2
710
|}
711
712
713
Firstly, the conclusion that MatDEM follows the law of energy conversation is confirmed again by the energy-time steps curves. For both two cases, the total energies are almost constant during the computation, no matter whether the slope is stable or not. However, the total energy of Case 2 has a larger value than that of Case 1, which can be explained as follows. In the MatDEM modeling, the gravitational potential energy is evaluated by the relative vertical displacement of particles. It is initialized to be zero for a particle, and becomes negative when the particle falls while positive when the particle rises. The kinetic energy is zero for a particle unless the motion of the particle happens. The heat induced by the viscous damping, the spring breaking and the friction is zero in the initial state. Therefore, the total energy should be equal to the elastic potential energy initially. Since the elastic potential energy is computed as Eq. (4), <math>K_n</math> and <math>K_s</math> dominate the value of the elastic potential energy. Because the values of <math>K_n</math> and <math>K_s</math> in Case 2 is greater than that in Case 1, the total energy of Case 2 has a larger value than that of Case 1.
714
715
Secondly, the variations of the kinetic energy and the gravitational potential energy in [[#img-8|Figure 8]]b and [[#img-9|Figure 9]]b are apparently different. With the time step increasing, the kinetic energy increases but the gravitational potential energy decreases in [[#img-8|Figure 8]]b. But, the gravitational potential energy and the kinetic energy are nearly constant in [[#img-9|Figure 9]]b. [[#img-10|Figure 10]]  plots the variations of the kinetic energy and the gravitational potential energy with time steps for all three cases.
716
717
<div id='img-10'></div>
718
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;width:auto;" 
719
|-style="background:white;"
720
|
721
{|
722
|+
723
|-
724
|style="text-align: center;padding:10px;"| [[Image:Review_852894090255-image32.png|306px]]
725
|-
726
|style="text-align: center;font-size: 75%;padding-bottom:10px;"|(a) Kinetic energy
727
|-
728
|style="text-align: center;padding:10px;"| [[Image:Review_852894090255-image33.png|306px]]
729
|-
730
|style="text-align: center;font-size: 75%;padding-bottom:10px;"|(b) Gravitational potential energy
731
|}
732
|-
733
| style="background:#efefef;text-align:left;padding:10px;font-size: 75%;"| '''Figure 10'''. Variations of the kinetic energy and the gravitational potential energy with time steps for all three cases
734
|}
735
736
737
For Case 3, the kinetic energy is observed, and the gravitational potential energy is positive. That may be caused by the slope deformation. The kinetic energy and the gravitational potential energy caused by the deformation can be easily distinguished from that caused by the failure. Above all, the magnitude of the energy induced by the deformation is tiny compared to that induced by the failure. For Case 3, the kinetic energy and the gravitational potential energy have peak values of 2800 J and 4700 J, respectively. For Case 2, they have peak values of 116 J and 945 J, respectively. Next, when the two energies are caused by deformation, their variations show the convergent trend. But, when the two energies are induced by the failure, their variations show the divergent trend. Since the clear difference in the variation of the kinetic energy and the gravitational potential energy can be observed between a stable slope and an instable slope, they have the potential to be used in the evaluation of the slope stability.
738
739
Finally, the heat is investigated. The variations of heat are illegible in [[#img-8|Figure 8]]b and [[#img-9|Figure 9]]b because the magnitude of heat is too small. [[#img-11|Figure 11]]a plots the variations of the heat with time steps for all three cases. Obviously, the heats in Case 2 and Case 3 are much greater than that in Case 1. But it is incorrect to have a conclusion that the heat in a stable MatDEM slope model is more than that in an instable MatDEM slope model. The heat in the system is composed of the contributions of the viscous damping, the spring breaking and the friction force, in which the latter two sources are calculated based on the micro particle parameters. In [[#img-11|Figure 11]]a, the heat curve for Case 1 is situated in the lowest position, which maybe resulted by the micro particle parameters in [[#tab-3|Table 3]]. Thus, the magnitude of the heat is unreasonable to be chosen as a reference in the slope stability assessment. But, the shapes of the three heat curves are quite different. For a better demonstration, the heat curve is repainted for Case 1 in [[#img-11|Figure 11]]b. The heat curves of Case 2 and Case 3 are roughly convex, but the heat curve of Case 1 is concave. If the slope is stable, the position change of the particles induced by the deformation will stop rapidly, which promises the heat to be stable after several time steps. If the slope is sliding, the positions of the particles are still continuously changed, which results in the on-going increasing of the heat. Therefore, the shape of the heat curve may be effective in the evaluation of the slope stability.
740
741
<div id='img-11'></div>
742
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;width:auto;" 
743
|-style="background:white;"
744
|
745
{|
746
|+
747
|-
748
|style="text-align: center;padding:10px;"| [[Image:Review_852894090255-image34.png|312px]]
749
|-
750
|style="text-align: center;font-size: 75%;padding-bottom:10px;"|(a) Three cases
751
|-
752
|style="text-align: center;padding:10px;"| [[Image:Review_852894090255-image35.png|300px]]
753
|-
754
|style="text-align: center;font-size: 75%;padding-bottom:10px;"|(b) Case 1
755
|}
756
|-
757
| style="background:#efefef;text-align:left;padding:10px;font-size: 75%;"| '''Figure 11'''. Variations of the heat with time steps [68]
758
|}
759
760
===4.2 A criterion based on the kinetic energy of the MatDEM model===
761
762
As mentioned in the Section 4.1, the variation of the kinetic energy and the gravitational potential energy and the shape of the heat curve are verified to be different between a stable slope and an instable slope when performing the simulation by MatDEM. Their validities are compared in this section from a view of practical application.
763
764
When performing the slope stability analysis by using SRM in MatDEM, the micro parameters of the particles will be modified according to the reduction factor. Both the gravitational potential energy and the kinetic energy are independent of the micro particle parameters. However, the heat is calculated based on the micro particle parameters. When the energy variations resulted by using different micro particle parameters are collected together, it is more suitable to compare the results of the gravitational potential energy and the kinetic energy than that of the heat. Additionally, the number of the time steps has a great influence on the reliability of the conclusion if estimating the stable state of a slope by using the shape of the heat curves. As shows in [[#img-11|Figure 11]]a, the overall trend of the heat curve is convergent for Case 2, while a slight turning happens at the sixth time step. A large number of time steps seem to be feasible to overcome the effect of the local abnormality, but it will cause another issue. In fact, the heat curve must be convergent ultimately no matter whether the slope is stable or not. In case of an instable slope, the position of the particles will be relocated finally when the sliding is terminated, and thus the heat will be also stable. Therefore, an incorrect conclusion may be drawn based on the shape of the heat curve when the number of time steps is inappropriate. In short, the gravitational potential energy and the kinetic energy are more suitable to be compared than the heat when estimating the stable state of a slope.
765
766
There is an essential association between the gravitational potential energy and the kinetic energy. The motion of the particles results in the change of their positions, so the kinetic energy always emerges before the gravitational potential energy in MatDEM. This deduction can be verified by the variations of the kinetic energy and the gravitational potential energy for Case 3. As shown in [[#img-10|Figure 10]] , for Case 3 it is the fourth time step that the kinetic energy begins to decrease, while it is the seventh time step that the gravitational potential energy begins to decrease. This phenomenon indicates that the convergence of the kinetic energy should be earlier than that of the gravitational potential energy for a stable slope. Moreover, the variation of the gravitational potential energy is noticed to be small when the time step ranges from four to seven for Case 1. The reason is that the motion of the particles in the model maybe temporarily dominated by the horizontal motion during the failure, but only the vertical motion of the particles induces the variation of the gravitational potential energy. Thus, the gravitational potential energy may have a more complicated curve than the kinetic energy. Based on the above analysis, the kinetic energy is taken as the main reference to estimate the stable state of a slope.
767
768
When performing the slope stability analysis by using SRM in MatDEM, the variations of the kinetic energy are collected under different reduction factors. Then, the variation tendencies and the magnitude of the kinetic energies are compared. The abrupt in the variation tendencies and the magnitude is considered as a sign for the emergence of the critical point, and the safe factor of the slope is determined finally. At this point, a criterion has been proposed to evaluate the slope stability based on the kinetic energy.
769
770
==5 Case study: Baijiabao landslide==
771
772
===5.1 Geological background===
773
774
Baijiabao landslide is situated in Guizhou village of Zigui County, Hubei Province of China on the right bank of the Xiangxi River (30<sup>o</sup>58′59.9′′N, 110<sup>o</sup>45′33.4′′E). A photograph of the landslide that was taken from the opposite bank of the Xiangxi River is shown in [[#img-12|Figure 12]]. The elevation of the landslide ranges from 125 m to 265 m above MSL, and the slope of the landslide ranges from 10<sup>o</sup> to 20<sup>o</sup>. The upper boundary of the landslide is defined by the interface between the bedrock and the soil. The left and right boundaries are defined by two natural gullies that contain many surface cracks. The toe of the landslide varies in elevation between 125 m to 135 m. The primary sliding direction of Baijiabao landslide varies between 75<sup>o</sup> and 85<sup>o</sup> (SW-NE). This landslide is 550 m long and 400 m wide, covers an area of 20 ha, and has an estimated volume of nearly 1 million m<sup>3</sup>. The landslide mass is primarily composed of loose Quaternary deposits. Its sliding zone is mainly silty clay. And the underlying bedrock contains quartz sandstone and argillaceous siltstone of the Early Jurassic Xiangxi Formation, which dips into the hill at an angle of 30<sup>o</sup> to 40<sup>o</sup> ([[#img-13|Figure 13]]).
775
776
<div id='img-12'></div>
777
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
778
|-
779
|style="padding:10px;"| [[Image:Review_852894090255-image36.png|600px]]
780
|- style="text-align: center; font-size: 75%;"
781
| colspan="1" style="padding:10px;"| '''Figure 12'''. A photograph of Baijiabao landslide, taken from the opposite bank of the Xiangxi River
782
|}
783
784
785
<div id='img-13'></div>
786
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
787
|-
788
|  style="text-align: center;padding:10px;"| [[Image:Review_852894090255-image37.png|414px]]
789
|-
790
|  style="text-align: center;font-size: 75%;"|(a) Engineering geological map    
791
|-
792
|  style="text-align: center;padding:10px;"| [[Image:Review_852894090255-image38.png|414px]]
793
|-
794
|  style="text-align: center;font-size: 75%;"|(b) Geological cross-section
795
|- style="text-align: center; font-size: 75%;"
796
| colspan="1" style="padding:10px;"| '''Figure 13'''. Geological background of Baijiabao landslide [46]
797
|}
798
799
800
Baijiabao landslide poses significant threats to public safety. The Zi-Xing road crosses the central part of the landslide mass as [[#img-13|Figure 13]]a shows. Before Three Gorges Reservoir (TGR) impoundment, 165 residents lived in the landslide area, while only 20 residents live there today. The deformation rate of the landslide was observed to increase significantly after TGR impoundment. Monitoring system installed in October 2006 consists of four GPS monitoring sites, which are noted as ZG323, ZG324, ZG325 and ZG326 in [[#img-13|Figure 13]]a. The detailed deformation history of Baijiabao landslide can be referred to the related literatures [46-48].
801
802
===5.2 Landslide model in the MatDEM modeling===
803
804
This subsection describes the establishment of the 3D MatDEM model for the stability analysis of Baijiabao landslide. Firstly, the geomorphy of the landslide is simulated by using ArcGIS software based on the contour line measured in 2006 ([[#img-14|Figure 14]]a). Then, a 3D model of Baijiabao landslide is built in MatDEM. The model is 564 m long, 497 m wide and 186 m height, covering the whole landslide mass ([[#img-14|Figure 14]]b). Next, geotechnical materials are allocated according to the engineering geologic investigation data. The sliding zone in the model ([[#img-14|Figure 14]]c) is positioned by using the interpolation of locations where the sliding zone is encountered in the drilling exploration. The landslide mass has a thickness varying between 20 m to 30 m in the front area, and has a thickness varying between 10 m to 40 m in the trailing area.
805
806
<div id='img-14'></div>
807
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
808
|-
809
|  style="text-align: center;padding:10px;"| [[Image:Review_852894090255-image39.png|198px]]
810
|-
811
|  style="text-align: center;font-size: 75%;"|(a) Geomorphy model
812
|-
813
|  style="text-align: center;padding:10px;"| [[Image:Review_852894090255-image40.png|360px]]
814
|-
815
|  style="text-align: center;font-size: 75%;"|(b) 3D MatDEM model
816
|-
817
|  style="text-align: center;padding:10px;"| [[Image:Review_852894090255-image41.png|318px]]
818
|-
819
|  style="text-align: center;font-size: 75%;"|(c) An enlargement to exhibit the sliding zone
820
|- style="text-align: center; font-size: 75%;"
821
| colspan="1" style="padding:10px;"| '''Figure 14'''. 3D MatDEM model for the stability analysis of Baijiabao landslide
822
|}
823
824
825
Finally, the 3D MatDEM model of Baijiabao landslide is composed by 1,207,985 particles, in which the landslide mass and the sliding zone are represented by 733,104 removable particles, and the bedrock is represented by 474,881 fixed particles. The radius of particles varies between 0.55 m and 0.83 m. In [[#img-14|Figure 14]]c, the landslide mass, the sliding zone and the bedrock are presented as brown particles, green particles and blue particles, respectively.
826
827
===5.3 Mechanical parameters and micro particle parameters===
828
829
According to the report of engineering geological exploration, the unit weight and shear strengths are listed in [[#tab-4|Table 4]] for the geotechnical materials in the landslide mass, the sliding zone and the bedrock, respectively. All parameters are measured in the native state, and taken their statistical mean value. It is noticed that the stability of Baijiabao landslide is actually affected by the rainfall and the reservoir water fluctuation. In view that the main purpose of this study is to improve the ability of MatDEM in 3D slope stability evaluation, the rainfall and the reservoir water fluctuation are ignored temporarily.
830
831
<div class="center" style="font-size: 75%;">'''Table 4'''. Unit weight and shear strengths of geotechincal materials in different regions</div>
832
833
<div id='tab-4'></div>
834
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;font-size:85%;width:auto;" 
835
|-style="text-align:center"
836
! Region !! Unit weight <math>\gamma</math> (kN/m<sup>3</sup>) !! <math>C</math> (kPa) !! <math>\varphi</math> (<sup>o</sup>)
837
|-style="text-align:center"
838
|  Landslide mass
839
|  18.0
840
|  17.1
841
|  20.2
842
|-style="text-align:center"
843
|  Sliding zone
844
|  20.5
845
|  24.6
846
|  15.2
847
|-style="text-align:center"
848
|  Bedrock
849
|  26.0
850
|  1,200
851
|  32.1
852
|}
853
854
855
Using the BP neural network-based calibration method established in Section 3, micro parameters of particles representing various geotechnical materials are obtained and listed in [[#tab-5|Table 5]]. Because the particles have a much larger radius than the particles used in Section 3, the size of the virtual 3D direct shear test model is adjusted. Previous studies [21,49] verified the effect of the particle size on the macro properties in DEM modeling, and declared that the ratio <math>L</math>/<math>R</math> is a critical variable to be carefully selected when performing virtual strength test. Here <math>L</math> is the size of the strength test model, and <math>R</math> denotes the average radius of particles. According to Su’s suggestion [49], <math>L</math>/<math>R</math> should be more than 200 to avoid the effect of the particle size on the macro properties of the particle assembly in DEM. Finally, the specimen to be used in virtual 3D direct shear test has a radius of 21 m and a height of 14 m.
856
857
<div class="center" style="font-size: 75%;">'''Table 5'''. Resulted micro particle parameters and corresponding shear strangths</div>
858
859
<div id='tab-1'></div>
860
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;font-size:85%;width:auto;" 
861
|-style="text-align:center"
862
! Region !! <math>K_n</math> (N/m) !! <math>K_s</math> (N/m) !! <math>X_b</math> (m) !! <math>F_{s0}</math> (Pa) !! <math>\mu_p</math> !! <math>C</math> (kPa) !! <math>\varphi</math> (<sup>o</sup>)
863
|-style="text-align:center"
864
|  Landslide mass
865
|  <math>3.27\times 10^7 </math>
866
|  <math>5.27\times 10^6</math>
867
|  <math>1.95\times 10^{-5}</math>
868
|  <math>1.25\times 10^{5}</math>
869
|  <math>0.220</math>
870
|  <math>17.0</math>
871
|  <math>20.2</math>
872
|-style="text-align:center"
873
|  Sliding zone
874
|  <math>2.49\times 10^7 </math>
875
|  <math>1.84\times 10^6</math>
876
| <math> 2.72\times 10^{-5}</math>
877
|  <math>3.67\times 10^{5}</math>
878
|  <math>0.120</math>
879
|  <math>24.6</math>
880
|  <math>15.2</math>
881
|-style="text-align:center"
882
|  Bedrock
883
|  <math>4.27\times 10^7 </math>
884
|  <math>2.09\times 10^6</math>
885
|  <math>2.04\times 10^{-5}</math>
886
|  <math>5.95\times 10^{5}</math>
887
|  <math>0.441</math>
888
|  <math>37.5</math>
889
|  <math>32.0</math>
890
|}
891
892
893
With respect to the bedrock, the micro particle parameters in [[#tab-5|Table 5]] resulted in the shear strengths different to the actual shear strengths in [[#tab-4|Table 4]]. The reason is that it is hard to represent the three geotechnical materials by using micro parameters of the same order of magnitude. The cohesion of the bedrock is nearly 100 times that of the landslide mass and the sliding zone. Because all the particles representing the bedrock are fixed in the model, the difference between the cohesion of the bedrock in [[#tab-4|Table 4]] and that in [[#tab-5|Table 5]] will have little influence on the computation. Moreover, to reduce the workload, the micro particle parameters of the bedrock will not be changed when the stability of the landslide is evaluated by using SRM.
894
895
===5.4 Computational results===
896
897
As concluded in Section 4, the abrupt in the variation tendencies and the magnitude of the kinetic energy can be utilized to determine the FOS when estimating the slope stability by using SRM in MatDEM. Under different reduction factors, the variations of the kinetic energy with time steps are recorded for the landslide model and plotted in [[#img-15|Figure 15]]. The kinetic energy curve resulted by a greater reduction factor always has a higher situation in [[#img-15|Figure 15]]. But, if the reduction factor is less than 1.06, the difference in the magnitude of the kinetic energy between different reduction factors is relatively small at each time step. When the reduction factor is increasing from 1.06 to 1.07, the magnitude difference has a dramatically increasing. Meanwhile, the kinetic energy curves present a convergent trend when the reduction factor is less than 1.06, but present a divergent trend when the reduction factor arrives at 1.07. Thus, the FOS of Baijiabao landslide is considered as 1.06.
898
899
<div id='img-15'></div>
900
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
901
|-
902
|style="padding:10px;"| [[Image:Review_852894090255-image42.png|318px]]
903
|- style="text-align: center; font-size: 75%;"
904
| colspan="1" style="padding:10px;"| '''Figure 15'''. Variation of kinetic energy with time steps under different reduction factors
905
|}
906
907
908
By using the proposed BP neural network-based calibration method, the micro particle parameters are constantly adjusted during the computation. [[#tab-6|Table 6]] lists the micro particle parameters and the shear strengths under different reduction factors.
909
910
<div class="center" style="font-size: 75%;">'''Table 6'''. Micro particle parameters and corresponding shear strangths under different reduction factors</div>
911
912
<div id='tab-6'></div>
913
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;font-size:85%;width:auto;" 
914
|-style="text-align:center"
915
! Region !! FOS !! <math>K_n</math> (N/m)!! <math>K_s</math> (N/m) !! <math>X_b</math> (m) !!<math>F_{s0}</math> (Pa) !! <math>\mu_p</math> !! <math>C</math> (kPa) !! <math>\varphi</math> (<sup>o</sup>)
916
|-style="text-align:center"
917
|  rowspan='7' |Landslide mass
918
|  <math>1.00</math>
919
|  <math>3.27\times 10^7 </math>
920
|  <math>5.27\times 10^6</math>
921
|  <math>1.95\times 10^{-5}</math>
922
|  <math>1.25\times 10^{5}</math>
923
|  <math>0.220</math>
924
|  <math>17.0</math>
925
|  <math>20.2</math>
926
|-style="text-align:center"
927
|  <math>1.10</math>
928
|  <math>3.35\times 10^7 </math>
929
|  <math>4.89\times 10^6</math>
930
|  <math>1.61\times 10^{-5}</math>
931
|  <math>1.03\times 10^{5}</math>
932
|  <math>0.175</math>
933
|  <math>15.5</math>
934
|  <math>18.5</math>
935
|-style="text-align:center"
936
|  <math>1.03</math>
937
|  <math>3.29\times 10^7 </math>
938
|  <math>5.13\times 10^6</math>
939
|  <math>1.86\times 10^{-5}</math>
940
|  <math>1.18\times 10^{5}</math>
941
|  <math>0.207</math>
942
|  <math>16.5</math>
943
|  <math>19.7</math>
944
|-style="text-align:center"
945
|  <math>1.04</math>
946
|  <math>3.18\times 10^7 </math>
947
|  <math>4.97\times 10^6</math>
948
|  <math>1.82\times 10^{-5}</math>
949
|  <math>1.16\times 10^{5}</math>
950
|  <math>0.202</math>
951
|  <math>16.3</math>
952
|  <math>19.5</math>
953
|-style="text-align:center"
954
|  <math>1.05</math>
955
|  <math>3.22\times 10^7 </math>
956
|  <math>5.05\times 10^6</math>
957
|  <math>1.77\times 10^{-5}</math>
958
|  <math>1.14\times 10^{5}</math>
959
|  <math>0.198</math>
960
|  <math>16.2</math>
961
|  <math>19.3</math>
962
|-style="text-align:center"
963
|  <math>1.06</math>
964
|  <math>3.53\times 10^7 </math>
965
|  <math>5.02\times 10^6</math>
966
|  <math>1.75\times 10^{-5}</math>
967
|  <math>1.12\times 10^{5}</math>
968
|  <math>0.194</math>
969
|  <math>16.0</math>
970
|  <math>19.1</math>
971
|-style="text-align:center"
972
|  <math>1.07</math>
973
|  <math>3.43\times 10^7 </math>
974
|  <math>4.95\times 10^6</math>
975
|  <math>1.71\times 10^{-5}</math>
976
|  <math>1.09\times 10^{5}</math>
977
|  <math>0.190</math>
978
|  <math>15.9</math>
979
|  <math>19.0</math>
980
|-style="text-align:center"
981
|  rowspan='7' | Sliding zone
982
|  <math>1.00</math>
983
|  <math>2.49\times 10^7 </math>
984
|  <math>1.84\times 10^6</math>
985
|  <math>2.72\times 10^{-5}</math>
986
|  <math>3.67\times 10^{5}</math>
987
|  <math>0.120</math>
988
|  <math>24.6</math>
989
|  <math>15.2</math>
990
|-style="text-align:center"
991
|  <math>1.10</math>
992
|  <math>2.45\times 10^7 </math>
993
|  <math>1.81\times 10^6</math>
994
|  <math>2.36\times 10^{-5}</math>
995
|  <math>3.05\times 10^{5}</math>
996
|  <math>0.098</math>
997
|  <math>22.4</math>
998
|  <math>13.9</math>
999
|-style="text-align:center"
1000
|  <math>1.03</math>
1001
|  <math>2.69\times 10^7 </math>
1002
|  <math>1.86\times 10^6</math>
1003
|  <math>2.61\times 10^{-5}</math>
1004
|  <math>3.48\times 10^{5}</math>
1005
|  <math>0.113</math>
1006
|  <math>23.9</math>
1007
|  <math>14.8</math>
1008
|-style="text-align:center"
1009
|  <math>1.04</math>
1010
|  <math>2.38\times 10^7 </math>
1011
|  <math>1.83\times 10^6</math>
1012
|  <math>2.57\times 10^{-5}</math>
1013
|  <math>3.42\times 10^{5}</math>
1014
|  <math>0.111</math>
1015
|  <math>23.7</math>
1016
|  <math>14.6</math>
1017
|-style="text-align:center"
1018
|  <math>1.05</math>
1019
|  <math>2.33\times 10^7 </math>
1020
|  <math>1.96\times 10^6</math>
1021
|  <math>2.54\times 10^{-5}</math>
1022
|  <math>3.36\times 10^{5}</math>
1023
|  <math>0.109</math>
1024
| <math>23.4</math>
1025
|  <math>14.5</math>
1026
|-style="text-align:center"
1027
|  <math>1.06</math>
1028
|  <math>2.41\times 10^7 </math>
1029
|  <math>1.92\times 10^6</math>
1030
|  <math>2.50\times 10^{-5}</math>
1031
|  <math>3.29\times 10^{5}</math>
1032
|  <math>0.107</math>
1033
|  <math>23.2</math>
1034
|  <math>14.4</math>
1035
|-style="text-align:center"
1036
|  <math>1.07</math>
1037
|  <math>2.46\times 10^7 </math>
1038
|  <math>1.87\times 10^6</math>
1039
|  <math>2.47\times 10^{-5}</math>
1040
|  <math>3.23\times 10^{5}</math>
1041
|  <math>0.104</math>
1042
|  <math>23.0</math>
1043
|  <math>14.2</math>
1044
|}
1045
1046
1047
[[#img-16|Figure 16]] illustrates the final displacement of Baijiabao landslide model when the reduction factor is taken as 1.06. Some particles in the front area of the model have the displacements of nearly 60 m, and most particles have the displacements less than 10 m. Thus, the FOS determined by the abrupt in the variation tendencies and the kinetic energy magnitude in fact reveals the landslide stability at the global level. The displacements when the landslide is close to its limit equilibrium state may be valuable in the landslide failure prediction. The field data [46] verified that at the monitoring period from January 2007 to December 2019, the ZG326 site ([[#img-13|Figure 13]]a) in the middle area of Baijiabao landslide recorded the max cumulative displacement of 1.8 m. Its displacements are still observed when the reservoir water level changes and the rainfall happen, while Baijiabao landslide hasn’t slipped entirely till now. Therefore, the surface movements before its failure are uncertain. In views that this study aims to establish a micro particle parameter calibration method and a criterion for the application of SRM in MatDEM to analysis 3D slope stability, the influence of the reservoir water and the rainfall on the stability of Baijiabao landslide is left out of consideration. Hence, the displacements in the limit equilibrium state obtained in this study are not suitable for the failure prediction for Baijiabao landslide.
1048
1049
<div id='img-16'></div>
1050
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
1051
|-
1052
|style="padding:10px;"|  [[Image:Review_852894090255-image43-c.png|324px]]
1053
|- style="text-align: center; font-size: 75%;"
1054
| colspan="1" style="padding:10px;"| '''Figure 16'''. The displacement of Baijiabao landslide model with the reduction factor 1.06
1055
|}
1056
1057
==6. Conclusions==
1058
1059
By taking MatDEM as an example, a calibration method based on BP neural network is proposed to acquire the micro particle parameters corresponding to the given shear strengths in this paper, and an energy criterion is proposed for determining the FOS of the slope when applying SRM in DEM. And the proposed enhancements are applied in the 3D stability analysis of Baijiabao landslide. The following conclusions can be made:
1060
1061
(1) BP neural network can well represent the complicated relationship between the micro particle parameters and the shear strengths of the particle aggregate in MatDEM. And the proposed BP neural network-based micro particle parameters calibration method is applicable to provide the micro particle parameters when the shear strength is constantly reduced.
1062
1063
(2) The variations of the kinetic energy, the gravitational potential energy and the heat are noticed to be quite different between the stable and instable slope. However, the abrupt in the variation tendency and the magnitude of the kinetic energy is more suitable for determining the FOS of the slope.
1064
1065
(3) When determining the FOS by using the kinetic energy criterion, the displacements of particles in the slope model are permitted, so the FOS actually reveals the slope stability at the global level.
1066
1067
==Acknowledgments==
1068
1069
The work is supported by the Open Research Programme of the Hubei Key Laboratory of Disaster Prevention and Mitigation (China Three Gorges University) (no. 2022KJZ07), the Open Research Fund of Key Laboratory of Geological Hazards on Three Gorges Reservoir Area of Ministry of Education (2020KDZ10), National Natural Science Funds (Project No. 52079070), and the Open Research Fund of State Key Laboratory of Simulation and Regulation of Water Cycle in River Basin (China Institute of Water Resources and Hydropower Research), Grant NO:IWHR-SKL-202020.
1070
1071
==References==
1072
1073
<div class="auto" style="text-align: left;width: auto; margin-left: auto; margin-right: auto;font-size: 85%;">
1074
1075
[1] Sun G., Jiang W., Cheng S.,  Zheng H. Optimization model for determining safety factor and thrust line in landslide assessments. International Journal of Geomechanics, 17(4):04016091, 2017.
1076
1077
[2] Deng D.P., Li L.,  Zhao L.H. LEM for stability analysis of 3D slopes with general-shaped slip surfaces. International Journal of Geomechanics, 17(10):06017017, 2017.
1078
1079
[3] Firincioglu B.S.,  Ercanoglu M. Insights and perspectives into the limit equilibrium method from 2D and 3D analyses. Engineering Geology, 281:105968, 2021.
1080
1081
[4] Lv Q., Liu Y.,  Yang Q. Stability analysis of earthquake-induced rock slope based on back analysis of shear strength parameters of rock mass. Engineering Geology, 228:39-49, 2017.
1082
1083
[5] Zheng H., Sun G.,  Liu D. A practical procedure for searching critical slip surfaces of slopes based on the strength reduction technique. Computers and Geotechnics, 36(1-2):1-5, 2009.
1084
1085
[6] Espada M., Muralha J., Lemos J.V., Jiang Q., Feng X.-T., Fan Q.,  Fan Y. Safety analysis of the left bank excavation slopes of Baihetan arch dam foundation using a discrete element model. Rock Mechanics and Rock Engineering, 51(8):2597-2615, 2018.
1086
1087
[7] Lu, Y., Tan, Y. and Li, X. Stability analyses on slopes of clay-rock mixtures using discrete element method. Engineering Geology, 244: 116-124, 2018.
1088
1089
[8] Nie Z., Zhang Z., Zheng H.,  Lin S. Stability analysis of landslides using BEM and variational inequality based contact model. Computers and Geotechnics, 123:103575, 2020.
1090
1091
[9] Sun G., Lin S., Zheng H., Tan Y., Sui T. The virtual element method strength reduction technique for the stability analysis of stony soil slopes. Computers and Geotechnics, 119:103349, 2020.
1092
1093
[10] Xu D., Wu A., Yang Y., Lu B., Liu F., Zheng H. A new contact potential based three-dimensional discontinuous deformation analysis method. International Journal of Rock Mechanics and Mining Sciences, 127:104206, 2020.
1094
1095
[11] Yang Y., Xia Y., Zheng H.,  Liu Z. Investigation of rock slope stability using a 3D nonlinear strength-reduction numerical manifold method. Engineering Geology, 292:106285, 2021.
1096
1097
[12] Wang B., Vardon P.J.,  Hicks M.A. Rainfall-induced slope collapse with coupled material point method. Engineering Geology, 239:1-12, 2018.
1098
1099
[13] Bui H.H., Fukagawa R., Sako K.,  Wells J.C. Slope stability analysis and discontinuous slope failure simulation by elasto-plastic smoothed particle hydrodynamics (SPH). Géotechnique, 61(7):565-574, 2011.
1100
1101
[14] He L., Tian Q., Zhao Z., Zhao X., Zhang Q., Zhao J. Rock Slope stability and stabilization analysis using the coupled DDA and FEM method: NDDA approach. International Journal of Geomechanics, 18(6):04018044, 2018.
1102
1103
[15] Sun L., Liu Q., Abdelaziz A., Tang X.,  Grasselli G. Simulating the entire progressive failure process of rock slopes using the combined finite-discrete element method. Computers and Geotechnics, 141:104557, 2022.
1104
1105
[16] Bao Y., Sun X., Chen J., Zhang W., Han X., Zhan J. Stability assessment and dynamic analysis of a large iron mine waste dump in Panzhihua, Sichuan, China. Environmental Earth Sciences, 78(2):48, 2019.
1106
1107
[17] Jiang M., Niu M., Zhang F., Wang H., Liao Z. Instability analysis of jointed rock slope subject to rainfall using DEM strength reduction technique. European Journal of Environmental and Civil Engineering, 26(10):4664-4686, 2021.
1108
1109
[18] Shi C., Yang B., Zhang Y., Yang J. Application of discrete-element numerical simulation for calculating the stability of dangerous rock mass: A case study. International Journal of Geomechanics, 20(12):04020231, 2020.
1110
1111
[19] Zienkiewicz O.C., Humpheson C.,  Lewis R.W. Associated and non-associated visco-plasticity and plasticity in soil mechanics. Geotechnique, 25(4):671-689, 1975.
1112
1113
[20] Griffiths D.V.,  Lane P.A. Slope stability analysis by finite elements. Géotechnique, 49(3):387-403, 1999.
1114
1115
[21] Yang B., Jiao Y., Lei S. A study on the effects of microparameters on macroproperties for specimens created by bonded particles. Engineering Computations, 23(6):607-631, 2006.
1116
1117
[22] Zhou Z.H., Wang H.N., Jiang M.J. Macro- and micro-mechanical relationship of the anisotropic behaviour of a bonded ellipsoidal particle assembly in the elastic stage. Acta Geotechnica, 16(12):3899-3921, 2021.
1118
1119
[23] Cheng Y., Wong L.N.Y. A study on mechanical properties and fracturing behavior of Carrara marble with the flat‐jointed model. International Journal for Numerical and Analytical Methods in Geomechanics, 44(6):803-822, 2020.
1120
1121
[24] Liu C., Pollard D.D.,  Shi B. Analytical solutions and numerical tests of elastic and failure behaviors of close-packed lattice for brittle rocks and crystals. Journal of Geophysical Research: Solid Earth, 118(1):71-82, 2013.
1122
1123
[25] Rackl M., Hanley K.J. A methodical calibration procedure for discrete element models. Powder Technology, 307:73-83, 2017.
1124
1125
[26] Hou J., Zhang M., Chen Q., Wang D., Javadi A., Zhang S.  Failure-mode analysis of loose deposit slope in Ya’an-Kangding Expressway under seismic loading using particle flow code. Granular Matter, 21(1):8, 2019.
1126
1127
[27] Su H., Fu Z., Gao A., Wen Z. Numerical simulation of soil levee slope instability using particle-flow code method. Natural Hazards Review, 20(2):04019001, 2019.
1128
1129
[28] Xu G.J., Zhong K.Z., Fan J.W., Zhu Y.J.,  Zhang Y.Q. Stability analysis of cohesive soil embankment slope based on discrete element method. Journal of Central South University, 27(7):1981-1991, 2020.
1130
1131
[29] Wang M.,  Cao P. Calibrating the micromechanical parameters of the PFC2D (3D) models using the improved simulated annealing algorithm. Mathematical Problems in Engineering, 2017:6401835, 2017.
1132
1133
[30] Cheng H., Shuku T., Thoeni K.,  Yamamoto H. Probabilistic calibration of discrete element simulations using the sequential quasi-Monte Carlo filter. Granular Matter, 20(1):11, 2018.
1134
1135
[31] Do H.Q., Aragón A.M., Schott D.L. A calibration framework for discrete element model parameters using genetic algorithms. Advanced Powder Technology, 29(6):1393-1403, 2018.
1136
1137
[32] Benvenuti L., Kloss C.,  Pirker S. Identification of DEM simulation parameters by Artificial Neural Networks and bulk experiments. Powder Technology, 291:456-465, 2016.
1138
1139
[33] He P., Fan Y., Pan B., Zhu Y., Liu J.,  Zhu D. Calibration and verification of dynamic particle flow parameters by the back-propagation neural network based on the genetic algorithm: recycled polyurethane powder. Materials, 12(20):3350, 2019.
1140
1141
[34] Ye F., Wheeler C., Chen B., Hu J., Chen K.,  Chen W. Calibration and verification of DEM parameters for dynamic particle flow conditions using a backpropagation neural network. Advanced Powder Technology, 30(2):292-301, 2019.
1142
1143
[35] Zhou, Y., Wu S., Jiao J.,  Zhang X. Research on mesomechanical parameters of rock and soil mass based on BP neural network. Rock Soil Mech, 32(12):3821-3826, 2011.
1144
1145
[36] Wang H., Zhang B., Mei G., Xu N. A statistics-based discrete element modeling method coupled with the strength reduction method for the stability analysis of jointed rock slopes. Engineering Geology, 264:105247, 2020.
1146
1147
[37] Deluzarche R.,  Cambou B. Discrete numerical modelling of rockfill dams. International Journal for Numerical and Analytical Methods in Geomechanics, 30(11):1075-1096, 2006.
1148
1149
[38] Shen H., Abbas S.M. Rock slope reliability analysis based on distinct element method and random set theory. International Journal of Rock Mechanics and Mining Sciences, 61:15-22, 2013.
1150
1151
[39] Liu C., Xu Q., Shi B., Deng S., Zhu H. Mechanical properties and energy conversion of 3D close-packed lattice model for brittle rocks. Computers & Geosciences, 103:12-20, 2017.
1152
1153
[40] Liu Y., Zhang D., Wang G.-y., Liu C., Zhang Y. Discrete element method-based prediction of areas prone to buried hill-controlled earth fissures. Journal of Zhejiang University-SCIENCE A, 20(10):794-803, 2019.
1154
1155
[41] Luo H., Xing A., Jin K., Xu S., Zhuang Y. Discrete element modeling of the Nayong rock avalanche, Guizhou, China constrained by dynamic parameters from seismic signal inversion. Rock Mechanics and Rock Engineering, 54(4):1629-1645, 2021.
1156
1157
[42] Scaringi G., Fan X., Xu Q., Liu C., Ouyang C., Domènech G., Yang F.,  Dai L. Some considerations on the use of numerical methods to simulate past landslides and possible new failures: the case of the recent Xinmo landslide (Sichuan, China). Landslides, 15(7):1359-1375, 2018.
1158
1159
[43] Zhan Q., Wang S., Wang L., Guo F., Zhao D., Yan J. Analysis of failure models and deformation evolution process of geological hazards in Ganzhou city, China. Frontiers in Earth Science, 9:731447, 2021.
1160
1161
[44] Graziani A., Rossini C.,  Rotonda T. Characterization and DEM modeling of shear zones at a large dam foundation. International Journal of Geomechanics, 12(6):648-664, 2012.
1162
1163
[45] Peña J.M., Lozano J.A.,  Larrañaga P. An empirical comparison of four initialization methods for the K-Means algorithm. Pattern Recognition Letters, 20(10):1027-1040, 1999.
1164
1165
[46] Yao W.M., Li C.D., Guo Y.C., Criss R.E., Zuo Q.J.,  Zhan H.B. Short-term deformation characteristics, displacement prediction, and kinematic mechanism of Baijiabao landslide based on updated monitoring data. Bulletin of Engineering Geology and the Environment, 81(9):393, 2022.
1166
1167
[47] Xu S.L., Niu R.Q. Displacement prediction of Baijiabao landslide based on empirical mode decomposition and long short-term memory neural network in Three Gorges area, China. Computers & Geosciences, 111:87-96, 2018.
1168
1169
[48] Yao W.M., Li C.D., Zuo Q.J., Zhan H.B., Criss R.E. Spatiotemporal deformation characteristics and triggering factors of Baijiabao landslide in Three Gorges Reservoir region, China. Geomorphology, 343:34-47, 2019.
1170
1171
[49] Su, H., Yang J., Hu B., Gao X., Ma H. Study of particle size effect of rock model based on particle discrete element method. Rock and Soil Mechanics, 39(12):4642-4650, 2018.
1172

Return to Jiang et al 2022b.

Back to Top

Document information

Published on 19/01/23
Accepted on 13/01/23
Submitted on 03/12/22

Volume 39, Issue 1, 2023
DOI: 10.23967/j.rimni.2023.01.003
Licence: CC BY-NC-SA license

Document Score

0

Views 190
Recommendations 0

Share this document

claim authorship

Are you one of the authors of this document?