Pulse fracture simulation in shale rock reservoirs. DEM and FEM-DEM approaches

<!-- metadata commented in wiki content
====

''', , José Manuel González, <math>\cdot </math>Francisco Zárate, <math>\cdot </math>Eugenio Oñate'''
-->

==Abstract==

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.

'''Keywords''': discrete element method, finite element method, FEM-DEM, pulse fracture, shale rock.

==1 Introduction==

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.

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]]].

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]]].

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.

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.

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.

==2 Pulse fracturing problem==

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.

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.

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.

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.

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.

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.

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.

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]]].

==3 DEM numerical simulation==

===3.1 DEM model overview===

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.

====3.1.1 Kinematic equations and integration scheme====

The translation and rotation of the discrete particles is governed by the standard dynamics equations for rigid bodies,

<span id="eq-1"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| 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>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (1)
|}

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]]). 

<div id='img-1'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_773076048-Fig2con327.png|350px|Motion of a rigid particle]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 1:''' Motion of a rigid particle
|}

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]])

<span id="eq-2"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| 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>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (2)
|}

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.

<div id='img-2a'></div>
<div id='img-2b'></div>
<div id='img-2'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_773076048-contact-interface.png|340px|]]
|[[Image:Draft_Samper_773076048-Fig6con327.png|200px|]]
|- style="text-align: center; font-size: 75%;"
| (a) 
| (b) 
|- style="text-align: center; font-size: 75%;"
| 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
|}

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

<span id="eq-3"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| 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>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
|}

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.

<span id="eq-4"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\Delta t\le {\Delta t}_{cr}=\frac{2}{{\omega }^{max}}\left(\sqrt{1+{\xi }^2}-\xi \right) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
|}

where <math display="inline">\xi </math> is a fraction of the critical damping [ ].

====3.1.2 Forces acting over the discrete element====

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]]).

The normal component of the interaction forces is calculated as,

<span id="eq-5"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| 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>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (5)
|}

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,

<span id="eq-6"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\sigma }_n=E{\varepsilon }_n+c{\dot{\varepsilon }}_n </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (6)
|}

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

<span id="eq-7"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| 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>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (7)
|}

where <math display="inline">K_n</math> is the normal stiffness parameter (see Eq.[[#eq-9|9]]).

The normal strain and strain rate values are computed from the kinematic variables as,

<span id="eq-8"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| 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>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
|}

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).

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

<span id="eq-9"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| 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>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (9)
|}

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]].

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]]]

<span id="eq-10"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| 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>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (10)
|}

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,

<span id="eq-11"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol{u}}^{ij}_s={\boldsymbol{u}}^{ij}-\left({\boldsymbol{u}}^{ij} , {\boldsymbol{n}}^{ij}\right){\boldsymbol{n}}^{ij} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (11)
|}

Equation [[#eq-10|10]] leads to a general expression for the shear stiffness (assumed to be the same for both shear directions), as

<span id="eq-12"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>K_s= G \frac{{\alpha }_{ij}{ A}_{ij}}{d_{ij}}=\frac{K_n}{2\left(1+\nu \right)} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (12)
|}

where <math display="inline">\nu </math> is the Poisson's ratio of the material.

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]]]:

<span id="eq-13.a"></span>
<span id="eq-13.b"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| 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>
| style="width: 5px;text-align: right;white-space: nowrap;" | (13.a)
|-
| 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>
| style="width: 5px;text-align: right;white-space: nowrap;" | (13.b)
|}
|}

More details of the local constitutive model for the DEM can be found in <span id='citeF-17'></span>[[#cite-17|[17]]].

===3.2 Normal and shear failure===

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

<span id="eq-14"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>F_{n_t}  \ge {\cal F}_{n_t}\quad , \quad F_s  \ge {\cal F}_s </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (14)
|}

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]]).

The interface strengths are defined as

<span id="eq-15"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| 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>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (15)
|}

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>).

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

<span id="eq-16"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\tau ^f = \beta (\sigma _{n_c}^f)_{UCS}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (16)
|}

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>.

Following tension failure, the constitutive behavior in the shear direction is governed by the standard Coulomb law

<span id="eq-17"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| 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>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (17)
|}

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.

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.

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]]].

<div id='img-3'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_773076048-Fig8a.png|250px|Failure line in terms of normal and shear forces. Uncoupled failure model]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 3:''' Failure line in terms of normal and shear forces. Uncoupled failure model
|}

<div id='img-4'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_773076048-Fig8b.png|500px|Undamaged and damaged elastic moduli under  tension (a) and  shear (b) forces]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 4:''' Undamaged and damaged elastic moduli under  tension (a) and  shear (b) forces
|}

===3.3 Elasto-plastic model for compression forces===

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.

The ''elasto-plastic relationships'' in the normal compressive direction are defined as<br/>

'''Loading path'''
<span id="eq-18"></span>

<span id="eq-18.a"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>dF_{n_c}=K_{T_n} du_n  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (18.a)
|}

'''Unloading path'''

<span id="eq-18.b"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>dF_{n_c}=K_{n_0}  du_n  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (18.b)
|}

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

<span id="eq-19"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>K_{T_n}=\frac{E_T}{E_0}K_{n_0}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (19)
|}

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]]).

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]]).

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.

<div id='img-5'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[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]]
|- style="text-align: center; font-size: 75%;"
| 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
|}

===3.4 Calibration of the material parameters===

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]]].


{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
|+ style="font-size: 75%;" |<span id='table-1'></span>Table. 1 DEM constitutive parameters analysis of  shale rock

|- style="border-top: 2px solid;"
| <math display="inline">\rho </math> (g/cc)
| <math>\mu _1</math>
| <math>\mu _2</math>
| <math>E_0</math>(GPa)
| <math>\nu </math>
| <math>\sigma ^f_{t}</math>(MPa) 
| <math>\tau _f</math>(MPa) 

|- style="border-top: 2px solid;"
|  2.55 
| 0.7 
| 0.6 
| 30 
| 0.2 
| 5.0 
| 25

|- style="border-top: 2px solid;"
|  LCS1 
| LCS2 
| LCS3
| YRC1 
| YRC2 
| YRC3 
| <math>\alpha </math>
|-
| (MPa) 
| (MPa) 
| (MPa) 
| 
| 
| 
| 
|- style="border-top: 2px solid;border-bottom: 2px solid;"
| 20 
| 30 
| 40 
| 2
| 5 
| 10
| 1.0

|}

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.

===3.5 UCS and BTS tests on shale rock material===

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]]].

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]]].

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]]].

<div id='img-6'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[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]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 6:''' Axial stress-axial strain curve for  UCS test in shale rock material. DEM results using 37000 spheres
|}

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.

<div id='img-7'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_773076048-Fig28.png|450px|Tensile stress-time curve for BTS test in shale rock material. DEM results using 27000 spheres]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 7:''' Tensile stress-time curve for BTS test in shale rock material. DEM results using 27000 spheres
|}

===3.6 DEM simulation of pulse fracture of a shale rock mass===

We present the application of the DEM technique previously described to the fracture of a shale mass under a pulse load.

====3.6.1 Geometry and mesh====

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.

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.

The analysis domain has been discretized with the DEM using 882693 circular discrete elements of 25 mm radius (average).

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.

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.

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:

*  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.

*  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.

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.

<div id='img-8'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_773076048-geometry.png|400px|Numerical model. Geometry and applied loads]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 8:''' Numerical model. Geometry and applied loads
|}

====3.6.2 Pressure pulse application====

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

<span id="eq-20"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| 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>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (20)
|}

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>.

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.

<div id='img-9'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_773076048-loading-pulse.png|450px|Loading pulse pressure function]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 9:''' Loading pulse pressure function
|}

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]].


{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
|+ style="font-size: 75%;" |<span id='table-2'></span>Table. 2 Characterization of the pulse loading curves

|- style="border-top: 2px solid;"
| style="text-align: left;" |  '''Pulse''' 
| '''Peak'''  (MPa) 
| '''Loading rate''' (MPa/ms) 
| '''Peak time''' (ms) 


|- style="border-top: 2px solid;"
| style="text-align: left;" |  Reference 
| 100 
| 200 
| 0.50 
|-
| style="text-align: left;" | Peak 250 
| 250 
| 200 
| 1.25 


|}

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.

====3.6.3 Initial stress state====

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.

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.

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]]].

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.


{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
|+ style="font-size: 75%;" |<span id='table-3'></span>Table. 3 Initial stresses in the shale rock mass for different depths

|- style="border-top: 2px solid;"
| style="text-align: left;" |  '''Depth''' (ft) 
| <math>\sigma _{x}</math>(MPa) 
| <math>\sigma _{y}</math>(MPa) 


|- style="border-top: 2px solid;"
| style="text-align: left;" |  0 
| - 
| - 
|-
| style="text-align: left;" | 500 
| 1.0 
| 1.2 
|-
| style="text-align: left;" | 5 000 
| 10.0 
| 12.0 
|-
| style="text-align: left;" | 10 000 
| 20.0 
| 24.0 


|}

===3.7 DEM results===

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.

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.

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.


{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
|+ 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]]]

|- style="border-top: 2px solid;"
|  Depth  
| <math>\sigma _{x}</math>
| <math>\sigma _{y}</math>
| Pressure peak  
| Fracture length (ft) 
| Fracture length (ft)
| Figure
|-
| (ft) 
| (MPa)
| (MPa)
| (MPa)
| FEM result <span id='citeF-21'></span>[[#cite-21|[21]]]
| DEM result 
| 


|- style="border-top: 2px solid;"
|   0 
| 0 
| 0 
| 100 
| &#8211; 
| 40
| 10
|-
| 500  
| 1.0 
| 1.2 
| 100 
| 40 
| 33
| 11  
|-
| 5000  
| 10 
| 12 
| 100 
| 5 
| 5.4 
| 12
|-
| 10000 
| 20 
| 24 
| 100 
| 2.3 
| 3.2 
| 13
|-
| 5000 (weak anisotropy) 
| 5 
| 12 
| 100 
| &#8211; 
| 10.5 
| 14a
|-
| 5000 (strong anisotropy) 
| 2.5 
| 12 
| 100 
| &#8211; 
| 15,2 
| 14b
|-
| 5000  
| 10 
| 12 
| 250 
| 28 
| 30.2 
| 15


|}

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.

====3.7.1 Fracture vs initial stress state====

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.

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]]).

<div id='img-10a'></div>
<div id='img-10b'></div>
<div id='img-10'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_773076048-damage-no-initial-1.png|400px|]]
|- style="text-align: center; font-size: 75%;"
| (a)
|-
|[[Image:Draft_Samper_773076048-damage-no-initial-2.png|400px|]]
|- style="text-align: center; font-size: 75%;"
| (b)
|-
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 10:''' Damage distribution (a) and type of fracture (b). No initial stress. Peak pressure load: 100 MPa
|}

<div id='img-11a'></div>
<div id='img-11b'></div>
<div id='img-11c'></div>
<div id='img-11'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
| colspan="2"|[[Image:Draft_Samper_773076048-Fig11b_1.png|450px|]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | (a)
|-
|[[Image:Draft_Samper_773076048-Fig11b_2.png|450px|]]
|[[Image:Draft_Samper_773076048-Fig11b_3.png|450px|]] 
|- style="text-align: center; font-size: 75%;"
| (b) 
| (c)
|- style="text-align: center; font-size: 75%;"
| 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]]]
|}

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.

<div id='img-12a'></div>
<div id='img-12b'></div>
<div id='img-12c'></div>
<div id='img-12'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
| colspan="2"| [[Image:Draft_Samper_773076048-damage-depth5000-1.png|450px|]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | (a)
|-
| [[Image:Draft_Samper_773076048-Fig12b_1.png|450px|]]
| [[Image:Draft_Samper_773076048-Fig12b_2.png|450px|]]
|- style="text-align: center; font-size: 75%;"
| (b) 
| (c)
|- style="text-align: center; font-size: 75%;"
| 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]]]
|}


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.

<div id='img-13a'></div>
<div id='img-13b'></div>
<div id='img-13c'></div>
<div id='img-13'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
| colspan="2"|[[Image:Draft_Samper_773076048-damage-depth10000-1.png|450px|]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | (a)
|-
|[[Image:Draft_Samper_773076048-Fig13b_1.png|450px|]]
|[[Image:Draft_Samper_773076048-Fig13b_2.png|450px|]]
|- style="text-align: center; font-size: 75%;"
| (b) 
| (c)
|- style="text-align: center; font-size: 75%;"
| 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]]]
|}

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.

====3.7.2 Fracture vs anisotropy====

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.

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.

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]].

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.

<div id='img-14a'></div>
<div id='img-14b'></div>
<div id='img-14'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_773076048-Fig14a.png|500px|Weak anisotropy σₓ = 5 MPa, σ<sub>y</sub> =12 MPa]]
|[[Image:Draft_Samper_773076048-Fig14b.png|500px|Strong anisotropy σₓ = 2.5 MPa, σ<sub>y</sub> =12  MPa]]
|- style="text-align: center; font-size: 75%;"
| (a) Weak anisotropy <math display="inline">\sigma _{x} = 5</math> MPa, <math display="inline">\sigma _{y} =12</math> MPa
| (b) Strong anisotropy <math display="inline">\sigma _{x} = 2.5</math> MPa, <math display="inline">\sigma _{y} =12 </math> MPa
|- style="text-align: center; font-size: 75%;"
| 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
|}

===3.8 Fracture vs peak load===

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]]].

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]]].

<div id='img-15a'></div>
<div id='img-15b'></div>
<div id='img-15c'></div>
<div id='img-15'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_773076048-Fig15a.png|450px|]]
|[[Image:Draft_Samper_773076048-Fig15b_1.png|450px|]]
|- style="text-align: center; font-size: 75%;"
| (a) 
| (b) 
|-
| colspan="2"|[[Image:Draft_Samper_773076048-Fig15b_2.png|450px|]]
|- style="text-align: center; font-size: 75%;"
|  colspan="2" | (c) 
|- style="text-align: center; font-size: 75%;"
| 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]]]
|}

==5 FEM-DEM numerical simulation==

===5.1 Constitutive model overview===

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.

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.

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.

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.

<div id='img-16'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_773076048-equivalence-fem-dem.png|600px|Equivalence between stiffness matrix (FEM) and cohesive link (DEM)]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 16:''' Equivalence between stiffness matrix (FEM) and cohesive link (DEM)
|}

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.

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]].

<div id='img-17'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_773076048-equivalence-fem-dem-damaged.png|600px|Equivalence between stiffness matrix (FEM) and cohesive link (DEM) with a damaged edge]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 17:''' Equivalence between stiffness matrix (FEM) and cohesive link (DEM) with a damaged edge
|}

<div id='img-18'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_773076048-3-noded-triangle.png|300px|3-noded triangle with two sides damaged]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 18:''' 3-noded triangle with two sides damaged
|}

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,

<span id="eq-21"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol K}^{(e)}=\left[1-\frac{d_i+d_j}{2}\right]{\boldsymbol K}^{(e)}_0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (21)
|}

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.

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]]].

Once the discrete particles are created their interaction is governed according the contact and friction laws, as in the general DEM formulation described earlier.

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.

===5.2 Numerical model===

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.

<div id='img-19'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_773076048-model-general-view.png|600px|Numerical model. General view and wellbore detail]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 19:''' Numerical model. General view and wellbore detail
|}

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.

===5.3 Material parameters===

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.


{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
|+ style="font-size: 75%;" |<span id='table-5'></span>Table. 5 Constitutive parameters for shale rock. FEM-DEM analysis

|- style="border-top: 2px solid;"
| colspan='4' | '''Shale rock constitutive parameters. FEM-DEM model'''


|- style="border-top: 2px solid;"
| style="text-align: left;" |  Parameter 
| Notation 
| Value 
| Units 

|- style="border-top: 2px solid;"
| style="text-align: left;" |   Density 
| <math>\rho </math>
| 2550 
| kg/m<math display="inline">^{3}</math> 
|-
| style="text-align: left;" | Young modulus 
| E 
| 30 000 
| MPa 
|-
| style="text-align: left;" | Poison coefficient 
| <math>\nu </math>
| 0.20 
| [ - ] 
|-
| style="text-align: left;" | Yield strength 
| <math>\sigma _{0}</math>
| 5.0 
| MPa 
|-
| style="text-align: left;" | Fracture energy 
| G 
| 32 
| J/m<math display="inline">^{2}</math> 


|}

===5.4 FEM-DEM results of pulse fracture simulation of shale rock mass===

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.

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.

<div id='img-20'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[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]]]]]
|- style="text-align: center; font-size: 75%;"
| 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]]]
|}

<div id='img-21'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[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]]]]]
|- style="text-align: center; font-size: 75%;"
| 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]]]
|}

<div id='img-22'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[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]]]]]
|- style="text-align: center; font-size: 75%;"
| 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]]]
|}

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.


{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
|+ style="font-size: 75%;" |<span id='table-6'></span>Table. 6 Fracture length. DEM-FEM results vs reference values

|- style="border-top: 2px solid;"
| style="text-align: left;" |  Case 
| FEM (ft) 
| DEM-FEM(ft) 


|- style="border-top: 2px solid;"
| style="text-align: left;" |  Depth 500 ft 
| 40 
| 33 
|-
| style="text-align: left;" | Depth 5000 ft 
| 5 
| 5.4 
|-
| style="text-align: left;" | Depth 10000 
| 2.3 
| 3.2 


|}

==6 Concluding remarks==

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.

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.

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.

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.

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]]]

The extension of both numerical models to 3D problems is straightforward and will be reported in future publications.

==Acknowledgements==

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.

===BIBLIOGRAPHY===

<div id="cite-1"></div>
'''[[#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

<div id="cite-2"></div>
'''[[#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

<div id="cite-3"></div>
'''[[#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

<div id="cite-4"></div>
'''[[#citeF-4|[4]]]'''  Cundall PA, Strack ODL. (1979) A discrete numerical method for granular assemblies. Geotechnique 2:47&#8211;65.

<div id="cite-5"></div>
'''[[#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.

<div id="cite-6"></div>
'''[[#citeF-6|[6]]]'''  GiD The personal pre and postprocessor. Version 13.0 (2017) (www.gidhome.com).

<div id="cite-7"></div>
'''[[#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.

<div id="cite-8"></div>
'''[[#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.

<div id="cite-9"></div>
'''[[#citeF-9|[9]]]'''  Huang H. (1999) ''Discrete element modelling of tool-rock interaction''. Ph.D Thesis, University of Minnesota.

<div id="cite-10"></div>
'''[[#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.

<div id="cite-11"></div>
'''[[#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

<div id="cite-12"></div>
'''[[#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.

<div id="cite-13"></div>
'''[[#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.

<div id="cite-14"></div>
'''[[#citeF-14|[14]]]'''  Munjiza A. (2004) The combined finite-discrete element method. Wiley & Son.

<div id="cite-15"></div>
'''[[#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.

<div id="cite-16"></div>
'''[[#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

<div id="cite-17"></div>
'''[[#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.

<div id="cite-18"></div>
'''[[#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.

<div id="cite-19"></div>
'''[[#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.

<div id="cite-20"></div>
'''[[#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.

<div id="cite-21"></div>
'''[[#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.

<div id="cite-22"></div>
'''[[#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).

<div id="cite-23"></div>
'''[[#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.

<div id="cite-24"></div>
'''[[#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.

<div id="cite-25"></div>
'''[[#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.

<div id="cite-26"></div>
'''[[#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.

<div id="cite-27"></div>
'''[[#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.

<div id="cite-28"></div>
'''[[#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).