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
''', , José Manuel González, <math>\cdot </math>Francisco Zárate, <math>\cdot </math>Eugenio Oñate'''
5
-->
6
7
==Abstract==
8
9
In this paper we analyze the capabilities of two numerical techniques based on DEM and FEM-DEM approaches for the simulation of fracture in shale rock caused by a pulse of pressure. We have studied the evolution of fracture in several fracture scenarios related to the initial stress state in the soil or the pressure pulse peak. Fracture length and type of failure have been taken as reference for validating the models. The results obtained show a good approximation to FEM results from the literature.
10
11
'''Keywords''': discrete element method, finite element method, FEM-DEM, pulse fracture, shale rock.
12
13
==1 Introduction==
14
15
Fracking techniques have become increasingly popular for the exploitation of natural gas reservoirs in shale rock. The potential of these techniques is based in the efficiency of the cracking network created in the rock. The network depends on the material properties, the fracking technology employed and the pressure function applied. The most common fracking techniques involve explosives, gas or hydraulic pressure. The initial stress state, the pressure peack and the time of application of the pressure load are relevant variables that influence the final configuration of the cracks. A brittle failure is expected to reach a wide, complex and well-connected network.
16
17
The numerical simulation of fracture is a popular research area in computational mechanics. Several numerical approaches are based on the finite element method (FEM). Non-conventional FEM technique have been formulated to properly reproduce cracks in concrete and geomaterials <span id='citeF-1'></span><span id='citeF-2'></span><span id='citeF-3'></span><span id='citeF-7'></span><span id='citeF-10'></span><span id='citeF-15'></span><span id='citeF-16'></span>[[#cite-1|[1]],[[#cite-2|2]],[[#cite-3|3]],[[#cite-7|7]],[[#cite-10|10]],[[#cite-15|15]],[[#cite-16|16]]].
18
19
The discrete element method (DEM) is an alternative and efficient numerical technique to reproduce multifracture situations in a solid. The simulation of the cracking response of concrete and geomaterials have been extensively treated in the last years using the DEM <span id='citeF-5'></span><span id='citeF-8'></span><span id='citeF-9'></span><span id='citeF-13'></span><span id='citeF-14'></span><span id='citeF-17'></span><span id='citeF-18'></span><span id='citeF-20'></span><span id='citeF-23'></span><span id='citeF-24'></span><span id='citeF-26'></span><span id='citeF-27'></span>[[#cite-5|[5]],[[#cite-8|8]],[[#cite-9|9]],[[#cite-13|13]],[[#cite-14|14]],[[#cite-17|17]],[[#cite-18|18]],[[#cite-20|20]],[[#cite-23|23]],[[#cite-24|24]],[[#cite-26|26]],[[#cite-27|27]]],  or a combination of DEM and FEM procedures <span id='citeF-11'></span><span id='citeF-18'></span><span id='citeF-19'></span><span id='citeF-24'></span><span id='citeF-25'></span><span id='citeF-28'></span>[[#cite-11|[11]],[[#cite-18|18]],[[#cite-19|19]],[[#cite-24|24]],[[#cite-25|25]],[[#cite-28|28]]].
20
21
The objective of this work is to analyze the capabilities of the DEM procedure proposed by Oñate ''et al.'' <span id='citeF-17'></span>[[#cite-17|[17]]] and the FEM-DEM approach proposed by Zarate and Oñate <span id='citeF-19'></span>[[#cite-19|[19]]] to simulate the fracture behavior of a shale rock when a pulse of pressure is applied. In the DEM procedure the rock is modelled as a collection of circular particles (in 2D) or spheres (in 3D). In the FEM-DEM approach the shale rock domain is initially discretized with a finite element mesh which is progressively transformed into a collection of particles as cracking evolves. Emphasis is put in the simulation of the fracture pattern, the length of the cracks and the type of failure, as these are relevant parameters for shale gas prospection in order to optimize the shape and extend of the cracking network.
22
23
Reza ''et al.'' <span id='citeF-21'></span><span id='citeF-22'></span>[[#cite-21|[21]],[[#cite-21|22]]] have presented a rate dependent constitutive model to simulated the cracks in shale rock for different  pulse loads using the FEM. The results of this work are taken as a reference in this paper for comparing the cracking patterns and the types of failure for different initial stress states and pressure peaks using the DEM and the FEM-DEM techniques.
24
25
The layout of the paper is as follows. In section [[#2 Pulse fracturing problem|2]] the basic concept of pulse fracturing is presented. Section [[#3 DEM numerical simulation|3]] presents the kinematic and constitutive equations for the DEM. Next a 2D model for the shale rock region including the geometry, the mesh, the material properties, the loading and the boundary conditions is described. Then the results for the pulse fracture problem obtained with the DEM are presented. In the following section the formulation of the FEM-DEM technique and the results obtained with this method for the same pulse fracture problem are presented. Finally, some conclusions are outlined.
26
27
==2 Pulse fracturing problem==
28
29
The efficiency of fracking methods demands an extensive and stable fracture network to reach the smallest gas reservoir and to extract the gas from the shale formation.
30
31
The fracture network typically depends on the geotechnical properties of the soil, the in-situ stresses when fracture is produced, the texture and the natural fracture features of the rock and the available technology to apply the fracturing load usually located within the soil at different depths. There are some different techniques to fracture the soil depending on its properties.
32
33
Hydraulic fracture is applied using gas or water at high pressure under a relatively slow loading rate loading producing a one-wing or at the most bi-wing fracture scheme <span id='citeF-21'></span>[[#cite-21|[21]]]  aligned with the direction of the maximum main stress, which reduces the expansion area of the crack.
34
35
Explosives imply a very high rate of loading leading to the simultaneous propagation of multiple fractures. Due to the high stresses and heat flux, the material next to the blasting point reaches a plastic behavior. Hence, a compaction pattern with residual stresses in the soil remains after the explosion.
36
37
Pulse fracturing considers a peak load slightly exceeding the maximum and minimum in-situ stresses but lower than the plastic or compaction limit. Peak loads and loading rates can be customized according to the soil mechanical properties to reach an optimized cracking network to improve the shale gas production efficiency. Brittle fractures are optimum for that goal so the definition of the pressure pulse is compulsory.
38
39
Typically the pulse fracturing process involves two loading processes. First, the combustion of the propellant leads to a dynamic load characterized by high energy stresses and shock waves next to the wellbore. High amplitude compressive stresses waves propagate around the wellbore and within the rock mass. As consequence, the surrounding rock support tensile stresses. The result are radial cracks with different orientations.
40
41
Once the first radial cracks appear, gas penetrates into them increasing the width and depth of the existing cracks. This second phase is developed at a lower speed as a quasi-static process <span id='citeF-22'></span>[[#cite-22|[22]]]. The gas penetration can create cracks of significant length.
42
43
In our work the DEM and DEM-FEM techniques proposed in <span id='citeF-17'></span>[[#cite-17|[17]]] and <span id='citeF-19'></span>[[#cite-19|[19]]], respectively are used for modelling the initial radial cracks induced by the pressure pulse. The effect of gas penetration on the cracks is not taken into account in our work. Hence, the actual fracture pattern will exhibit extended fractures for all the cases studied. Gas penetration into the fracture network can be considered using  the computational strategy presented in <span id='citeF-28'></span>[[#cite-28|[28]]].
44
45
==3 DEM numerical simulation==
46
47
===3.1 DEM model overview===
48
49
The DEM was initially developed by Cundall ''et al.'' <span id='citeF-4'></span>[[#cite-4|[4]]] in the 1970's. It is based in the interaction of discrete elements (also called particles) &#8211; typically cylinders (in 2D) and spheres (in 3D) &#8211; to simulate the behavior of continuum and discontinuum domains. This interaction is governed by a set of kinematic equations involving the forces acting over the discrete elements and the displacements, velocities and accelerations of the elements. The forces acting over a discrete element are related to the stresses and strains according to a particular constitutive model. In our work we use the local constitutive model for the DEM for cohesive and non-cohesive materials proposed by Oñate et al. <span id='citeF-17'></span>[[#cite-17|[17]]]. In the following a brief description of this model is presented.
50
51
====3.1.1 Kinematic equations and integration scheme====
52
53
The translation and rotation of the discrete particles is governed by the standard dynamics equations for rigid bodies,
54
55
<span id="eq-1"></span>
56
{| class="formulaSCP" style="width: 100%; text-align: left;" 
57
|-
58
| 
59
{| style="text-align: left; margin:auto;width: 100%;" 
60
|-
61
| style="text-align: center;" | <math>m_i{\ddot{\boldsymbol{u}}}_i={\boldsymbol{F}}_{i} \qquad \qquad \boldsymbol{I}_i{\dot{\boldsymbol{\omega }}}_i={\boldsymbol{T}}_{{i}} </math>
62
|}
63
| style="width: 5px;text-align: right;white-space: nowrap;" | (1)
64
|}
65
66
where <math display="inline">{\boldsymbol u}_{i}</math> and <math display="inline">\boldsymbol{\omega }_{i}</math> are the <math display="inline">i</math>-th particle displacement and the angular velocity, <math display="inline">m_{i}</math> and <math display="inline">I_{i}</math> are mass and the inertia of the particle, and <math display="inline">{\boldsymbol F}_{i}</math> and <math display="inline">{\boldsymbol T}_{i}</math> are vector containing the forces and torques resultant respect to the central axes due to the interaction of a particle with its neighbors (Figure [[#img-1|1]]). 
67
68
<div id='img-1'></div>
69
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
70
|-
71
|[[Image:Draft_Samper_773076048-Fig2con327.png|350px|Motion of a rigid particle]]
72
|- style="text-align: center; font-size: 75%;"
73
| colspan="1" | '''Figure 1:''' Motion of a rigid particle
74
|}
75
76
The set of forces applied on a particle include external forces (<math display="inline">{\boldsymbol{F}}^{ext}_i</math>), damping forces (<math display="inline">{\boldsymbol{F}}^{damp}_i</math>) and interaction forces between neighbor particles (<math display="inline">{\boldsymbol{F}}^{ij}</math>) (Figure [[#img-2|2]])
77
78
<span id="eq-2"></span>
79
{| class="formulaSCP" style="width: 100%; text-align: left;" 
80
|-
81
| 
82
{| style="text-align: left; margin:auto;width: 100%;" 
83
|-
84
| style="text-align: center;" | <math>{\boldsymbol{F}}_i={\boldsymbol{F}}^{ext}_i+{\boldsymbol{F}}^{damp}_i+\sum \limits _{j=1}^{n_i}{\boldsymbol{F}}^{ij} </math>
85
|}
86
| style="width: 5px;text-align: right;white-space: nowrap;" | (2)
87
|}
88
89
where <math display="inline">n_i</math> is the number of particles that are in contact with the <math display="inline">i</math>th particle.
90
91
<div id='img-2a'></div>
92
<div id='img-2b'></div>
93
<div id='img-2'></div>
94
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
95
|-
96
|[[Image:Draft_Samper_773076048-contact-interface.png|340px|]]
97
|[[Image:Draft_Samper_773076048-Fig6con327.png|200px|]]
98
|- style="text-align: center; font-size: 75%;"
99
| (a) 
100
| (b) 
101
|- style="text-align: center; font-size: 75%;"
102
| colspan="2" | '''Figure 2:''' (a) Definition of contact interfaces. (b) Forces acting on a contact interface section <math>{A}^{ij}</math> along the normal and shear directions
103
|}
104
105
The expression for the torques can be derived from Eq.[[#eq-2|2]] <span id='citeF-18'></span>[[#cite-18|[18]]]. The dynamic equations (??) are integrated in time using an explicit scheme as shown in Eq.[[#eq-3|3]] for the translation motion
106
107
<span id="eq-3"></span>
108
{| class="formulaSCP" style="width: 100%; text-align: left;" 
109
|-
110
| 
111
{| style="text-align: left; margin:auto;width: 100%;" 
112
|-
113
| style="text-align: center;" | <math>{\ddot{\boldsymbol{u}}}^n_i=\frac{{\boldsymbol{F}}^n_i}{m_i} \qquad   \qquad         {\dot{\boldsymbol{u}}}^{n+1/2}_i={\dot{\boldsymbol{u}}}^{n-1/2}_i+{\ddot{\boldsymbol{u}}}^n_i\Delta t     \qquad  \qquad    {\boldsymbol{u}}^{n+1}_i={\boldsymbol{u}}^n_i+{\dot{\boldsymbol{u}}}^{n+1/2}_i\Delta t </math>
114
|}
115
| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
116
|}
117
118
The explicit time integration scheme is chosen due to the high computational cost of the DEM solution for large problems. However the stability of the scheme is conditioned to the time step value. The critical time step is related to the high frequency of the problem, (<math display="inline">\omega ^{max}</math>), i.e.
119
120
<span id="eq-4"></span>
121
{| class="formulaSCP" style="width: 100%; text-align: left;" 
122
|-
123
| 
124
{| style="text-align: left; margin:auto;width: 100%;" 
125
|-
126
| style="text-align: center;" | <math>\Delta t\le {\Delta t}_{cr}=\frac{2}{{\omega }^{max}}\left(\sqrt{1+{\xi }^2}-\xi \right) </math>
127
|}
128
| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
129
|}
130
131
where <math display="inline">\xi </math> is a fraction of the critical damping [ ].
132
133
====3.1.2 Forces acting over the discrete element====
134
135
The interaction forces at the contact interface between two particles <math>i</math> and <math>j</math> (<math display="inline">\mathbf{F}^{ij}</math>) are obtained from the normal (<math display="inline">\mathbf{F}_{n}^{ij}</math>) and tangential (<math display="inline">\mathbf{F}_{s}^{ij}</math>) components (Figure [[#img-2|2b]]).
136
137
The normal component of the interaction forces is calculated as,
138
139
<span id="eq-5"></span>
140
{| class="formulaSCP" style="width: 100%; text-align: left;" 
141
|-
142
| 
143
{| style="text-align: left; margin:auto;width: 100%;" 
144
|-
145
| style="text-align: center;" | <math>{{F}}^{ij}_n={\sigma }_n {\alpha }_{ij}{A}^{ij}\qquad \hbox{with}\qquad  {A}^{ij}=\pi r^2_c </math>
146
|}
147
| style="width: 5px;text-align: right;white-space: nowrap;" | (5)
148
|}
149
150
where <math display="inline">\sigma _{n}</math> is the normal stress at the contact interface, <math display="inline">r_{c}</math> is the minimum radius between the two particles (Figure [[#img-2|2]]a) and <math display="inline">{\alpha }_{ij}</math> is a parameter that depends on the number of contacts and the packing of the particles <span id='citeF-17'></span>[[#cite-17|[17]]]. In our work we have used a global definition of <math display="inline">\alpha _{ij} =\alpha=40 \frac{P}{N_c}</math>. where <math display="inline">N_c</math> and <math display="inline">P</math> are the average number of contacts per sphere and the average porosity for the whole particle assembly. The normal stress <math display="inline">\sigma _{n}</math> is calculated from the strain and the strain rate along the normal direction as,
151
152
<span id="eq-6"></span>
153
{| class="formulaSCP" style="width: 100%; text-align: left;" 
154
|-
155
| 
156
{| style="text-align: left; margin:auto;width: 100%;" 
157
|-
158
| style="text-align: center;" | <math>{\sigma }_n=E{\varepsilon }_n+c{\dot{\varepsilon }}_n </math>
159
|}
160
| style="width: 5px;text-align: right;white-space: nowrap;" | (6)
161
|}
162
163
where <math display="inline">E</math> is the Young modulus, <math display="inline">c</math> is a local damping parameter calculated as a fraction <math display="inline">\xi </math> of the critical damping <math display="inline">\bar c</math>  per unit of length, as
164
165
<span id="eq-7"></span>
166
{| class="formulaSCP" style="width: 100%; text-align: left;" 
167
|-
168
| 
169
{| style="text-align: left; margin:auto;width: 100%;" 
170
|-
171
| style="text-align: center;" | <math>c=\frac{\xi \overline{c}}{r_c}=2\frac{\xi }{r_c}\sqrt{m_{ij}K_n} ,  \qquad \hbox{with}\qquad m_{ij}=\frac{m_im_j}{m_i+m_j}   </math>
172
|}
173
| style="width: 5px;text-align: right;white-space: nowrap;" | (7)
174
|}
175
176
where <math display="inline">K_n</math> is the normal stiffness parameter (see Eq.[[#eq-9|9]]).
177
178
The normal strain and strain rate values are computed from the kinematic variables as,
179
180
<span id="eq-8"></span>
181
{| class="formulaSCP" style="width: 100%; text-align: left;" 
182
|-
183
| 
184
{| style="text-align: left; margin:auto;width: 100%;" 
185
|-
186
| style="text-align: center;" | <math>{\varepsilon }_n=\frac{{\boldsymbol{u}}_n}{d_{ij}}\qquad  {\dot{\varepsilon }}_n=\frac{{\dot{\boldsymbol{u}}}_n}{d_{ij}} </math>
187
|}
188
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
189
|}
190
191
where <math display="inline">u_n</math> and <math display="inline">\dot u_n</math> are the relative displacement and the relative velocity between two particles along the normal direction at the contact interface and <math display="inline">d_{ij}</math> is the distance between the centroids of the two particles (Figure&nbsp;[[#img-2|2]]a).
192
193
Equations [[#eq-5|5]]&#8211;[[#eq-8|8]] lead to a general relation between  the normal <math display="inline">{F}^{ij}_n</math> force and the kinematic variables (Figure&nbsp;[[#img-2|2]]b) as
194
195
<span id="eq-9"></span>
196
{| class="formulaSCP" style="width: 100%; text-align: left;" 
197
|-
198
| 
199
{| style="text-align: left; margin:auto;width: 100%;" 
200
|-
201
| style="text-align: center;" | <math>{{F}}^{ij}_n=\frac{{\alpha }_{ij}{ A}_{ij}}{d_{ij}}\left[Eu_n+2\frac{\xi }{r_c}\sqrt{m_{ij}K_n}{\dot{u}}_n\right]=K_nu_n+C_n{\dot{u}}_n </math>
202
|}
203
| style="width: 5px;text-align: right;white-space: nowrap;" | (9)
204
|}
205
206
where <math display="inline">K_{n}</math> and <math display="inline">C_{n}</math> are the normal stiffness and the normal viscous damping parameters that can be deduced from Eq.[[#eq-9|9]].
207
208
A similar approach leads to the constitutive expression for the shear forces in the two tangential directions as (Figure&nbsp;[[#img-2|2]]b) <span id='citeF-17'></span>[[#cite-17|[17]]]
209
210
<span id="eq-10"></span>
211
{| class="formulaSCP" style="width: 100%; text-align: left;" 
212
|-
213
| 
214
{| style="text-align: left; margin:auto;width: 100%;" 
215
|-
216
| style="text-align: center;" | <math>{\boldsymbol{F}}^{ij}_s= [{F}^{ij}_{s_1} ,{F}^{ij}_{s_2}]^T = {\tau }_s{\overline{A}}^{ij}=G{\gamma }_s  {\alpha }_{ij}{A}_{ij}=G {\alpha }_{ij}{A}_{ij} \frac{{\boldsymbol{u}}^{ij}_s}{d_{ij}} </math>
217
|}
218
| style="width: 5px;text-align: right;white-space: nowrap;" | (10)
219
|}
220
221
where vector <math display="inline">{\boldsymbol{u}}^{ij}_s</math> contains  the shear components of the relative displacements between particles. This vector is calculated as,
222
223
<span id="eq-11"></span>
224
{| class="formulaSCP" style="width: 100%; text-align: left;" 
225
|-
226
| 
227
{| style="text-align: left; margin:auto;width: 100%;" 
228
|-
229
| style="text-align: center;" | <math>{\boldsymbol{u}}^{ij}_s={\boldsymbol{u}}^{ij}-\left({\boldsymbol{u}}^{ij} , {\boldsymbol{n}}^{ij}\right){\boldsymbol{n}}^{ij} </math>
230
|}
231
| style="width: 5px;text-align: right;white-space: nowrap;" | (11)
232
|}
233
234
Equation [[#eq-10|10]] leads to a general expression for the shear stiffness (assumed to be the same for both shear directions), as
235
236
<span id="eq-12"></span>
237
{| class="formulaSCP" style="width: 100%; text-align: left;" 
238
|-
239
| 
240
{| style="text-align: left; margin:auto;width: 100%;" 
241
|-
242
| style="text-align: center;" | <math>K_s= G \frac{{\alpha }_{ij}{ A}_{ij}}{d_{ij}}=\frac{K_n}{2\left(1+\nu \right)} </math>
243
|}
244
| style="width: 5px;text-align: right;white-space: nowrap;" | (12)
245
|}
246
247
where <math display="inline">\nu </math> is the Poisson's ratio of the material.
248
249
The damping forces are computed from the application of a global damping over the set of particles. This damping component is characterized by translation (<math display="inline">\alpha ^{t}</math>) and rotation (<math display="inline">\alpha ^{r}</math>) damping parameters defined as a fraction of the stiffness parameters. In this work we have taken <math display="inline">\alpha ^{r} = \alpha ^{t} = 0,10</math>. The damping forces act in opposite direction to the motion of the particles according to the following expressions <span id='citeF-17'></span>[[#cite-17|[17]]]:
250
251
<span id="eq-13.a"></span>
252
<span id="eq-13.b"></span>
253
{| class="formulaSCP" style="width: 100%; text-align: left;" 
254
|-
255
| 
256
{| style="text-align: left; margin:auto;width: 100%;" 
257
|-
258
| style="text-align: center;" | <math>{\boldsymbol{F}}^{damp}_i =-{\alpha }^t\left|{\boldsymbol{F}}^{ext}_i+{\boldsymbol{F}}^{ij}\right| \frac{{\dot{\boldsymbol{u}}}_i}{\left|{\dot{\boldsymbol{u}}}_i\right|}</math>
259
| style="width: 5px;text-align: right;white-space: nowrap;" | (13.a)
260
|-
261
| style="text-align: center;" | <math> {\boldsymbol{T}}^{damp}_i =-{\alpha }^r\left|{\boldsymbol{T}}_i\right|\frac{{\dot{\boldsymbol{\omega }}}_i}{\left|{\dot{\omega }}_i\right|} </math>
262
| style="width: 5px;text-align: right;white-space: nowrap;" | (13.b)
263
|}
264
|}
265
266
More details of the local constitutive model for the DEM can be found in <span id='citeF-17'></span>[[#cite-17|[17]]].
267
268
===3.2 Normal and shear failure===
269
270
Cohesive bonds at a contact interface are assumed to start breaking when the interface strength is exceeded in the normal direction by the tensile  contact force, ''or'' in the tangential direction by the shear force. The ''uncoupled'' failure (decohesion) criterion for the normal and tangential directions at the contact interface between particles <math display="inline">i</math> and <math display="inline">j</math> is written as
271
272
<span id="eq-14"></span>
273
{| class="formulaSCP" style="width: 100%; text-align: left;" 
274
|-
275
| 
276
{| style="text-align: left; margin:auto;width: 100%;" 
277
|-
278
| style="text-align: center;" | <math>F_{n_t}  \ge {\cal F}_{n_t}\quad , \quad F_s  \ge {\cal F}_s </math>
279
|}
280
| style="width: 5px;text-align: right;white-space: nowrap;" | (14)
281
|}
282
283
where <math display="inline">{\cal F}_{n_t}</math> and <math display="inline">{\cal F}_s</math> are the interface strengths for pure tension and  shear-compression conditions, respectively, <math display="inline">F_{n_t}</math> is the normal tensile force and <math display="inline">{F}_s</math> is the modulus of the shear force vector defined in Eq.([[#eq-10|10]]).
284
285
The interface strengths are defined as
286
287
<span id="eq-15"></span>
288
{| class="formulaSCP" style="width: 100%; text-align: left;" 
289
|-
290
| 
291
{| style="text-align: left; margin:auto;width: 100%;" 
292
|-
293
| style="text-align: center;" | <math>{\cal F}_{n_t} = \sigma _t^f \bar A^{ij}  \quad ,\quad {\cal F}_s = \tau ^f \bar A^{ij} + \mu _1 \vert F_{n_c}\vert   </math>
294
|}
295
| style="width: 5px;text-align: right;white-space: nowrap;" | (15)
296
|}
297
298
where <math display="inline">\sigma _t^f</math> and <math display="inline">\tau ^f</math> are the tensile and shear strengths respectively, <math display="inline">F_{n_c}</math> is the compressive normal force at the contact interface and <math display="inline">\mu _1=\tan \phi _1</math> is a (static) friction parameter, where <math display="inline">\phi _1</math> is an internal friction angle. These values are assumed to be an intrinsic property of the material and are determined experimentally. In our work <math display="inline">\sigma _t^f</math> is  taken as the tensile strength of the material measured in a bending-tensile (BT) test (i.e. <math display="inline">\sigma ^f_t = (\sigma ^f_t)_{BT}</math>).
299
300
As for the shear strength <math display="inline">\tau ^f</math> we have estimated its value as a percentage of the maximum compressive stress in an UCS test, <math display="inline">(\sigma _{n_c}^f)_{UCS}</math>, as
301
302
<span id="eq-16"></span>
303
{| class="formulaSCP" style="width: 100%; text-align: left;" 
304
|-
305
| 
306
{| style="text-align: left; margin:auto;width: 100%;" 
307
|-
308
| style="text-align: center;" | <math>\tau ^f = \beta (\sigma _{n_c}^f)_{UCS}  </math>
309
|}
310
| style="width: 5px;text-align: right;white-space: nowrap;" | (16)
311
|}
312
313
where <math display="inline">\beta </math> is a parameter that is calibrated in numerical experiments via UCS and BTS tests. For shale rock materials we have used <math display="inline">\beta =0.50</math>.
314
315
Following tension failure, the constitutive behavior in the shear direction is governed by the standard Coulomb law
316
317
<span id="eq-17"></span>
318
{| class="formulaSCP" style="width: 100%; text-align: left;" 
319
|-
320
| 
321
{| style="text-align: left; margin:auto;width: 100%;" 
322
|-
323
| style="text-align: center;" | <math>{\boldsymbol F}_s= \mu _2 \vert F_{n_c}\vert \frac{{\boldsymbol u}_s}{|{\boldsymbol u}_s|}\quad \hbox{with} \quad \mu _2 = \tan \phi _2  </math>
324
|}
325
| style="width: 5px;text-align: right;white-space: nowrap;" | (17)
326
|}
327
328
where <math display="inline">\mu _2</math> is a dynamic Coulomb friction coefficient and <math display="inline">\phi _2</math> is the post-failure internal friction angle. The value of <math display="inline">\phi _2</math> is determined from experimental tests.
329
330
Figure [[#img-3|3]]a shows the graphical representation of the failure criterium described by Eqs.([[#eq-14|14]]), ([[#eq-15|15]])  and ([[#eq-17|17]]). This criterium assumes that the tension and shear forces contribute to the failure of the contact interface in a decoupled manner. On the other hand, shear failure under normal compressive forces follows a failure line that is a function of the shear failure stress, the compression force and the internal friction angle.
331
332
Elastic damage under tensile and shear conditions has been taken into account in this work by assuming a linear softening behaviour defined by the softening moduli <math display="inline">H_n</math> and <math display="inline">H_t</math> introduced into the force-displacement relationships in the normal (tensile) and shear directions, respectively. Figure [[#img-4|4]] shows the evolution of the normal tension force <math display="inline">F_{n_t}</math> and the shear force <math display="inline">F_s</math> at a contact interface until failure in terms of the relative normal and tangential displacements. The effect of damage in the two constitutive laws is also shown. For details see <span id='citeF-17'></span>[[#cite-17|[17]]].
333
334
<div id='img-3'></div>
335
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
336
|-
337
|[[Image:Draft_Samper_773076048-Fig8a.png|250px|Failure line in terms of normal and shear forces. Uncoupled failure model]]
338
|- style="text-align: center; font-size: 75%;"
339
| colspan="1" | '''Figure 3:''' Failure line in terms of normal and shear forces. Uncoupled failure model
340
|}
341
342
<div id='img-4'></div>
343
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
344
|-
345
|[[Image:Draft_Samper_773076048-Fig8b.png|500px|Undamaged and damaged elastic moduli under  tension (a) and  shear (b) forces]]
346
|- style="text-align: center; font-size: 75%;"
347
| colspan="1" | '''Figure 4:''' Undamaged and damaged elastic moduli under  tension (a) and  shear (b) forces
348
|}
349
350
===3.3 Elasto-plastic model for compression forces===
351
352
The compressive stress-strain behaviour in the normal direction at the contact interface for frictional cohesive materials, such as cement, rock and concrete, is typically governed by an initial elastic law followed by a non-linear constitutive equation that varies for each material. The compressive normal stress increases under  linear elastic conditions until it reaches the  limit normal compressive  stress <math display="inline">\sigma _{n_c}^l</math> (also called yield stress). This is defined as the axial stress level where the experimental curve relating the axial stress and the axial strain   starts to deviate from the linear elastic behaviour. After this point the material is assumed to yield under ''elastic-plastic'' conditions.
353
354
The ''elasto-plastic relationships'' in the normal compressive direction are defined as<br/>
355
356
'''Loading path'''
357
<span id="eq-18"></span>
358
359
<span id="eq-18.a"></span>
360
{| class="formulaSCP" style="width: 100%; text-align: left;" 
361
|-
362
| 
363
{| style="text-align: left; margin:auto;width: 100%;" 
364
|-
365
| style="text-align: center;" | <math>dF_{n_c}=K_{T_n} du_n  </math>
366
|}
367
| style="width: 5px;text-align: right;white-space: nowrap;" | (18.a)
368
|}
369
370
'''Unloading path'''
371
372
<span id="eq-18.b"></span>
373
{| class="formulaSCP" style="width: 100%; text-align: left;" 
374
|-
375
| 
376
{| style="text-align: left; margin:auto;width: 100%;" 
377
|-
378
| style="text-align: center;" | <math>dF_{n_c}=K_{n_0}  du_n  </math>
379
|}
380
| style="width: 5px;text-align: right;white-space: nowrap;" | (18.b)
381
|}
382
383
In Eqs.([[#eq-18|18]]) <math display="inline">dF_{n_c}</math> and <math display="inline">du_n</math> are respectively the increment of the normal compressive force and the normal (relative) displacement, <math display="inline"> K_{n_0}</math> is the initial (elastic) compressive stiffness for a value of <math display="inline">E=E_0</math> (Figure [[#img-5|5]]), and <math display="inline">K_{T_n}</math> is the tangent compressive stiffness given by
384
385
<span id="eq-19"></span>
386
{| class="formulaSCP" style="width: 100%; text-align: left;" 
387
|-
388
| 
389
{| style="text-align: left; margin:auto;width: 100%;" 
390
|-
391
| style="text-align: center;" | <math>K_{T_n}=\frac{E_T}{E_0}K_{n_0}  </math>
392
|}
393
| style="width: 5px;text-align: right;white-space: nowrap;" | (19)
394
|}
395
396
where <math display="inline">E_T</math> is the slope of the normal stress-strain curve in the elastoplastic branch (i.e. <math display="inline">K_T=E_1</math>, <math display="inline">E_2</math>, <math display="inline">E_3</math> in Figure [[#img-5|5]]).
397
398
Plasticity effects in the normal compressive direction also affect the evolution of the tangential forces at the interface, as the interface shear strength is related to the normal compression force by Eq.([[#eq-15|15]]).
399
400
Figure [[#img-5|5]] shows the diagram relating the compressive axial stress and the compressive axial strain used for modelling the elasto-plastic constitutive behaviour of shale rock at the contact interfaces. The form of each diagram is  obtained from experimental tests on cylindrical samples.
401
402
<div id='img-5'></div>
403
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
404
|-
405
|[[Image:Draft_Samper_773076048-Fig10b.png|350px|Compressive axial stress-compressive axial strain diagram for elastoplastic material. LCS1 =σ<sub>n<sub>c</sub></sub><sup>l</sup>: limit  compressive stress defining the onset of elastoplastic behaviour at the contact interface]]
406
|- style="text-align: center; font-size: 75%;"
407
| colspan="1" | '''Figure 5:''' Compressive axial stress-compressive axial strain diagram for elastoplastic material. LCS1 <math>=\sigma _{n_c}^l</math>: limit  compressive stress defining the onset of elastoplastic behaviour at the contact interface
408
|}
409
410
===3.4 Calibration of the material parameters===
411
412
Oñate ''et al.'' <span id='citeF-17'></span>[[#cite-17|[17]]] have proposed a methodology to calibrate the material parameters for shale rock for the local DEM model presented in the previous section, using experimental results from material tests. Table [[#table-1|1]] summarizes the DEM constitutive paramaters for the analysis of shale rock in this work. For more details see <span id='citeF-17'></span>[[#cite-17|[17]]].
413
414
415
{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
416
|+ style="font-size: 75%;" |<span id='table-1'></span>Table. 1 DEM constitutive parameters analysis of  shale rock
417
418
|- style="border-top: 2px solid;"
419
| <math display="inline">\rho </math> (g/cc)
420
| <math>\mu _1</math>
421
| <math>\mu _2</math>
422
| <math>E_0</math>(GPa)
423
| <math>\nu </math>
424
| <math>\sigma ^f_{t}</math>(MPa) 
425
| <math>\tau _f</math>(MPa) 
426
427
|- style="border-top: 2px solid;"
428
|  2.55 
429
| 0.7 
430
| 0.6 
431
| 30 
432
| 0.2 
433
| 5.0 
434
| 25
435
436
|- style="border-top: 2px solid;"
437
|  LCS1 
438
| LCS2 
439
| LCS3
440
| YRC1 
441
| YRC2 
442
| YRC3 
443
| <math>\alpha </math>
444
|-
445
| (MPa) 
446
| (MPa) 
447
| (MPa) 
448
| 
449
| 
450
| 
451
| 
452
|- style="border-top: 2px solid;border-bottom: 2px solid;"
453
| 20 
454
| 30 
455
| 40 
456
| 2
457
| 5 
458
| 10
459
| 1.0
460
461
|}
462
463
The DEM formulation and the constitutive model presented in the previous lines have been implemented in the DEMPack code (www.cimne.com/dempack). The geometry, the discretization mesh and the postprocessing of the results have been performed with the pre-postprocessor GiD <span id='citeF-6'></span>[[#cite-6|[6]]]. Both software codes have been developed at CIMNE.
464
465
===3.5 UCS and BTS tests on shale rock material===
466
467
We have simulated with the DEM  a UCS test and a BTS test on a shale rock material corresponding to a Middle Brown gaseous shale in Devonian formation from Lincoln County, West Virginia. The essential material parameters for the DEM simulations were taken from <span id='citeF-21'></span>[[#cite-21|[21]]].
468
469
The simulations were carried out in a cylindrical sample of dimensions <math display="inline">150\times 300</math> mm. The material properties used for the DEM analysis  are listed in Table [[#table-1|1]]. The value of <math display="inline">\tau _f =25</math> MPa was obtained using Eq.([[#eq-16|16]]) with <math display="inline">\beta =0.5</math> and <math display="inline">(\sigma _{n_f}^f)_{UCS}=50</math> MPa as reported in <span id='citeF-21'></span>[[#cite-21|[21]]].
470
471
The curve in Figure [[#img-6|6]]  shows the axial stress-axial strain curve for the UCS test. A maximum compressive stress of 48 MPa was obtained using a discretization of 37000 spheres. This yields a 4% error versus the experimental value of 50 MPa <span id='citeF-21'></span>[[#cite-21|[21]]].
472
473
<div id='img-6'></div>
474
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
475
|-
476
|[[Image:Draft_Samper_773076048-37k_UCS_rock.png|450px|Axial stress-axial strain curve for  UCS test in shale rock material. DEM results using 37000 spheres]]
477
|- style="text-align: center; font-size: 75%;"
478
| colspan="1" | '''Figure 6:''' Axial stress-axial strain curve for  UCS test in shale rock material. DEM results using 37000 spheres
479
|}
480
481
The curve in Figure [[#img-7|7]] shows the tensile stress - versus time for the BTS test obtained with the DEM using 27000 spheres. A failure tensile stress of <math display="inline">\sigma ^f_{t}=5.4</math>  MPa was obtained. This gives a 8% error versus the experimental value of <math display="inline">\sigma ^f_{t}=5</math> MPa.
482
483
<div id='img-7'></div>
484
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
485
|-
486
|[[Image:Draft_Samper_773076048-Fig28.png|450px|Tensile stress-time curve for BTS test in shale rock material. DEM results using 27000 spheres]]
487
|- style="text-align: center; font-size: 75%;"
488
| colspan="1" | '''Figure 7:''' Tensile stress-time curve for BTS test in shale rock material. DEM results using 27000 spheres
489
|}
490
491
===3.6 DEM simulation of pulse fracture of a shale rock mass===
492
493
We present the application of the DEM technique previously described to the fracture of a shale mass under a pulse load.
494
495
====3.6.1 Geometry and mesh====
496
497
The reference problem <span id='citeF-21'></span>[[#cite-21|[21]]]  considers a square geometry 200 ft per side modelled assuming horizontal and vertical symmetry. The geometry in the simulation is a 100 ft (30,48m) sided square, with the wellbore in the lower left corner.
498
499
In this work the pulse fracture problem has been simulated using a two-dimensional (2D) plain strain square geometry of 15x8 m, assuming horizontal symmetry.
500
501
The analysis domain has been discretized with the DEM using 882693 circular discrete elements of 25 mm radius (average).
502
503
A wellbore with 152,4 mm diameter (6 inch) has been introduced in the center of the lower part of the square domain where the pulse load is applied.
504
505
The boundary conditions applied consider that no displacement is observed far away enough from the wellbore. Hence, prescribed zero displacements are applied at every boundary.
506
507
The point of application of the pulse is located deep into the ground where a pre-existing stress state exists, depending on the depth. This stress state becomes an initial stress field for the pulse fracturing problem and is a relevant variable in this study to analyze the different cracking patterns. The simulation is then split into two phases:
508
509
*  A velocity is applied at the boundaries in order to reproduce the initial stress state. The velocity is approximated so that the stress distribution reaches the expected value and shows a high degree of uniformity. The duration of this first pulse is 100 ms in order to ensure the stability of the results. This first phase takes most of the calculation time due to the dynamic response of the problem.
510
511
*  Once the initial stress state is reached, zero displacements are prescribed at all boundaries and then, the pulse load is applied. The calculation time includes a time interval after the pulse to reach the stationary state after the pulse.
512
513
The analysis domain is shown in Figure [[#img-8|8]] with the location of the wellbore, the symmetry axis and representation of the initial stress field prescribed at the boundaries.
514
515
<div id='img-8'></div>
516
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
517
|-
518
|[[Image:Draft_Samper_773076048-geometry.png|400px|Numerical model. Geometry and applied loads]]
519
|- style="text-align: center; font-size: 75%;"
520
| colspan="1" | '''Figure 8:''' Numerical model. Geometry and applied loads
521
|}
522
523
====3.6.2 Pressure pulse application====
524
525
The pulse loading curve is characterized by a linear rising branch up to the peak load, and an exponential decay (Figure [[#img-9|9]]). The  expression for the pulse load is
526
527
<span id="eq-20"></span>
528
{| class="formulaSCP" style="width: 100%; text-align: left;" 
529
|-
530
| 
531
{| style="text-align: left; margin:auto;width: 100%;" 
532
|-
533
| style="text-align: center;" | <math>P\left(t\right)=\left\{\begin{array}{ll}\beta t, & 0\le t\le t_{peak} \\ P_{peak}e^{-D\left(t-t_{peak}\right)}, & t_{peak}\le t<\infty  \end{array} \right. </math>
534
|}
535
| style="width: 5px;text-align: right;white-space: nowrap;" | (20)
536
|}
537
538
where <math display="inline">\beta </math> is the rising loading rate, <math display="inline">P_{peak}</math> is the maximum loading pressure and <math display="inline">t_{peak}</math> is the time for the maximum value of the loading pressure (<math display="inline">P_{peak}</math>). The decay exponent <math>D</math> is taken as constant and equal to <math display="inline">D= 0.5 \hbox{ms}^{-1}</math>.
539
540
Two pulse loading curves have been considered. The reference pulse loading has a peak of 100 MPa at a time of 0.50 ms; this means a loading rate (<math display="inline">\beta </math>) of 200 MPa/ms. This reference pulse is applied for different initial stress states of the soil corresponding to different depths.
541
542
<div id='img-9'></div>
543
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
544
|-
545
|[[Image:Draft_Samper_773076048-loading-pulse.png|450px|Loading pulse pressure function]]
546
|- style="text-align: center; font-size: 75%;"
547
| colspan="1" | '''Figure 9:''' Loading pulse pressure function
548
|}
549
550
To study the capabilities of the model, a variation of the pulse load is considered by increasing the pressure peak to 250 MPa, keeping the same loading rate. The load curve parameters are shown in Table [[#table-2|2]].
551
552
553
{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
554
|+ style="font-size: 75%;" |<span id='table-2'></span>Table. 2 Characterization of the pulse loading curves
555
556
|- style="border-top: 2px solid;"
557
| style="text-align: left;" |  '''Pulse''' 
558
| '''Peak'''  (MPa) 
559
| '''Loading rate''' (MPa/ms) 
560
| '''Peak time''' (ms) 
561
562
563
|- style="border-top: 2px solid;"
564
| style="text-align: left;" |  Reference 
565
| 100 
566
| 200 
567
| 0.50 
568
|-
569
| style="text-align: left;" | Peak 250 
570
| 250 
571
| 200 
572
| 1.25 
573
574
575
|}
576
577
The pressure load is modeled by applying a point load over the particles located in the wellbore bound. This load is calculated taking into account the pressure value and the area of the particles affected by the load. The force is applied to each particle in the radial direction with respect to the wellbore center.
578
579
====3.6.3 Initial stress state====
580
581
Three different samples of shale rock at different depths are considered. The depth of these samples determines the initial stress state when the pressure pulse takes place.
582
583
Taking into account the horizontal stress and the pore pressure gradient considered in <span id='citeF-21'></span>[[#cite-21|[21]]] the initial horizontal stresses are calculated and details are given in Table [[#table-3|3]]. The initial stresses are assumed to be constant over the analysis domain. Three different depths have been considered leading to the same number of simulation cases. An additional case is considered for a sample extracted from the surface, i. e. with zero initial stress. This case shows the top value of the cracks and is taken as a reference.
584
585
The study includes the influence of the anisotropy in the cracking response of the shale material, because of its importance in the fracture pattern according to <span id='citeF-21'></span>[[#cite-21|[21]]].
586
587
For all the cases mentioned before, a reference pulse loading is applied, with a load peak of 100 MPa at a loading rate of 200 MPa/ms.
588
589
590
{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
591
|+ style="font-size: 75%;" |<span id='table-3'></span>Table. 3 Initial stresses in the shale rock mass for different depths
592
593
|- style="border-top: 2px solid;"
594
| style="text-align: left;" |  '''Depth''' (ft) 
595
| <math>\sigma _{x}</math>(MPa) 
596
| <math>\sigma _{y}</math>(MPa) 
597
598
599
|- style="border-top: 2px solid;"
600
| style="text-align: left;" |  0 
601
| - 
602
| - 
603
|-
604
| style="text-align: left;" | 500 
605
| 1.0 
606
| 1.2 
607
|-
608
| style="text-align: left;" | 5 000 
609
| 10.0 
610
| 12.0 
611
|-
612
| style="text-align: left;" | 10 000 
613
| 20.0 
614
| 24.0 
615
616
617
|}
618
619
===3.7 DEM results===
620
621
In this section simulation results for the pulse fracturing problem using the DEM formulation described earlier are presented. Fracture simulation results respect to depth of the sample, that is, the initial stress state in the soil when the pulse is applied, are presented in Figures [[#img-10|10]] to [[#img-13|13]]. They show the fracture pattern and the type of fracture with the numerical results taken as reference, which are depicted in the same figures.
622
623
The fracture pattern is represented as a damage distribution around the wellbore after the pulse loading has been applied. The damage parameter takes into account the number of broken bonds over the initial number for every single element, and so, it varies from zero to one for every particle.
624
625
The fracture length is measured as the line distance between the farther damaged points from the wellbore for a defined fracture. Due to the symmetry of the numerical simulation, the lengths depicted in the figures are half of the DEM results. The  fracture length computed with the DEM is given in Table [[#table-4|4]]  and compared to the FEM results presented in <span id='citeF-21'></span>[[#cite-21|[21]]]. DEM and FEM results agree reasonably well for the cases considered.
626
627
628
{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
629
|+ style="font-size: 75%;" |<span id='table-4'></span>Table. 4 Fracture length. DEM results compared to those in <span id='citeF-21'></span>[[#cite-21|[21]]]
630
631
|- style="border-top: 2px solid;"
632
|  Depth  
633
| <math>\sigma _{x}</math>
634
| <math>\sigma _{y}</math>
635
| Pressure peak  
636
| Fracture length (ft) 
637
| Fracture length (ft)
638
| Figure
639
|-
640
| (ft) 
641
| (MPa)
642
| (MPa)
643
| (MPa)
644
| FEM result <span id='citeF-21'></span>[[#cite-21|[21]]]
645
| DEM result 
646
| 
647
648
649
|- style="border-top: 2px solid;"
650
|   0 
651
| 0 
652
| 0 
653
| 100 
654
| &#8211; 
655
| 40
656
| 10
657
|-
658
| 500  
659
| 1.0 
660
| 1.2 
661
| 100 
662
| 40 
663
| 33
664
| 11  
665
|-
666
| 5000  
667
| 10 
668
| 12 
669
| 100 
670
| 5 
671
| 5.4 
672
| 12
673
|-
674
| 10000 
675
| 20 
676
| 24 
677
| 100 
678
| 2.3 
679
| 3.2 
680
| 13
681
|-
682
| 5000 (weak anisotropy) 
683
| 5 
684
| 12 
685
| 100 
686
| &#8211; 
687
| 10.5 
688
| 14a
689
|-
690
| 5000 (strong anisotropy) 
691
| 2.5 
692
| 12 
693
| 100 
694
| &#8211; 
695
| 15,2 
696
| 14b
697
|-
698
| 5000  
699
| 10 
700
| 12 
701
| 250 
702
| 28 
703
| 30.2 
704
| 15
705
706
707
|}
708
709
The type of fracture is represented in the figures for every single case to highlight the difference between tensile (green) and shear (red) failure, due to the importance of this result. Reference results <span id='citeF-21'></span>[[#cite-21|[21]]] are depicted over the DEM simulations, so as to ease the comparison.
710
711
====3.7.1 Fracture vs initial stress state====
712
713
The fracture length in the DEM simulation shows a branching distribution similar to the reference results <span id='citeF-21'></span>[[#cite-21|[21]]] although the fracture is not so clearly defined because of the DEM discretization. The fracture lengths are comparable as well as the type of failure, mainly due to tensile failure in both cases.
714
715
The results with zero initial stress state &#8211; depth of 0ft (Figure [[#img-10|10]]) show a similar fracture distribution and pattern to the results for 500 ft (Figure [[#img-11|11]]).
716
717
<div id='img-10a'></div>
718
<div id='img-10b'></div>
719
<div id='img-10'></div>
720
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
721
|-
722
|[[Image:Draft_Samper_773076048-damage-no-initial-1.png|400px|]]
723
|- style="text-align: center; font-size: 75%;"
724
| (a)
725
|-
726
|[[Image:Draft_Samper_773076048-damage-no-initial-2.png|400px|]]
727
|- style="text-align: center; font-size: 75%;"
728
| (b)
729
|-
730
|- style="text-align: center; font-size: 75%;"
731
| colspan="1" | '''Figure 10:''' Damage distribution (a) and type of fracture (b). No initial stress. Peak pressure load: 100 MPa
732
|}
733
734
<div id='img-11a'></div>
735
<div id='img-11b'></div>
736
<div id='img-11c'></div>
737
<div id='img-11'></div>
738
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
739
|-
740
| colspan="2"|[[Image:Draft_Samper_773076048-Fig11b_1.png|450px|]]
741
|- style="text-align: center; font-size: 75%;"
742
| colspan="2" | (a)
743
|-
744
|[[Image:Draft_Samper_773076048-Fig11b_2.png|450px|]]
745
|[[Image:Draft_Samper_773076048-Fig11b_3.png|450px|]] 
746
|- style="text-align: center; font-size: 75%;"
747
| (b) 
748
| (c)
749
|- style="text-align: center; font-size: 75%;"
750
| colspan="2" |'''Figure 11:''' Damage distribution (a) and type of fracture (b and c). Depth 500 ft. <math>\sigma _{x} = 1.0</math> MPa, <math>\sigma _{y} = 1.2 </math> MPa. Peak pressure load: 100 MPa. FEM results in box from <span id='citeF-21'></span>[[#cite-21|[21]]]
751
|}
752
753
The fractured area decreases as the initial stresses increase their value (Figures [[#img-11|11]]&#8211;&#8211;[[#img-13|13]]). The fracture length is comparable to the reference results <span id='citeF-21'></span>[[#cite-21|[21]]]. Side fractures appear although they are not so well defined so as to be compared with the reference results.
754
755
<div id='img-12a'></div>
756
<div id='img-12b'></div>
757
<div id='img-12c'></div>
758
<div id='img-12'></div>
759
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
760
|-
761
|[[Image:Draft_Samper_773076048-damage-depth5000-1.png|600px|]]
762
|[[Image:Draft_Samper_773076048-Fig12b_1.png|450px|]]
763
|- style="text-align: center; font-size: 75%;"
764
| (a) 
765
| (b) 
766
|-
767
| colspan="2"|[[Image:Draft_Samper_773076048-Fig12b_2.png|450px|]]
768
|- style="text-align: center; font-size: 75%;"
769
|  colspan="2" | (c) 
770
|- style="text-align: center; font-size: 75%;"
771
| colspan="2" | '''Figure 12:''' Damage distribution (a) and (b) and type of fracture (c). Depth 5000 ft. <math>\sigma _{x} = 10</math> MPa, <math>\sigma _{y} = 12</math> MPa. Peak pressure load: 100 MPa. FEM results in box from <span id='citeF-21'></span>[[#cite-21|[21]]]
772
|}
773
774
The type of fracture is mainly due to tensile forces in both cases, although in the DEM simulation a small fracture area due to shear effects appears.
775
776
<div id='img-13a'></div>
777
<div id='img-13b'></div>
778
<div id='img-13c'></div>
779
<div id='img-13'></div>
780
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
781
|-
782
|[[Image:Draft_Samper_773076048-damage-depth10000-1.png|600px|]]
783
|[[Image:Draft_Samper_773076048-Fig13b_1.png|468px|]]
784
|- style="text-align: center; font-size: 75%;"
785
| (a) 
786
| (b) 
787
|-
788
| colspan="2"|[[Image:Draft_Samper_773076048-Fig13b_2.png|468px|]]
789
|- style="text-align: center; font-size: 75%;"
790
|  colspan="2" | (c) 
791
|- style="text-align: center; font-size: 75%;"
792
| colspan="2" | '''Figure 13:''' Damage distribution (a) and (b) and type of fracture (c). Depth 10000 ft. <math>\sigma _{x} = 20</math> MPa, <math>\sigma _{y} = 24</math> MPa. Peak pressure load: 100 MPa. FEM results in box from <span id='citeF-21'></span>[[#cite-21|[21]]]
793
|}
794
795
The high initial stress level leads to a small fracture in the DEM simulation.  Length values are comparable to those in <span id='citeF-21'></span>[[#cite-21|[21]]], as well as the type of failure.
796
797
====3.7.2 Fracture vs anisotropy====
798
799
A further study of the cracking response has been carried out to verify the assumptions pointed out in <span id='citeF-21'></span>[[#cite-21|[21]]] in relation to the degree of anisotropy of the initial stress field. Fractures tend to follow the direction of the maximum stress at each point and this behavior increases with the level of stress anisotropy  in the soil.
800
801
Two scenarios have been studied to analyze the capability of the DEM to reproduce this phenomenon. The degree of anisotropy is introduced by decreasing the minimum horizontal stress (<math display="inline">\sigma _{x}</math>), while the maximum stress (<math display="inline">\sigma _{y}</math>) remains unchanged taking as reference the sample at 5000 ft. depth and the reference pulse loading. These simulation cases are not studied in <span id='citeF-21'></span>[[#cite-21|[21]]], so there are no reference values for these results.
802
803
Two stress fields have been considered &#8211; weak and strong anisotropy &#8211; by decreasing the minimum horizontal stress to a 50% and 25% respectively, of its reference value. The initial stresses considered in these two scenarios are shown in Table [[#table-4|4]].
804
805
Figure [[#img-14|14]] shows the results obtained for the damage distribution for both cases. The results show a clear trend for the fractures that follow a predominant vertical direction coincident with the maximum stress <math display="inline">\sigma _{y}</math>.   The fracture length increases with the degree of anisotropy, as expected. The results evidences show that the DEM has a remarkable capability to reproduce the anisotropic effects in the fracture pattern for a pulse load.
806
807
<div id='img-14a'></div>
808
<div id='img-14b'></div>
809
<div id='img-14'></div>
810
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
811
|-
812
|[[Image:Draft_Samper_773076048-Fig14a.png|450px|Weak anisotropy σₓ = 5 MPa, σ<sub>y</sub> =12 MPa]]
813
|[[Image:Draft_Samper_773076048-Fig14b.png|450px|Strong anisotropy σₓ = 2.5 MPa, σ<sub>y</sub> =12  MPa]]
814
|- style="text-align: center; font-size: 75%;"
815
| (a) Weak anisotropy <math display="inline">\sigma _{x} = 5</math> MPa, <math display="inline">\sigma _{y} =12</math> MPa
816
| (b) Strong anisotropy <math display="inline">\sigma _{x} = 2.5</math> MPa, <math display="inline">\sigma _{y} =12 </math> MPa
817
|- style="text-align: center; font-size: 75%;"
818
| colspan="2" | '''Figure 14:''' Damage distribution for weak (a) and strong (b) anisotropy (Table [[#table-4|4]]).  Depth 5000 ft. Peak pressure load: 100 MPa
819
|}
820
821
===3.8 Fracture vs peak load===
822
823
Finally, as the third significant variable, the effect of the pulse peak loading on the fracture pattern is studied, keeping the loading rate constant. Fracture patterns have been obtained for the depth of 5000 ft and two loading peaks: 100 and 250 MPa. Results are depicted in Figures [[#img-12|12]] and [[#img-15|15]] including the reference results <span id='citeF-21'></span>[[#cite-21|[21]]].
824
825
The higher value of the peak pressure load the larger is the fracture zone and so is the fracture length. An increase of the shear failure area is appreciated both in the DEM results and the reference work (Figure [[#img-15|15]]). The fracture length obtained in this work is in agreement with that in <span id='citeF-21'></span>[[#cite-21|[21]]].
826
827
<div id='img-15a'></div>
828
<div id='img-15b'></div>
829
<div id='img-15c'></div>
830
<div id='img-15'></div>
831
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
832
|-
833
|[[Image:Draft_Samper_773076048-Fig15a.png|450px|]]
834
|[[Image:Draft_Samper_773076048-Fig15b_1.png|450px|]]
835
|- style="text-align: center; font-size: 75%;"
836
| (a) 
837
| (b) 
838
|-
839
| colspan="2"|[[Image:Draft_Samper_773076048-Fig15b_2.png|450px|]]
840
|- style="text-align: center; font-size: 75%;"
841
|  colspan="2" | (c) 
842
|- style="text-align: center; font-size: 75%;"
843
| colspan="2" | '''Figure 15:''' Damage distribution (a) and failure criteria (b). Depth 5000 ft. <math>\sigma _{x} = 10</math> MPa, <math>\sigma _{y} = 12</math> MPa. Peak pressure load: 250 MPa. FEM results in box from <span id='citeF-21'></span>[[#cite-21|[21]]]
844
|}
845
846
==5 FEM-DEM numerical simulation==
847
848
===5.1 Constitutive model overview===
849
850
The DEM is a flexible method to simulate non-continuum media, in particular the propagation of initial cracks. The contact properties between particles are a powerful degree of freedom to control the onset and propagation of cracks.
851
852
On the other hand, these contact properties are defined in the micro scale, while the material properties usually refer to experimental results in the macro scale. The step between both scales is not easy and needs a calibration task. The FEM otherwise is based in a continuum formulation involving the macro properties of the material. This feature avoids the calibration task and leads to the current stress state at any calculation point. For that reason the FEM allows one to establish failure criteria compatible with the energy balance equations, which makes it consistent and easy to apply for different materials.
853
854
The distance between DEM and FEM approaches is wide. Extensive research has been carried out in last years <span id='citeF-11'></span><span id='citeF-18'></span><span id='citeF-24'></span>[[#cite-11|[11,18,24]]] to combine the FEM and the DEM procedures, taking profit from the advantages of both numerical methods. Zárate and Oñate <span id='citeF-19'></span>[[#cite-19|[19]]]  have recently presented a coupled FEM-DEM formulation for the numerical simulation of cracks starting from a continuum discretization of the domain. This formulation is used here for simulating the fracture of shale rock under a pulse load.
855
856
In the 2D FEM-DEM formulation considered here the continuum is initially discretized using linear 3-noded triangles whose nodes define the eventual future location of a discrete element. The DEM only takes part in the simulation process when cracks appear. The normal contact forces between discrete elements are calculated by integrating the stiffness matrix of the linear triangle along its sides. Each side connects the discrete particles as shown in Figure [[#img-16|16]] <span id='citeF-27'></span>[[#cite-27|[27]]]. The mechanical problem in the crack-free region is solved using the standard FEM.
857
858
<div id='img-16'></div>
859
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
860
|-
861
|[[Image:Draft_Samper_773076048-equivalence-fem-dem.png|600px|Equivalence between stiffness matrix (FEM) and cohesive link (DEM)]]
862
|- style="text-align: center; font-size: 75%;"
863
| colspan="1" | '''Figure 16:''' Equivalence between stiffness matrix (FEM) and cohesive link (DEM)
864
|}
865
866
A damage model has been implemented in the FEM formulation to simulate the response of the continuum domain. The onset of a crack depends on the damage level in each triangular element  calculated as an internal variable (d) at each point according to a failure criteria. The stress over the edge is computed as a mean value of the stresses in the elements sharing that edge. The stress state leads to a damage level calculated using the Mohr-Coulomb failure criteria.
867
868
A smoothing technique is employed to determine the damage over the element edges. Once the damage limit is reached, a stiffness loss is assumed in the triangle element associated to the area determined by the centroid of the triangle and the damaged edge as shown in Figure [[#img-17|17]].
869
870
<div id='img-17'></div>
871
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
872
|-
873
|[[Image:Draft_Samper_773076048-equivalence-fem-dem-damaged.png|600px|Equivalence between stiffness matrix (FEM) and cohesive link (DEM) with a damaged edge]]
874
|- style="text-align: center; font-size: 75%;"
875
| colspan="1" | '''Figure 17:''' Equivalence between stiffness matrix (FEM) and cohesive link (DEM) with a damaged edge
876
|}
877
878
<div id='img-18'></div>
879
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
880
|-
881
|[[Image:Draft_Samper_773076048-3-noded-triangle.png|300px|3-noded triangle with two sides damaged]]
882
|- style="text-align: center; font-size: 75%;"
883
| colspan="1" | '''Figure 18:''' 3-noded triangle with two sides damaged
884
|}
885
886
The largest damage parameter in two of the three sides of the element (Figure [[#img-18|18]]) determines the stiffness matrix of the element that is recalculated at every time step as follows,
887
888
<span id="eq-21"></span>
889
{| class="formulaSCP" style="width: 100%; text-align: left;" 
890
|-
891
| 
892
{| style="text-align: left; margin:auto;width: 100%;" 
893
|-
894
| style="text-align: center;" | <math>{\boldsymbol K}^{(e)}=\left[1-\frac{d_i+d_j}{2}\right]{\boldsymbol K}^{(e)}_0 </math>
895
|}
896
| style="width: 5px;text-align: right;white-space: nowrap;" | (21)
897
|}
898
899
where <math display="inline">{\boldsymbol K}_0</math> is the initial stiffness matrix of the undamaged element and <math display="inline">d_i</math> and <math display="inline">d_j</math> are the two maximum values of the damage parameter for the three element sides.
900
901
When a cohesive bond is removed, the side nodes of the element are disconnected and two discrete particles are introduced precisely at the same nodal position. Their radius and mass are defined as the maximum ones to guarantee the contact without overlapping other discrete elements so as to avoid spurious forces from the numerical approximation. There are several algorithms to define these discrete elements in the literature <span id='citeF-12'></span>[[#cite-12|[12]]].
902
903
Once the discrete particles are created their interaction is governed according the contact and friction laws, as in the general DEM formulation described earlier.
904
905
A relevant point in this FEM-DEM approach is its low computational cost. Most of the cost in a DEM simulation is due to the contact searching algorithms. In the FEM-DEM technique the number of discrete elements is in general much lower than the number of nodes because the fractured area is expected to be smaller than the whole calculation domain. Hence, the computational time consumed by the searching algorithms is much lower than in a standard DEM solution.
906
907
===5.2 Numerical model===
908
909
The pulse fracture problem described in Section [[#4.6 DEM simulation of pulse fracture of a shale rock mass|4.6]] has been solved with the FEM-DEM technique presented above assuming a double symmetry respect to the wellbore. A 2D square domain of 8 m side length has been discretized using the finite element mesh of 22128 3-noded linear triangles and 11237 nodes as depicted in Figure [[#img-19|19]]. The size of the elements is variable from a minimum of 12 mm in the finer bound next to the wellbore to a maximum of 900 mm far from the wellbore where fracture is not expected.
910
911
<div id='img-19'></div>
912
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
913
|-
914
|[[Image:Draft_Samper_773076048-model-general-view.png|600px|Numerical model. General view and wellbore detail]]
915
|- style="text-align: center; font-size: 75%;"
916
| colspan="1" | '''Figure 19:''' Numerical model. General view and wellbore detail
917
|}
918
919
The pressure pulse has been applied in the round bound of the wellbore by means of a distributed load with the time evolution function described in Section 3.3.2. The effect of the gas pressure within the cracks has been modeled by applying the pressure force (deduced from the values of Figure [[#img-9|9]]) at the side of the elements adjacent to the crack. The problem is assumed to be dynamic.  An explicit scheme has been employed for the time integration of the discretized dynamic equations in both the FEM and DEM domains.
920
921
===5.3 Material parameters===
922
923
The material properties of the shale rock correspond to the constitutive parameters presented in <span id='citeF-21'></span>[[#cite-21|[21]]]. Table [[#table-5|5]] list the values used in the numerical simulation.
924
925
926
{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
927
|+ style="font-size: 75%;" |<span id='table-5'></span>Table. 5 Constitutive parameters for shale rock. FEM-DEM analysis
928
929
|- style="border-top: 2px solid;"
930
| colspan='4' | '''Shale rock constitutive parameters. FEM-DEM model'''
931
932
933
|- style="border-top: 2px solid;"
934
| style="text-align: left;" |  Parameter 
935
| Notation 
936
| Value 
937
| Units 
938
939
|- style="border-top: 2px solid;"
940
| style="text-align: left;" |   Density 
941
| <math>\rho </math>
942
| 2550 
943
| kg/m<math display="inline">^{3}</math> 
944
|-
945
| style="text-align: left;" | Young modulus 
946
| E 
947
| 30 000 
948
| MPa 
949
|-
950
| style="text-align: left;" | Poison coefficient 
951
| <math>\nu </math>
952
| 0.20 
953
| [ - ] 
954
|-
955
| style="text-align: left;" | Yield strength 
956
| <math>\sigma _{0}</math>
957
| 5.0 
958
| MPa 
959
|-
960
| style="text-align: left;" | Fracture energy 
961
| G 
962
| 32 
963
| J/m<math display="inline">^{2}</math> 
964
965
966
|}
967
968
===5.4 FEM-DEM results of pulse fracture simulation of shale rock mass===
969
970
Three different cases of pulse fracturing of shale rock mass have been simulated using the FEM-DEM approach for different depths as described in Table 3 for the DEM simulation, with the corresponding initial stress states. For every case the pattern and the fracture length have been measured and compared to the reference values in <span id='citeF-21'></span>[[#cite-21|[21]]]. Results have been treated to show the whole fracture by a symmetric extension of the results obtained in the quarter domain.
971
972
Figures [[#img-20|20]] to  [[#img-22|22]] show the crack pattern obtained with the FEM-DEM formulation compared to the reference results <span id='citeF-21'></span>[[#cite-21|[21]]] depicted in a box for the three cases considered.
973
974
<div id='img-20'></div>
975
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
976
|-
977
|[[Image:Draft_Samper_773076048-crack-simulation-depth500.png|600px|Crack simulation with FEM-DEM model. Depth 500 ft. FEM results in box from <span id='citeF-21'></span>[[#cite-21|[21]]]]]
978
|- style="text-align: center; font-size: 75%;"
979
| colspan="1" | '''Figure 20:''' Crack simulation with FEM-DEM model. Depth 500 ft. FEM results in box from <span id='citeF-21'></span>[[#cite-21|[21]]]
980
|}
981
982
<div id='img-21'></div>
983
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
984
|-
985
|[[Image:Draft_Samper_773076048-crack-simulation-depth5000.png|600px|Crack simulation with FEM-DEM model. Depth 5000 ft. FEM results in box from <span id='citeF-21'></span>[[#cite-21|[21]]]]]
986
|- style="text-align: center; font-size: 75%;"
987
| colspan="1" | '''Figure 21:''' Crack simulation with FEM-DEM model. Depth 5000 ft. FEM results in box from <span id='citeF-21'></span>[[#cite-21|[21]]]
988
|}
989
990
<div id='img-22'></div>
991
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
992
|-
993
|[[Image:Draft_Samper_773076048-crack-simulation-depth10000.png|600px|Crack simulation with FEM-DEM model. Depth 10000 ft. FEM results in box from <span id='citeF-21'></span>[[#cite-21|[21]]]]]
994
|- style="text-align: center; font-size: 75%;"
995
| colspan="1" | '''Figure 22:''' Crack simulation with FEM-DEM model. Depth 10000 ft. FEM results in box from <span id='citeF-21'></span>[[#cite-21|[21]]]
996
|}
997
998
The lengths of the cracks obtained with the FEM-DEM approach are shown in Table [[#table-6|6]]. Good agreement with the results reported in  <span id='citeF-21'></span>[[#cite-21|[21]]] is obtained.
999
1000
1001
{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
1002
|+ style="font-size: 75%;" |<span id='table-6'></span>Table. 6 Fracture length. DEM-FEM results vs reference values
1003
1004
|- style="border-top: 2px solid;"
1005
| style="text-align: left;" |  Case 
1006
| FEM (ft) 
1007
| DEM-FEM(ft) 
1008
1009
1010
|- style="border-top: 2px solid;"
1011
| style="text-align: left;" |  Depth 500 ft 
1012
| 40 
1013
| 33 
1014
|-
1015
| style="text-align: left;" | Depth 5000 ft 
1016
| 5 
1017
| 5.4 
1018
|-
1019
| style="text-align: left;" | Depth 10000 
1020
| 2.3 
1021
| 3.2 
1022
1023
1024
|}
1025
1026
==6 Concluding remarks==
1027
1028
Two different numerical methods based on DEM and the FEM-DEM techniques have been presented and applied to the simulation of cracking  in a shale rock reservoir under a pressure pulse.
1029
1030
The DEM is able to reproduce the multicracking pattern in the shale rock induced by the pressure pulse. The cracking pattern and its length is consistent with the reference FEM results.
1031
1032
A FEM-DEM technique has also been applied to the same problem. The cracks obtained with this model reproduce the crack pattern and length obtained with the FEM (and also the standard DEM) with a lower computational cost.
1033
1034
The results presented show the capabilities of the DEM and the FEM-DEM approaches here described to simulate the cracking behavior of shale rock and other geomaterials under a pressure pulse.
1035
1036
The effect of the gas pressure within the cracks can be accounted for by coupling the FEM-DEM approach with an embedded FEM solver for a compressible gas. The coupling strategy solves the equation for the gas flow in the finite element mesh that fills the spaces created by cracks. Details are given in <span id='citeF-28'></span>[[#cite-28|[28]]]
1037
1038
The extension of both numerical models to 3D problems is straightforward and will be reported in future publications.
1039
1040
==Acknowledgements==
1041
1042
This work has been carried out with the financial support from Advanced grant projects COMDESMAT and ICEBREAKER of the European Research Council and the BALAMED project (BIA2012-39172) of MINECO (Spain). The support of CIMNE for making available the DEMPack code (www.cimne.com/dempack) and the GiD pre-postprocessor (www.gidhome.com) is gratefully acknowledged.
1043
1044
===BIBLIOGRAPHY===
1045
1046
<div id="cite-1"></div>
1047
'''[[#citeF-1|[1]]]'''  Cervera M., Chiumenti M., Codina R. (2010) Mixed stabilized finite element methods in nonlinear solid mechanics Part I: Formulation. Computer Methods in Applied Mechanics and Engineering 199:2559&#8211;2570
1048
1049
<div id="cite-2"></div>
1050
'''[[#citeF-2|[2]]]'''  Cervera M., Chiumenti M., Codina R. (2010) Mixed stabilized finite element methods in nonlinear solid mechanics Part II: Strain localization. Computer Methods in Applied Mechanics and Engineering 199:2571&#8211;2589
1051
1052
<div id="cite-3"></div>
1053
'''[[#citeF-3|[3]]]'''  Cervera M., Chiumenti M., Codina R. (2011) Mesh objective modelling of cracks using continuous linear strain and displacements interpolations. Int. J Num. Meth. Eng.  87:962&#8211;987
1054
1055
<div id="cite-4"></div>
1056
'''[[#citeF-4|[4]]]'''  Cundall PA, Strack ODL. (1979) A discrete numerical method for granular assemblies. Geotechnique 2:47&#8211;65.
1057
1058
<div id="cite-5"></div>
1059
'''[[#citeF-5|[5]]]'''  Donzé F., Richelieu F., Magnier S. (2009) Advances in discrete element method applied to soil, rock and concrete mechanics in state of the art of geotechnical engineering. Electro J Geotechnical Eng. 8:1-44.
1060
1061
<div id="cite-6"></div>
1062
'''[[#citeF-6|[6]]]'''  GiD The personal pre and postprocessor. Version 13.0 (2017) (www.gidhome.com).
1063
1064
<div id="cite-7"></div>
1065
'''[[#citeF-7|[7]]]'''  Hernandez J.A., Oliver J., Cante J.C., Weyler R. (2011)  Numerical modeling of crack formation in powder forming processes.  International Journal of Solids and Structures. 48(2):292–316.
1066
1067
<div id="cite-8"></div>
1068
'''[[#citeF-8|[8]]]'''  Hsiegh Y-M, Li H-H, Huang T-H, Jeng F-S. (2008) Interpretations on how the macroscopic mechanical behavior of sandstone affected by macroscopic properties revealed by bonded-particle model. Eng. Geol. 99(1):1-10.
1069
1070
<div id="cite-9"></div>
1071
'''[[#citeF-9|[9]]]'''  Huang H. (1999) ''Discrete element modelling of tool-rock interaction''. Ph.D Thesis, University of Minnesota.
1072
1073
<div id="cite-10"></div>
1074
'''[[#citeF-10|[10]]]'''  Johnson P.R., Petrinic N.,  Süli E. (2005) Element-splitting for simulation of fracture in 3D Solid continua. VIII International Conference on Computational Plasticity, Barcelona.
1075
1076
<div id="cite-11"></div>
1077
'''[[#citeF-11|[11]]]'''  Katagiri S.,  Takada S. (2003) Development of FEM-DEM combined method for fracture analysis of a continuous media. Memoirs of the Graduate School of Science and Technology, Kobe University. 20A, 65-79. Kobe University
1078
1079
<div id="cite-12"></div>
1080
'''[[#citeF-12|[12]]]'''  Labra C., Oñate E. (2009) High-density sphere packing for discrete element method simulations Commun. Numer. Meth. Engng. 25(7):837&#8211;849.
1081
1082
<div id="cite-13"></div>
1083
'''[[#citeF-13|[13]]]'''  Labra C., Rojek J., Oñate E., Zarate F. (2008) Advances in discrete element modelling of underground excavations. Acta Geotech 3(4):317&#8211;322.
1084
1085
<div id="cite-14"></div>
1086
'''[[#citeF-14|[14]]]'''  Munjiza A. (2004) The combined finite-discrete element method. Wiley & Son.
1087
1088
<div id="cite-15"></div>
1089
'''[[#citeF-15|[15]]]'''  Oliver J., Caicedo M., Roubin E., Huespe A.E., Hernandez J.A. (2015) Continuum approach to computational multiscale modeling of propagating fracture. Comput. Meth. in Appl.Mech. and Engng, 294:384&#8211;427.
1090
1091
<div id="cite-16"></div>
1092
'''[[#citeF-16|[16]]]'''  Oliver J., Huespe A.E., Dias I.F. (2012) Strain localization, strong discontinuities and material fracture: Matches and mismatches Computer Methods in Applied Mechanics and Engineering. Volume 241-244, 1 October 2012, Pages 323-336
1093
1094
<div id="cite-17"></div>
1095
'''[[#citeF-17|[17]]]'''  Oñate E., Zarate F., Miquel J., Santasusana M., Celigueta M.A., Arrufat F., Gankikota R., Valiullin K.,  Ring L. (2015) A local constitutive model for the discrete element method. Application to geomaterials and concrete. Computational Particle Mechanics. 2:139&#8211;160.
1096
1097
<div id="cite-18"></div>
1098
'''[[#citeF-18|[18]]]'''  Oñate T., Rojek J. (2004) Combination of dicrete element and finite element methods for dynamic analysis of geomechanics problems. Comput. Methods Appl. Mech. Eng. 193:3087-3128.
1099
1100
<div id="cite-19"></div>
1101
'''[[#citeF-19|[19]]]'''  Zarate F.,  Oñate E. (2015) A simple FEM-DEM technique for fracture prediction in materials and structures Computational Particle Mechanics, 2(3):301-314.
1102
1103
<div id="cite-20"></div>
1104
'''[[#citeF-20|[20]]]'''  Potyondy D., Cundall P. (2004) A bonded-particle model for rock. Int J Rock Mech Min Sci 41(8):1329-1364. Rock Mechanics Results from the Underground Research Laboratory, Canada.
1105
1106
<div id="cite-21"></div>
1107
'''[[#citeF-21|[21]]]'''  Reza Safari M., Gandikota R., Mutlu U., Ji M., Glaville J.,  Abass H. (2013) Pulsed fracturing in shale reservoirs: geomechanical aspects, ductile-brittle transition and field implications. Unconventional Resources Tecnology Conference (URTeC), Denver, CO, USA, 12-14  Aug. 2013, pp. 448&#8211;461.
1108
1109
<div id="cite-22"></div>
1110
'''[[#citeF-22|[22]]]'''  Reza Safari M., Huang J.,  Mutlu U. (2013) Ductile to brittle transition, generation of complex fracture networks and engineering implications. Applied Geoscience Conference. Houston (Texas).
1111
1112
<div id="cite-23"></div>
1113
'''[[#citeF-23|[23]]]'''  Rojek, J., Labra C., Su O.,  Oñate E. (2012) Comparative study of different micromechanical parameters. Int J Solids Struc 49(13):1497-1517.
1114
1115
<div id="cite-24"></div>
1116
'''[[#citeF-24|[24]]]'''  Rojek J., Oñate E. (2007) Multiscale analysis using a coupled discrete/finite element model. Interact. Multiscale Mech 1(1):1-31.
1117
1118
<div id="cite-25"></div>
1119
'''[[#citeF-25|[25]]]'''  Santasusana M., Irazábal J., Oñate E., Carbonell J.M. (2016) The Double Hierarchy Method. A parallel 3D contact method for the interaction of spherical particles with rigid FE boundaries using the DEM. Comp. Part. Mech. 3:407&#8211;428.
1120
1121
<div id="cite-26"></div>
1122
'''[[#citeF-26|[26]]]'''  Tran V-T., Donzé F-V., Marin P. (2011) A discrete element model of concrete under high triaxial loading. Cem Concr Compos 33(9):936-948.
1123
1124
<div id="cite-27"></div>
1125
'''[[#citeF-27|[27]]]'''  Ubach P.A., Arrufat F., Ring L., Gandikota R., Zarate F.,  Oñate E. (2016)  Application of an enhanced discrete element method to oil and gas drilling processes. Comp. Part. Mech. 3(1):29&#8211;41.
1126
1127
<div id="cite-28"></div>
1128
'''[[#citeF-28|[28]]]'''  Zárate F., Löhner R.,  Oñate E. (2017) Modeling of multifracture in solids induced by bleat load accounting for coupled gas-structure interaction effects. Comput. Particle Mechanis, (In preparation).
1129

Return to Gonzalez et al 2018b.

Back to Top