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
==COMBINATION OF DISCRETE ELEMENT AND FINITE ELEMENT METHODS FOR DYNAMIC ANALYSIS OF GEOMECHANICS PROBLEMS==
3
4
'''E. Oñate<math>^\dagger </math> and J. Rojek<math>^{\dagger \ddagger }</math>'''
5
6
{|
7
|-
8
|<math>^\dagger </math>International Centre for Numerical Methods in Engineering
9
|-
10
| Universidad Politécnica de Cataluña 
11
|-
12
| Gran Capitán s/n, 08034 Barcelona, Spain 
13
|-
14
| E&#8211;mail: [mailto:onate@cimne.upc.es onate@cimne.upc.es]
15
|-
16
| Web page: [http://www.cimne.upc.es www.cimne.upc.es]
17
|-
18
| <math>^{\ddagger }</math>Institute of Fundamental Technological Research
19
|-
20
| Polish Academy of Sciences 
21
|-
22
| Swietokrzyska 21, 00-049 Warsaw, Poland 
23
|-
24
| E&#8211;mail: [mailto:jrojek@ippt.gov.pl jrojek@ippt.gov.pl]
25
|-
26
| 
27
|}
28
-->
29
30
==Abstract==
31
32
The paper presents combination of Discrete Element Method (DEM) and Finite Element Method (FEM) for dynamic analysis of geomechanics problems. Combined models can employ spherical (or cylindrical in 2D) rigid elements and  finite elements in the discretization of different parts of the system. The DEM is a suitable tool to model soil/rock materials while the FEM in many cases can be a better choice to model other parts of the system considered. A typical example can be an idealization of rock cutting with a tool discretized with finite elements and rock or soil samples modelled with discrete elements. The FEM presented in the paper allows  large elasto-plastic deformations in the solid regions. Such problems require the use of stabilized FEM  to deal with the incompressibility constraint in order  to eliminate the volumetric locking defect, especially when using triangular and tetrahedra elements with equal order interpolation for the displacement and the pressure variables. In the  paper a stabilization based on the Finite Calculus (FIC) approach is used. Both theoretical algorithms of DEM and stabilized FEM are implemented in an explicit dynamic code. The paper presents some details of both formulations. A combined numerical algorithm is described finally. Selected numerical results illustrate the possibilities and performance of discrete/finite element analysis in geomechanics problems.
33
34
==1 Introduction==
35
36
The discrete element method (DEM) is nowadays acknowledged to be an effective procedure for analysis of granular and rock materials under static and dynamic loads. Typically in the DEM the soil/rock masses are modelled by a collection of spheres interacting with each other. The normal and tangential forces between spheres govern the motion of the discrete continuum. Modelling of discontinuities by the DEM is quite simple as they are  produced by loss of adhesive contact between the spheres. Examples of application of the DEM to geomechanics and similar problems can be found in <span id='citeF-1'></span><span id='citeF-2'></span>[[#cite-1|[1,2]]].
37
38
The coupling of the finite element method with the DEM is an effective approach for problems where standard continuum  finite elements are adequate to model a part of the analysis domain, whereas the DEM is more adequate to treat other areas. Examples of such problems can be found in the interaction of solids and structures with granular media. The use of simple triangular and tetrahedra elements to model the “continuum” domain is very attractive due to the simplicity of their geometry in order to treat  frictional contact conditions. Standard linear triangles and tetrahedra can suffer however from serious drawbacks in the presence of material nonlinearities inducing a near-incompressible behaviour as the numerical solution locks in these cases. The volumetric locking defect can be eliminated by using adequate mixed formulations with different interpolations for the displacement and pressure variables <span id='citeF-3'></span>[[#cite-3|[3]]]. The use of an equal order interpolation for these variables is very attractive from the computational point of view and this is  possible if an adequate stabilized formulation is chosen.
39
40
In this paper a stabilization method based on the  finite calculus (FIC) approach  <span id='citeF-4'></span><span id='citeF-5'></span>[[#cite-4|[4,5]]] is used for deriving locking-free linear triangles and tetrahedra which can be effectively used in conjunction with the DEM.
41
42
One of the motivations of the present research is the simulation of the process of rock cutting, estimating at the same time   the wear of the cutting tool. The main physical phenomenon considered is the interaction of the tool with a rock leading  to failure of the rock and tool wear. The latter process is usually slower, although in some cases visible changes of the tool shape can be appreciated after a few minutes. The tool is assumed to be rigid in our numerical model. Two kinds of tool discretization have been used. The first one employs finite elements to represent the tool. In the second type of discretization the tool is also discretised with discrete  elements. This  allows to  modify the shape of the tool easily  by removing “worn” particles.
43
44
The formulation of the discrete element method has been extended for thermal and thermo-mechanical coupled problems in order to take into account temperature  as one of the important factors influencing the tool wear process. Heat is generated during cutting  due to friction dissipation. In the model friction between the tool and rock particles as well as friction between rock particles themselves is considered. Heat conduction through the tool and rock is analysed and the temperature distribution is obtained. Temperature of the tool surface lowers its hardness and increases  wear.
45
46
The content of the paper is the following. In the next sections the main features of the DEM used are described. Details of different models derived to treat the frictional contact conditions are given.
47
48
In the second part of the paper the formulation of volumetric locking-free linear triangles and tetrahedra is presented using a FIC approach. Next an algorithm for dynamic analysis of geomechanics problems combining DEM and locking-free FEM is described. Examples of the applications of coupled DEM-FEM formulation to a number of geomechanics problems are presented.
49
50
==2 DISCRETE ELEMENT METHOD FORMULATION==
51
52
The numerical model of rock/soil media adopted in the present study is based on the spherical discrete element method (DEM) which is widely recognized as a suitable tool to model geomaterials <span id='citeF-6'></span><span id='citeF-7'></span><span id='citeF-8'></span><span id='citeF-9'></span>[[#cite-6|[6,7,8,9]]]. Formulation of spherical discrete elements (called also distinct elements)  following the main assumptions of Cundall <span id='citeF-6'></span><span id='citeF-10'></span>[[#cite-6|[6,10]]] has been developed by Rojek et al. in <span id='citeF-11'></span>[[#cite-11|[11]]] and implemented in an explicit dynamic formulation. Within the DEM, it is assumed that the solid material  can be represented as a collection of rigid particles (spheres or balls in 3D and discs in 2D) interacting among themselves in the normal and tangential directions.  Material deformation is assumed to be concentrated at the contact points. Appropriate contact laws allow us to obtain desired macroscopic material properties. The contact law used for rock modelling takes into account cohesive bonds between rigid particles. The contact law can be seen as the formulation of the material model at the microscopic level. Cohesive bonds can be broken, thus allowing  to simulate fracture of material and its propagation.
53
54
===2.1 Equations of motion===
55
56
The translational and rotational motion of rigid spherical or cylindrical elements (particles) (Fig. [[#img-1|1]]) is described by means of the standard equations of rigid body dynamics. For the <math display="inline">i</math>-th element we have
57
58
<span id="eq-1"></span>
59
{| class="formulaSCP" style="width: 100%; text-align: left;" 
60
|-
61
| 
62
{| style="text-align: left; margin:auto;width: 100%;" 
63
|-
64
| style="text-align: center;" | <math>m_i \ddot{\textbf{u}}_i= {\mathbf{F}}_i\,, </math>
65
|}
66
| style="width: 5px;text-align: right;white-space: nowrap;" | (1)
67
|}
68
69
<span id="eq-2"></span>
70
{| class="formulaSCP" style="width: 100%; text-align: left;" 
71
|-
72
| 
73
{| style="text-align: left; margin:auto;width: 100%;" 
74
|-
75
| style="text-align: center;" | <math>I_i \dot{{\omega }}_i= {\mathbf{T}}_i\,, </math>
76
|}
77
| style="width: 5px;text-align: right;white-space: nowrap;" | (2)
78
|}
79
80
where <math>{u}</math>is the element centroid displacement in a fixed (inertial) coordinate frame <math>\mathbf{X}</math>, <math>\omega </math>&#8211; the angular velocity, <math display="inline">m</math> &#8211; the element mass, <math display="inline">I</math> &#8211; the moment of inertia, <math>\mathbf{F}</math>&#8211; the resultant force, and <math>\mathbf{T}</math>&#8211; the resultant moment about the central axes. Vectors <math>\mathbf{F}</math> and <math>\mathbf{T}</math> are sums of all forces and moments applied to the <math display="inline">i</math>-th element due to external load, contact interactions with neighbouring spheres and other obstacles, as well as forces resulting from damping in the system. The form of the rotational equation ([[#eq-2|2]]) is valid for spheres and cylinders (in 2D) and is simplified with respect to a general form for an arbitrary rigid body with the rotational inertial properties represented by a second order tensor. In the general case it is more convenient to describe the rotational motion with respect to a co-rotational frame <math>\mathbf{x}</math> which is embedded at each element since in this frame the tensor of inertia is constant. The tensor of inertia for a sphere or cylinder (in 2D) does not change in the fixed global coordinate system <math>\mathbf{X}</math>, and hence  the rotational motion can be easily described in this system.
81
82
<div id='img-1'></div>
83
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
84
|-
85
|[[Image:Draft_Samper_102318540-motion.png|396px|Motion of a rigid particle]]
86
|- style="text-align: center; font-size: 75%;"
87
| colspan="1" | '''Figure 1:''' Motion of a rigid particle
88
|}
89
90
Equations of motion ([[#eq-1|1]]) and ([[#eq-2|2]]) are integrated in time using a central difference scheme. The time integration operator for the translational motion at the <math display="inline">n</math>-th time step is as follows:
91
92
<span id="eq-3"></span>
93
{| class="formulaSCP" style="width: 100%; text-align: left;" 
94
|-
95
| 
96
{| style="text-align: left; margin:auto;width: 100%;" 
97
|-
98
| style="text-align: center;" | <math>\ddot{{\textbf u}}_i^{n}=\frac{{{\textbf F}}_i^n}{m_i}\,,  </math>
99
|}
100
| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
101
|}
102
103
<span id="eq-4"></span>
104
{| class="formulaSCP" style="width: 100%; text-align: left;" 
105
|-
106
| 
107
{| style="text-align: left; margin:auto;width: 100%;" 
108
|-
109
| style="text-align: center;" | <math>\dot{{\textbf u}}_i^{n+1/2}=\dot{{\textbf u}}_i^{n-1/2}+\ddot{{\textbf u}}_i^{n}{\Delta t}\,,  </math>
110
|}
111
| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
112
|}
113
114
<span id="eq-5"></span>
115
{| class="formulaSCP" style="width: 100%; text-align: left;" 
116
|-
117
| 
118
{| style="text-align: left; margin:auto;width: 100%;" 
119
|-
120
| style="text-align: center;" | <math>{{\textbf u}}_i^{n+1}={{\textbf u}}_i^{n}+\dot{{\textbf u}}_i^{n+1/2}{\Delta t}\,. </math>
121
|}
122
| style="width: 5px;text-align: right;white-space: nowrap;" | (5)
123
|}
124
125
The first two steps in the integration scheme for the rotational motion are identical to those given by Eqs.([[#eq-3|3]]) and ([[#eq-4|4]]):
126
127
<span id="eq-6"></span>
128
{| class="formulaSCP" style="width: 100%; text-align: left;" 
129
|-
130
| 
131
{| style="text-align: left; margin:auto;width: 100%;" 
132
|-
133
| style="text-align: center;" | <math>\dot{{\omega }}_i^{n}=\frac{{{T}}_i^n}{I_i}\,, </math>
134
|}
135
| style="width: 5px;text-align: right;white-space: nowrap;" | (6)
136
|}
137
138
<span id="eq-7"></span>
139
{| class="formulaSCP" style="width: 100%; text-align: left;" 
140
|-
141
| 
142
{| style="text-align: left; margin:auto;width: 100%;" 
143
|-
144
| style="text-align: center;" | <math>{{\omega }}_i^{n+1/2}={{\omega }}_i^{n-1/2}+\dot{{\omega }}_i^{n}{\Delta t}\,. </math>
145
|}
146
| style="width: 5px;text-align: right;white-space: nowrap;" | (7)
147
|}
148
149
The vector of incremental rotation <math display="inline">{\Delta \theta }=\{ \Delta \theta _x\,\Delta \theta _y\,\Delta \theta _z\} ^{\mathrm{T}}</math> is calculated as
150
151
<span id="eq-8"></span>
152
{| class="formulaSCP" style="width: 100%; text-align: left;" 
153
|-
154
| 
155
{| style="text-align: left; margin:auto;width: 100%;" 
156
|-
157
| style="text-align: center;" | <math>{{\Delta \theta }}_i={{\omega }}_i^{n+1/2}{\Delta t}\,, </math>
158
|}
159
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
160
|}
161
162
Knowledge of the incremental rotation suffices to update the tangential contact forces. It is also possible  to track the rotational position of particles, if  necessary <span id='citeF-12'></span>[[#cite-12|[12]]]. Then the rotation matrices between the moving frames embedded in  the particles and the fixed global frame must be updated incrementally using an adequate multiplicative  scheme <span id='citeF-11'></span>[[#cite-11|[11]]].
163
164
===2.2 Contact search algorithm===
165
166
Changing contact pairs of elements during the analysis process must be automatically detected. The simple approach to identify interaction pairs by checking every sphere against every other sphere would be very inefficient, as the computational time is proportional to <math display="inline">n^2</math>, where <math display="inline">n</math> is the number of elements. In our formulation the search is based on  quad-tree and oct-tree structures. In this case the computation time of the contact search is proportional to <math display="inline">n\ln{n}</math>, which allows to solve large frictional contact systems.
167
168
===2.3 Evaluation of contact forces===
169
170
====2.3.1 <span id='lb-2.3.1'></span>Decomposition ~of ~the ~contact ~force====
171
172
Once  contact between a pair of elements has been detected, the forces occurring at the contact point are calculated. The interaction between the two interacting bodies can be represented by the contact forces <math display="inline">{F}_1</math> and <math display="inline">{F}_2</math>, which by the Newton's third law satisfy the following relation:
173
174
{| class="formulaSCP" style="width: 100%; text-align: left;" 
175
|-
176
| 
177
{| style="text-align: left; margin:auto;width: 100%;" 
178
|-
179
| style="text-align: center;" | <math>{F}_{1}=-{F}_2\,. </math>
180
|}
181
| style="width: 5px;text-align: right;white-space: nowrap;" | (9)
182
|}
183
184
We take <math display="inline">{F}={F}_1</math> and decompose <math>{F}</math>into the normal and tangential components, <math display="inline">{F}_n</math>and <math display="inline">{F}_T</math>, respectively (Fig. [[#img-2|2]])
185
186
{| class="formulaSCP" style="width: 100%; text-align: left;" 
187
|-
188
| 
189
{| style="text-align: left; margin:auto;width: 100%;" 
190
|-
191
| style="text-align: center;" | <math> {F}={F}_n+{F}_T =F_n{n}+{F}_T \,, </math>
192
|}
193
|}
194
195
where <math>{n}</math>is the unit vector normal to the particle surface at the contact point (this implies that it lies along the line connecting the centers of the two particles) and directed outwards from the particle 1.
196
197
The contact forces <math display="inline">F_n</math> and <math display="inline">{F}_T</math> are obtained using a constitutive model formulated for the contact between two rigid spheres (or discs in 2D) (Fig. [[#img-3|3]]). The contact interface in our formulation  is characterized by the normal and tangential stiffness <math display="inline">k_n</math> and <math display="inline">k_T</math>, respectively, the Coulomb friction coefficient <math display="inline">\mu </math>, and the contact damping coefficient <math display="inline">c_n</math>.
198
199
<div id='img-2'></div>
200
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
201
|-
202
|[[Image:Draft_Samper_102318540-towspheres.png|246px|Decomposition of the contact force into the normal and tangential components]]
203
|- style="text-align: center; font-size: 75%;"
204
| colspan="1" | '''Figure 2:''' Decomposition of the contact force into the normal and tangential components
205
|}
206
207
<div id='img-3'></div>
208
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
209
|-
210
|[[Image:Draft_Samper_102318540-contact-point.png|384px|Model of the contact interface]]
211
|- style="text-align: center; font-size: 75%;"
212
| colspan="1" | '''Figure 3:''' Model of the contact interface
213
|}
214
215
====2.3.2 <span id='lb-2.3.2'></span>Normal ~contact ~force====
216
217
The normal contact force <math display="inline">F_n</math> is decomposed into the elastic part <math display="inline">F_{ne}</math> and the damping contact force <math display="inline">F_{nd}</math>
218
219
<span id="eq-10"></span>
220
{| class="formulaSCP" style="width: 100%; text-align: left;" 
221
|-
222
| 
223
{| style="text-align: left; margin:auto;width: 100%;" 
224
|-
225
| style="text-align: center;" | <math>F_n = F_{ne} + F_{nd}\,. </math>
226
|}
227
| style="width: 5px;text-align: right;white-space: nowrap;" | (10)
228
|}
229
230
The damping part is used  to decrease oscillations of the contact forces and to dissipate kinetic energy.
231
232
The elastic part of the normal contact force <math display="inline">F_{ne}</math> is proportional to the normal stiffness <math display="inline">k_n</math> and to the penetration of the two particle surfaces <math display="inline">u_{rn}</math>, i.e.
233
234
<span id="eq-11"></span>
235
{| class="formulaSCP" style="width: 100%; text-align: left;" 
236
|-
237
| 
238
{| style="text-align: left; margin:auto;width: 100%;" 
239
|-
240
| style="text-align: center;" | <math>F_{ne}=k_n u_{rn}\,. </math>
241
|}
242
| style="width: 5px;text-align: right;white-space: nowrap;" | (11)
243
|}
244
245
The penetration is calculated as
246
247
{| class="formulaSCP" style="width: 100%; text-align: left;" 
248
|-
249
| 
250
{| style="text-align: left; margin:auto;width: 100%;" 
251
|-
252
| style="text-align: center;" | <math>u_{rn} = d - r_1 - r_2\,, </math>
253
|}
254
| style="width: 5px;text-align: right;white-space: nowrap;" | (12)
255
|}
256
257
where <math display="inline">d</math> is the distance between the particle centres, and <math display="inline">r_1</math>, <math display="inline">r_2</math> their radiae. In the  present study no cohesion is allowed, so  no tensile normal contact forces are allowed and hence
258
259
<span id="eq-13"></span>
260
{| class="formulaSCP" style="width: 100%; text-align: left;" 
261
|-
262
| 
263
{| style="text-align: left; margin:auto;width: 100%;" 
264
|-
265
| style="text-align: center;" | <math>F_{ne}\le 0\,.   </math>
266
|}
267
| style="width: 5px;text-align: right;white-space: nowrap;" | (13)
268
|}
269
270
If <math display="inline">u_{rn}<0</math>, Eq. ([[#eq-11|11]]) holds, otherwise <math display="inline">F_{ne}= 0</math>.
271
272
The contact damping force is assumed to be of viscous type and given by
273
274
{| class="formulaSCP" style="width: 100%; text-align: left;" 
275
|-
276
| 
277
{| style="text-align: left; margin:auto;width: 100%;" 
278
|-
279
| style="text-align: center;" | <math>F_{nd}=c_n v_{rn} </math>
280
|}
281
| style="width: 5px;text-align: right;white-space: nowrap;" | (14)
282
|}
283
284
where <math display="inline">v_{rn}</math> is the normal relative velocity of the centres of the two particles in contact defined by
285
286
<span id="eq-15"></span>
287
{| class="formulaSCP" style="width: 100%; text-align: left;" 
288
|-
289
| 
290
{| style="text-align: left; margin:auto;width: 100%;" 
291
|-
292
| style="text-align: center;" | <math>{v}_{rn}=(\dot{{u}}_2-\dot{{u}}_1)\cdot { n}\,. </math>
293
|}
294
| style="width: 5px;text-align: right;white-space: nowrap;" | (15)
295
|}
296
297
The damping coefficient <math display="inline">c_n</math> can be taken as a fraction of the critical damping <math display="inline">C_{cr}</math> for the system of two rigid bodies with masses <math display="inline">m_1</math> and <math display="inline">m_2</math>, connected with a spring  of stiffness <math display="inline">k_n</math> (<span id='citeF-13'></span>[[#cite-13|[13]]]) with
298
299
<span id="eq-16"></span>
300
{| class="formulaSCP" style="width: 100%; text-align: left;" 
301
|-
302
| 
303
{| style="text-align: left; margin:auto;width: 100%;" 
304
|-
305
| style="text-align: center;" | <math>C_{cr}=2\sqrt{\frac{m_1m_2k_n}{m_1+m_2} }\,. </math>
306
|}
307
| style="width: 5px;text-align: right;white-space: nowrap;" | (16)
308
|}
309
310
====2.3.3 <span id='lb-2.3.3'></span>Tangential ~frictional ~contact====
311
312
In the absence of cohesion (if the particles were not bonded at all or  the initial cohesive bond has been broken) the tangential reaction <math display="inline">{F}_T</math>appears by friction opposing the relative motion at the contact point. The relative tangential velocity at the contact point <math display="inline">{v}_{rT}</math> is calculated from the following relationship
313
314
<span id="eq-17"></span>
315
{| class="formulaSCP" style="width: 100%; text-align: left;" 
316
|-
317
| 
318
{| style="text-align: left; margin:auto;width: 100%;" 
319
|-
320
| style="text-align: center;" | <math>{v}_{rT}= {v}_{r}-{v}_{r}\cdot {n}\,, </math>
321
|}
322
| style="width: 5px;text-align: right;white-space: nowrap;" | (17)
323
|}
324
325
{| class="formulaSCP" style="width: 100%; text-align: left;" 
326
|-
327
| 
328
{| style="text-align: left; margin:auto;width: 100%;" 
329
|-
330
| style="text-align: center;" | <math>with {v}_r=(\dot{{u}}_{2} +{\omega }_{2} \times {r}_{c2})- (\dot{{u}}_{1} +{\omega }_1 \times {r}_{c1})\,, </math>
331
|}
332
| style="width: 5px;text-align: right;white-space: nowrap;" | (18)
333
|}
334
335
where <math display="inline">\dot{{u}}_{1}</math>, <math display="inline">\dot{{u}}_{2}</math>, and <math display="inline">{\omega }_1</math>, <math display="inline">{\omega }_2</math> are the translational and rotational velocities of the particles, and <math display="inline">{ r}_{c1}</math> and <math display="inline">{r}_{c2}</math> are the vectors connecting particle centres with contact points.
336
337
<div id='img-4'></div>
338
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
339
|-
340
|[[Image:Draft_Samper_102318540-coulomb.png|462px|]]
341
|[[Image:Draft_Samper_102318540-coulombreg.png|462px|Friction force vs. relative tangential displacement a) Coulomb law, b) regularized Coulomb law]]
342
|- style="text-align: center; font-size: 75%;"
343
| colspan="2" | '''Figure 4:''' Friction force vs. relative tangential displacement a) Coulomb law, b) regularized Coulomb law
344
|}
345
346
The relationship between the friction force <math display="inline">\| {F}_T\| </math>and the relative tangential displacement <math display="inline">u_{rT}</math> for the classical Coulomb model (for a constant normal force <math display="inline">F_n</math>) is shown in Fig. [[#img-4|4]]a. This relationship would produce non physical oscillations of the friction force in the numerical solution due to possible changes of the direction of sliding velocity. To prevent this, the Coulomb friction model must be regularized. The regularization procedure chosen involves decomposition of the tangential relative velocity into  reversible and irreversible parts, <math display="inline">{v}_{rT}^{\mathrm{r}}</math> and <math display="inline">{v}_{rT}^{\mathrm{ir}}</math>, respectively as:
347
348
{| class="formulaSCP" style="width: 100%; text-align: left;" 
349
|-
350
| 
351
{| style="text-align: left; margin:auto;width: 100%;" 
352
|-
353
| style="text-align: center;" | <math>{v}_{rT}= {v}_{rT}^{\mathrm{r}}+{v}_{r}^{\mathrm{ir}}\,. </math>
354
|}
355
| style="width: 5px;text-align: right;white-space: nowrap;" | (19)
356
|}
357
358
This is equivalent to formulating the frictional contact as a problem analogous to that of elastoplasticity, which can be seen clearly from the friction force-tangential displacement relationship in Fig. [[#img-4|4]]b. This analogy allows us to calculate the friction force employing the standard radial return algorithm []. First a trial state is calculated
359
360
<span id="eq-20"></span>
361
{| class="formulaSCP" style="width: 100%; text-align: left;" 
362
|-
363
| 
364
{| style="text-align: left; margin:auto;width: 100%;" 
365
|-
366
| style="text-align: center;" | <math>{F}_T^{\mathrm{trial}}={F}_T^{\mathrm{old}}-k_T {v}_{rT} \Delta t\,,   </math>
367
|}
368
| style="width: 5px;text-align: right;white-space: nowrap;" | (20)
369
|}
370
371
and then the slip condition is checked
372
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>\phi ^{\mathrm{trial}}=\| {F}_T^{\mathrm{trial}}\|{-\mu}|F_n|\,. </math>
379
|}
380
| style="width: 5px;text-align: right;white-space: nowrap;" | (21)
381
|}
382
383
If <math display="inline">\phi ^{\mathrm{trial}}\le{0}</math>, we have stick contact and the friction force is assigned the trial value
384
385
{| class="formulaSCP" style="width: 100%; text-align: left;" 
386
|-
387
| 
388
{| style="text-align: left; margin:auto;width: 100%;" 
389
|-
390
| style="text-align: center;" | <math>{F}_T^{\mathrm{new}}={F}_T^{\mathrm{trial}}\,, </math>
391
|}
392
| style="width: 5px;text-align: right;white-space: nowrap;" | (22)
393
|}
394
395
otherwise (slip contact) a return mapping is performed giving
396
397
{| class="formulaSCP" style="width: 100%; text-align: left;" 
398
|-
399
| 
400
{| style="text-align: left; margin:auto;width: 100%;" 
401
|-
402
| style="text-align: center;" | <math>{F}_T^{\mathrm{new}}=\mu |F_n|\frac{{F}_T^{\mathrm{trial}}}{\| {F}_T^{\mathrm{trial}} \| }\,. </math>
403
|}
404
| style="width: 5px;text-align: right;white-space: nowrap;" | (23)
405
|}
406
407
===2.4 Background damping===
408
409
A quasi-static state of equilibrium for the assembly of particles can be achieved by application of adequate damping. The contact damping previously described is a function of the relative velocity of the contacting bodies. It is sometimes necessary to apply damping for non-contacting particles to dissipate their energy. We have considered two types of such damping of viscous  and non-viscous type, referred here as background damping. In both cases damping terms <math display="inline">{{F}}_i^{\mathrm{damp}}</math> and <math display="inline">{\mathbf{T}}_{i}^{\mathrm{damp}}</math> are  added to equations of motion ([[#eq-1|1]]) and ([[#eq-2|2]])
410
411
<span id="eq-24"></span>
412
{| class="formulaSCP" style="width: 100%; text-align: left;" 
413
|-
414
| 
415
{| style="text-align: left; margin:auto;width: 100%;" 
416
|-
417
| style="text-align: center;" | <math>m_i \ddot{{u}}_i= {{F}}_i+{{F}}_i^{\mathrm{damp}}\,, </math>
418
|}
419
| style="width: 5px;text-align: right;white-space: nowrap;" | (24)
420
|}
421
422
<span id="eq-25"></span>
423
{| class="formulaSCP" style="width: 100%; text-align: left;" 
424
|-
425
| 
426
{| style="text-align: left; margin:auto;width: 100%;" 
427
|-
428
| style="text-align: center;" | <math>I_i \dot{{\omega }}_i= {\mathbf{T}}_i+{\mathbf{T}}_i^{\mathrm{damp}}\,. </math>
429
|}
430
| style="width: 5px;text-align: right;white-space: nowrap;" | (25)
431
|}
432
433
where
434
435
*  for viscous damping
436
437
<span id="eq-26"></span>
438
{| class="formulaSCP" style="width: 100%; text-align: left;" 
439
|-
440
| 
441
{| style="text-align: left; margin:auto;width: 100%;" 
442
|-
443
| style="text-align: center;" | <math>
444
445
{{F}}_{i}^{\mathrm{damp}}=-\alpha ^{vt}m_{i}\dot{{ u}}_{i}\,,   </math>
446
|}
447
| style="width: 5px;text-align: right;white-space: nowrap;" | (26)
448
|}
449
450
<span id="eq-27"></span>
451
{| class="formulaSCP" style="width: 100%; text-align: left;" 
452
|-
453
| 
454
{| style="text-align: left; margin:auto;width: 100%;" 
455
|-
456
| style="text-align: center;" | <math>
457
458
{\mathbf{T}}_{i}^{\mathrm{damp}}=-\alpha ^{vr}I_{i}{\omega }_{i}\,,   </math>
459
|}
460
| style="width: 5px;text-align: right;white-space: nowrap;" | (27)
461
|}
462
463
*  for non-viscous damping
464
465
<span id="eq-28"></span>
466
{| class="formulaSCP" style="width: 100%; text-align: left;" 
467
|-
468
| 
469
{| style="text-align: left; margin:auto;width: 100%;" 
470
|-
471
| style="text-align: center;" | <math>
472
473
{{F}}_{i}^{\mathrm{damp}}=-\alpha ^{nvt}\Vert {{F}}_{i}\Vert \frac{\dot{{u}}_{i}}{\Vert \dot{{u}}_{i}\Vert }\,,   </math>
474
|}
475
| style="width: 5px;text-align: right;white-space: nowrap;" | (28)
476
|}
477
478
<span id="eq-29"></span>
479
{| class="formulaSCP" style="width: 100%; text-align: left;" 
480
|-
481
| 
482
{| style="text-align: left; margin:auto;width: 100%;" 
483
|-
484
| style="text-align: center;" | <math>
485
486
{\mathbf{T}}_{i}^{\mathrm{damp}}=-\alpha ^{nvr}\Vert {\mathbf{T}}_{i}\Vert \frac{{\omega }_{i}}{\Vert {\omega }_{i}\Vert }\,.   </math>
487
|}
488
| style="width: 5px;text-align: right;white-space: nowrap;" | (29)
489
|}
490
491
where <math display="inline">\alpha ^{vt}</math>, <math display="inline">\alpha ^{vr}</math>, <math display="inline">\alpha ^{nvt}</math> and <math display="inline">\alpha ^{nvr}</math> are damping constants. It can be seen from Eqs. ([[#eq-26|26]])&#8211;([[#eq-29|29]]) that both the non-viscous and viscous damping terms are opposite to the velocity and the difference lays in the evaluation of  the damping force. Viscous damping is proportional to the magnitude of velocity, while non-viscous damping is proportional to the magnitude of the resultant force and the bending moment.
492
493
===2.5 Numerical stability===
494
495
Explicit integration in time yields high computational efficiency and it enables the solution of large models. The known disadvantage of the explicit integration scheme is its conditional numerical stability imposing the limitation on the time step <math display="inline">{\Delta t}</math>, i.e.
496
497
{| class="formulaSCP" style="width: 100%; text-align: left;" 
498
|-
499
| 
500
{| style="text-align: left; margin:auto;width: 100%;" 
501
|-
502
| style="text-align: center;" | <math>{\Delta t} \le {\Delta t}_{cr} </math>
503
|}
504
| style="width: 5px;text-align: right;white-space: nowrap;" | (30)
505
|}
506
507
where <math>{\Delta t}_{cr}</math> is a critical time step determined by the highest natural frequency of the system <math>\omega _{ max}</math>
508
509
<span id="eq-31"></span>
510
{| class="formulaSCP" style="width: 100%; text-align: left;" 
511
|-
512
| 
513
{| style="text-align: left; margin:auto;width: 100%;" 
514
|-
515
| style="text-align: center;" | <math>{\Delta t}_{cr}=\frac{2}{\omega _{max}}\,. </math>
516
|}
517
| style="width: 5px;text-align: right;white-space: nowrap;" | (31)
518
|}
519
520
If damping exists, the critical time increment is given by
521
522
<span id="eq-32"></span>
523
{| class="formulaSCP" style="width: 100%; text-align: left;" 
524
|-
525
| 
526
{| style="text-align: left; margin:auto;width: 100%;" 
527
|-
528
| style="text-align: center;" | <math>{\Delta t}_{cr}=\frac{2}{\omega _{max}}\left(\sqrt{1+\xi ^2}-\xi \right)\,, </math>
529
|}
530
| style="width: 5px;text-align: right;white-space: nowrap;" | (32)
531
|}
532
533
where <math display="inline">\xi </math> is the fraction of the critical damping corresponding to the highest frequency <math>\omega _{max}</math>. Exact calculation of the highest frequency <math>\omega _{max}</math>would require solution of the eigenvalue problem defined for the whole system of connected rigid particles. In an approximate solution procedure, an eigenvalue problem can be defined separately for every rigid particle using the linearized equations of motion
534
535
<span id="eq-33"></span>
536
{| class="formulaSCP" style="width: 100%; text-align: left;" 
537
|-
538
| 
539
{| style="text-align: left; margin:auto;width: 100%;" 
540
|-
541
| style="text-align: center;" | <math>{\mathbf{m}}_i \ddot{\mathbf{r}}_i+ {\mathbf{k}}_i {\mathbf{r}}_i= \mathbf{0}\,, </math>
542
|}
543
| style="width: 5px;text-align: right;white-space: nowrap;" | (33)
544
|}
545
546
where
547
548
<span id="eq-34"></span>
549
{| class="formulaSCP" style="width: 100%; text-align: left;" 
550
|-
551
| 
552
{| style="text-align: left; margin:auto;width: 100%;" 
553
|-
554
| style="text-align: center;" | <math>{\mathbf{m}}_i = \{ m_i\; m_i\; m_i\; I_i\; I_i\; I_i\} ^{\mathrm{T}}\;,\quad  {\mathbf{r}}_i = \{ (u_x)_i\; (u_y)_i\; (u_z)_i\; (\theta _x)_i\; (\theta _y)_i\; (\theta _z)_i\} ^{\mathrm{T}}\,, </math>
555
|}
556
| style="width: 5px;text-align: right;white-space: nowrap;" | (34)
557
|}
558
559
and <math display="inline">{\mathbf{k}}_i</math> is the stiffness matrix accounting for the contributions from all the penalty constraints active for the <math display="inline">i</math>-th particle. Equation ([[#eq-34|34]]) defines  vectors <math display="inline">{\mathbf{m}}_i</math> and <math display="inline">{\mathbf{r}}_i</math> for a spherical particle in a 3D space. For a cylindrical particle in a 2D  model the respective vectors are defined as
560
561
<span id="eq-35"></span>
562
{| class="formulaSCP" style="width: 100%; text-align: left;" 
563
|-
564
| 
565
{| style="text-align: left; margin:auto;width: 100%;" 
566
|-
567
| style="text-align: center;" | <math>{\mathbf{m}}_i = \{ m_i\; m_i\; I_i\} ^{\mathrm{T}}\;,\quad {\mathbf{r}}_i = \{  (u_x)_i\; (u_y)_i\;(\theta _z)_i\} ^{\mathrm{T}}\,. </math>
568
|}
569
| style="width: 5px;text-align: right;white-space: nowrap;" | (35)
570
|}
571
572
Equation ([[#eq-33|33]]) leads to the eigenproblem
573
574
<span id="eq-36"></span>
575
{| class="formulaSCP" style="width: 100%; text-align: left;" 
576
|-
577
| 
578
{| style="text-align: left; margin:auto;width: 100%;" 
579
|-
580
| style="text-align: center;" | <math>{\mathbf{k}}_i {\mathbf{r}}_i= \lambda _j {\mathbf{m}}_i {\mathbf{r}}_i\,, </math>
581
|}
582
| style="width: 5px;text-align: right;white-space: nowrap;" | (36)
583
|}
584
585
whose eigenvalues <math display="inline">\lambda _j</math> (<math display="inline">j=1,\ldots , 6</math> in 3D case, and <math display="inline">j=1,2,3</math> for 2D case) are the squared frequencies of the free vibration frequencies:
586
587
<span id="eq-37"></span>
588
{| class="formulaSCP" style="width: 100%; text-align: left;" 
589
|-
590
| 
591
{| style="text-align: left; margin:auto;width: 100%;" 
592
|-
593
| style="text-align: center;" | <math>\lambda _j=\omega _j^2\,. </math>
594
|}
595
| style="width: 5px;text-align: right;white-space: nowrap;" | (37)
596
|}
597
598
In a 3D problem, three of  six frequencies <math display="inline">\omega _j</math> are translational, and the other three &#8211; rotational.
599
600
In the algorithm implemented, a further simplification is assumed. The maximum frequency is estimated as the maximum of natural frequencies of mass&#8211;spring systems defined for all the particles with one translational and one rotational degree of freedom. The translational and rotational free vibrations are governed by the following equations:
601
602
<span id="eq-38"></span>
603
{| class="formulaSCP" style="width: 100%; text-align: left;" 
604
|-
605
| 
606
{| style="text-align: left; margin:auto;width: 100%;" 
607
|-
608
| style="text-align: center;" | <math>m_i \ddot{u}_n + k_n u_n=0\,, </math>
609
|}
610
| style="width: 5px;text-align: right;white-space: nowrap;" | (38)
611
|}
612
613
<span id="eq-39"></span>
614
{| class="formulaSCP" style="width: 100%; text-align: left;" 
615
|-
616
| 
617
{| style="text-align: left; margin:auto;width: 100%;" 
618
|-
619
| style="text-align: center;" | <math>I_i \ddot{\theta } + k_\theta \theta=0\,, </math>
620
|}
621
| style="width: 5px;text-align: right;white-space: nowrap;" | (39)
622
|}
623
624
where it is assumed that the translational motion <math display="inline">u_n</math>is due to the contact interaction along the normal direction (the spring stiffness <math display="inline">k_n</math> represents the penalty stiffness in the normal direction) and the rotational stiffness is due to the contact stiffness (penalty) in the tangential direction. Given the tangential penalty <math display="inline">k_T</math>, the rotational stiffness <math display="inline">k_\theta </math> can be obtained as
625
626
<span id="eq-40"></span>
627
{| class="formulaSCP" style="width: 100%; text-align: left;" 
628
|-
629
| 
630
{| style="text-align: left; margin:auto;width: 100%;" 
631
|-
632
| style="text-align: center;" | <math>k_{\theta }=k_T r^2\,, </math>
633
|}
634
| style="width: 5px;text-align: right;white-space: nowrap;" | (40)
635
|}
636
637
where <math display="inline">r</math> is the length of the vector connecting the mass centre to the contact point.
638
639
The natural frequency of the translational vibrations is given by:
640
641
<span id="eq-41"></span>
642
{| class="formulaSCP" style="width: 100%; text-align: left;" 
643
|-
644
| 
645
{| style="text-align: left; margin:auto;width: 100%;" 
646
|-
647
| style="text-align: center;" | <math>\omega _n=\sqrt{\frac{k_n}{m_i}}\,, </math>
648
|}
649
| style="width: 5px;text-align: right;white-space: nowrap;" | (41)
650
|}
651
652
while the rotational frequency <math display="inline">\omega _{\theta }</math> can be obtained from
653
654
<span id="eq-42"></span>
655
{| class="formulaSCP" style="width: 100%; text-align: left;" 
656
|-
657
| 
658
{| style="text-align: left; margin:auto;width: 100%;" 
659
|-
660
| style="text-align: center;" | <math>\omega _{\theta }=\sqrt{\frac{k_{\theta }}{I_i}}\,. </math>
661
|}
662
| style="width: 5px;text-align: right;white-space: nowrap;" | (42)
663
|}
664
665
where <math display="inline">I_i</math> is the rotational inertia of a sphere
666
667
{| class="formulaSCP" style="width: 100%; text-align: left;" 
668
|-
669
| 
670
{| style="text-align: left; margin:auto;width: 100%;" 
671
|-
672
| style="text-align: center;" | <math>I= \frac{2}{5} mr^2 </math>
673
|}
674
| style="width: 5px;text-align: right;white-space: nowrap;" | (43)
675
|}
676
677
and <math display="inline">k_\theta </math> is given by Eq. ([[#eq-40|40]]). Substituting these expressions into Eq. ([[#eq-42|42]]) gives  the rotational frequency  as
678
679
{| class="formulaSCP" style="width: 100%; text-align: left;" 
680
|-
681
| 
682
{| style="text-align: left; margin:auto;width: 100%;" 
683
|-
684
| style="text-align: center;" | <math>\omega _{\theta }=\sqrt{\frac{5 k_{T}}{2 m_i}}\,. </math>
685
|}
686
| style="width: 5px;text-align: right;white-space: nowrap;" | (44)
687
|}
688
689
If <math display="inline">k_T=k_n</math>, the rotational frequency <math display="inline">\omega _{\theta }</math> is considerably higher than the translational frequency <math display="inline">\omega _n</math> obtained from Eq. ([[#eq-41|41]]), which results in a smaller critical time increment (see Eq. ([[#eq-31|31]])). To reduce the influence of the rotational frequencies in the value of the critical time step, the rotational inertia terms are scaled adequately. The concept of scaling rotational inertia terms is commonly used for shell elements <span id='citeF-14'></span>[[#cite-14|[14]]].
690
691
==3 MICROMECHANICAL CONSTITUTIVE MODELS==
692
693
===3.1 Elastic perfectly brittle model===
694
695
The elastic perfectly brittle contact model is characterized by linear elastic behaviour when  cohesive bonds are active. An instantaneous breakage of these bonds occurs when the interface strength is exceeded. When two particles are bonded the contact forces in both normal and tangential directions are calculated from the linear constitutive relationships:
696
697
<span id="eq-45"></span>
698
<span id="eq-46"></span>
699
{| class="formulaSCP" style="width: 100%; text-align: left;" 
700
|-
701
| 
702
{| style="text-align: left; margin:auto;width: 100%;" 
703
|-
704
| style="text-align: center;" | <math>\sigma \!=\!k_n u_n\,, </math>
705
| style="width: 5px;text-align: right;white-space: nowrap;" | (45)
706
|-
707
| style="text-align: center;" | <math> \tau \!=\!k_t \, u_t\,,  </math>
708
| style="width: 5px;text-align: right;white-space: nowrap;" | (46)
709
|}
710
|}
711
712
where <math display="inline">\sigma </math> and <math display="inline">\tau </math> are the  normal and tangential contact force, respectively, <math display="inline">k_n</math> and <math display="inline">k_t</math> are the interface stiffness in the normal and tangential directions and <math display="inline">u_n</math> and <math display="inline">u_t</math> the normal and tangential relative displacements, respectively.
713
714
Cohesive bonds are broken instantaneously when the interface strength is exceeded in the tangential direction by the tangential contact force or in the normal direction by the tensile contact force. The failure (decohesion) criterion is written (for 2D) as:
715
716
<span id="eq-47"></span>
717
<span id="eq-48"></span>
718
{| class="formulaSCP" style="width: 100%; text-align: left;" 
719
|-
720
| 
721
{| style="text-align: left; margin:auto;width: 100%;" 
722
|-
723
| style="text-align: center;" | <math>\sigma \!\le \!R_n\,, </math>
724
| style="width: 5px;text-align: right;white-space: nowrap;" | (47)
725
|-
726
| style="text-align: center;" | <math> |\tau | \!\le \!R_t\,,  </math>
727
| style="width: 5px;text-align: right;white-space: nowrap;" | (48)
728
|}
729
|}
730
731
where <math display="inline">R_n</math> and <math display="inline">R_t</math> are the interface strengths in the normal and tangential directions, respectively.
732
733
In the absence of cohesion the normal contact force can be compressive only, i.e.
734
735
<span id="eq-49"></span>
736
{| class="formulaSCP" style="width: 100%; text-align: left;" 
737
|-
738
| 
739
{| style="text-align: left; margin:auto;width: 100%;" 
740
|-
741
| style="text-align: center;" | <math>\sigma \le  0  </math>
742
|}
743
| style="width: 5px;text-align: right;white-space: nowrap;" | (49)
744
|}
745
746
and the (positive) tangential contact force is given by
747
748
<span id="eq-50"></span>
749
{| class="formulaSCP" style="width: 100%; text-align: left;" 
750
|-
751
| 
752
{| style="text-align: left; margin:auto;width: 100%;" 
753
|-
754
| style="text-align: center;" | <math>\tau  = \mu |\sigma |  </math>
755
|}
756
| style="width: 5px;text-align: right;white-space: nowrap;" | (50)
757
|}
758
759
if <math display="inline">\sigma{<0}</math>  or zero otherwise. The friction force is given by Eq. ([[#eq-50|50]]) expressing the Coulomb friction law, with <math display="inline">\mu </math>  being the Coulomb friction coefficient.
760
761
<div id='img-5'></div>
762
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
763
|-
764
|[[Image:Draft_Samper_102318540-brittlen.png|312px|Normal contact force in the elastic perfectly brittle model]]
765
|- style="text-align: center; font-size: 75%;"
766
| colspan="1" | '''Figure 5:''' Normal contact force in the elastic perfectly brittle model
767
|}
768
769
<div id='img-6'></div>
770
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
771
|-
772
|[[Image:Draft_Samper_102318540-brittlet.png|312px|Tangential contact force in the elastic perfectly brittle model (tensile normal contact force)]]
773
|- style="text-align: center; font-size: 75%;"
774
| colspan="1" | '''Figure 6:''' Tangential contact force in the elastic perfectly brittle model (tensile normal contact force)
775
|}
776
777
<div id='img-7'></div>
778
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
779
|-
780
|[[Image:Draft_Samper_102318540-brittlesurf.png|273px|Failure surface for the elastic perfectly brittle model]]
781
|- style="text-align: center; font-size: 75%;"
782
| colspan="1" | '''Figure 7:''' Failure surface for the elastic perfectly brittle model
783
|}
784
785
Contact laws for the normal and tangential directions for the elastic perfectly     brittle model are shown in Figs. [[#img-5|5]] and [[#img-6|6]], respectively. The failure surface for     the elastic perfectly brittle model defined by conditions ([[#eq-47|47]]) and ([[#eq-48|48]])      is shown in Fig.  [[#img-7|7]].
786
787
===3.2 Elasto-plastic contact model with linear softening===
788
789
A contact law relates the contact force acting between two spheres to their relative displacement. The model developed is based on the analogous elasto-plastic relationships established for the normal and tangential directions:
790
791
<span id="eq-51"></span>
792
{| class="formulaSCP" style="width: 100%; text-align: left;" 
793
|-
794
| 
795
{| style="text-align: left; margin:auto;width: 100%;" 
796
|-
797
| style="text-align: center;" | <math>\sigma _i=k_i(u_i-u_i^p)\,, \quad i=n,t  </math>
798
|}
799
| style="width: 5px;text-align: right;white-space: nowrap;" | (51)
800
|}
801
802
where <math display="inline">\sigma _i</math> is the normal or tangential contact force, <math display="inline">u_i</math> the total relative displacement, <math display="inline">u_i^p</math> the plastic component of the relative displacement, <math display="inline">k_i</math> the elastic parameters which characterize the contact interface stiffness and <math display="inline">n, t</math> the  normal and tangential directions.
803
804
A general yield criterion has been assumed in the following form:
805
806
<span id="eq-52"></span>
807
{| class="formulaSCP" style="width: 100%; text-align: left;" 
808
|-
809
| 
810
{| style="text-align: left; margin:auto;width: 100%;" 
811
|-
812
| style="text-align: center;" | <math>f(\sigma ,\alpha )=|\sigma | -[\sigma _Y -H \alpha ]\,,  </math>
813
|}
814
| style="width: 5px;text-align: right;white-space: nowrap;" | (52)
815
|}
816
817
where <math display="inline">H</math>  is  the plastic softening modulus (assumed to be positive) and <math display="inline">\alpha </math>  is  the softening  parameter represented by the plastic part of the relative displacement <math display="inline">u^p</math>.
818
819
An elasto-plastic contact law with linear softening is assumed for the shear and normal tensile contact forces, cf. Figs. [[#img-8|8]] and [[#img-9|9]]. An elastic linear law is assumed for compression.
820
821
Dependence of the tangential force on the normal force is expressed by the application of the 2D failure criterion shown in Fig. [[#img-10|10]].  Yielding of contact bonds can occur under a combined tension and shear  action (if the normal contact force is tensile) or due to shear only (if the normal contact force is compressive). After the contact bonds are broken due to yielding, standard frictional contact can occur between the spherical elements.
822
823
''Flow rule''
824
825
The plastic strain increments for both  formulations (i.e. the normal and tangential motions) are calculated from the following flow rule:
826
827
<span id="eq-53"></span>
828
{| class="formulaSCP" style="width: 100%; text-align: left;" 
829
|-
830
| 
831
{| style="text-align: left; margin:auto;width: 100%;" 
832
|-
833
| style="text-align: center;" | <math>\dot{u}^p_i=\gamma \, \mbox{sign}(\sigma _i) \,,  </math>
834
|}
835
| style="width: 5px;text-align: right;white-space: nowrap;" | (53)
836
|}
837
838
where <math display="inline">\gamma </math> is the absolute value of the slip rate (<math display="inline">\gamma{>0}</math>), <math display="inline">\dot{u}^p_i</math> the derivative of the plastic relative displacement (normal or tangential) and <math display="inline">\sigma _i</math> the contact force in the normal or tangential direction.
839
840
<div id='img-8'></div>
841
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
842
|-
843
|[[Image:Draft_Samper_102318540-elpln.png|486px|Elasto-plastic contact law with softening for the normal direction]]
844
|- style="text-align: center; font-size: 75%;"
845
| colspan="1" | '''Figure 8:''' Elasto-plastic contact law with softening for the normal direction
846
|}
847
848
<div id='img-9'></div>
849
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
850
|-
851
|[[Image:Draft_Samper_102318540-elplt.png|522px|Elasto-plastic contact law with softening for the tangential direction]]
852
|- style="text-align: center; font-size: 75%;"
853
| colspan="1" | '''Figure 9:''' Elasto-plastic contact law with softening for the tangential direction
854
|}
855
856
<div id='img-10'></div>
857
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
858
|-
859
|[[Image:Draft_Samper_102318540-yielsurf.png|231px|Yield surface for elasto-plastic model]]
860
|- style="text-align: center; font-size: 75%;"
861
| colspan="1" | '''Figure 10:''' Yield surface for elasto-plastic model
862
|}
863
864
===3.3 Contact model with elastic damage===
865
866
The elastic perfectly brittle model can be easily extended to account for elastic damage by assuming a softening behaviour  defined by a softening modulus <math display="inline">H</math>  introduced into the force-displacement relationship (Fig. [[#img-11|11]]). <div id='img-11'></div>
867
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
868
|-
869
|[[Image:Draft_Samper_102318540-damagen.png|486px|Normal contact force in the contact model with elastic damage]]
870
|- style="text-align: center; font-size: 75%;"
871
| colspan="1" | '''Figure 11:''' Normal contact force in the contact model with elastic damage
872
|}
873
874
The constitutive relationship for 1D elastic damage is given by:
875
876
<span id="eq-54"></span>
877
{| class="formulaSCP" style="width: 100%; text-align: left;" 
878
|-
879
| 
880
{| style="text-align: left; margin:auto;width: 100%;" 
881
|-
882
| style="text-align: center;" | <math>\sigma =k_n^D \,u_n =(1-\omega )k_n u_n  </math>
883
|}
884
| style="width: 5px;text-align: right;white-space: nowrap;" | (54)
885
|}
886
887
where <math display="inline">k_n^D </math> is the elastic damaged secant modulus and <math display="inline">\omega </math> the scalar  damage variable.
888
889
The scalar damage variable <math display="inline">\omega </math> is a measure of material damage. For the undamaged state <math display="inline">\omega=0</math> and for a damaged state <math display="inline">0<\omega \le 1</math>. The scalar damage variable <math display="inline">\omega </math> can be written in the following form:
890
891
<span id="eq-55"></span>
892
{| class="formulaSCP" style="width: 100%; text-align: left;" 
893
|-
894
| 
895
{| style="text-align: left; margin:auto;width: 100%;" 
896
|-
897
| style="text-align: center;" | <math>\omega = \frac{\psi (u_n)-1}{\psi (u_n)}  </math>
898
|}
899
| style="width: 5px;text-align: right;white-space: nowrap;" | (55)
900
|}
901
902
where <math display="inline">\psi (u_n)</math>  is a function of total relative displacement. For linear strain softening <math display="inline">\psi (u_n)</math>   is defined by
903
904
<span id="eq-56"></span>
905
{| class="formulaSCP" style="width: 100%; text-align: left;" 
906
|-
907
| 
908
{| style="text-align: left; margin:auto;width: 100%;" 
909
|-
910
| style="text-align: center;" | <math>\psi (u_n)=\left\{ \begin{array}{lll}1 & {for} & \displaystyle u_n\le \frac{R_n}{k_n}\\[2ex] \displaystyle \frac{k^2_n u_n}{HR_n +k_n R_n-Hk_n u_n} & {for}                   & \displaystyle \frac{R_n}{k_n}\le u_n\le \frac{R_n}{k_n}+\frac{R_n}{H}\\[2ex] \displaystyle \infty  & {for}                   & \displaystyle u_n\ge \frac{R_n}{k_n}+\frac{R_n}{H} \end{array} \right.  </math>
911
|}
912
| style="width: 5px;text-align: right;white-space: nowrap;" | (56)
913
|}
914
915
where <math display="inline">R_n</math> is the initial tensile strength and <math display="inline">H</math> is the softening modulus (taken as positive).
916
917
A similar contact force&#8211;displacement law with damage (Fig. [[#img-12|12]]) can be introduced for the tangential direction. Then the bond decohesion can occur either due to tension or shear.
918
919
<div id='img-12'></div>
920
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
921
|-
922
|[[Image:Draft_Samper_102318540-damaget.png|546px|Tangential contact force in the contact model with elastic damage]]
923
|- style="text-align: center; font-size: 75%;"
924
| colspan="1" | '''Figure 12:''' Tangential contact force in the contact model with elastic damage
925
|}
926
927
An optional failure criterion has been implemented  with a failure criterion based on the tensile contact force only. This criterion models better fracture of brittle materials, as the macroscopic failure is due to brittle rupture of atomic bonds in tension. This microscopic mechanism explains the macroscopic strain softening behaviour both under compression and tension states.
928
929
After the contact bonds are broken due to damage (<math display="inline">\omega=0</math>) standard frictional contact can occur between the spherical elements.
930
931
===3.4 Contact model with friction, wear and heat generation===
932
933
Frictional contact, wear and heat generation must be considered in order  to accurately model tool-rock/soil interaction. Cohesion is not included into the formulation of this model. Hence, the normal contact forces can be compressive only. The normal contact force <math display="inline">\sigma </math>  is calculated from the linear relationships:
934
935
<span id="eq-57"></span>
936
{| class="formulaSCP" style="width: 100%; text-align: left;" 
937
|-
938
| 
939
{| style="text-align: left; margin:auto;width: 100%;" 
940
|-
941
| style="text-align: center;" | <math>\sigma = k_n u_n  </math>
942
|}
943
| style="width: 5px;text-align: right;white-space: nowrap;" | (57)
944
|}
945
946
where <math display="inline">u_n\le 0</math>  is the normal relative displacement (penetration of one sphere against another). The frictional contact force is evaluated using the regularized Coulomb law
947
948
<span id="eq-58"></span>
949
{| class="formulaSCP" style="width: 100%; text-align: left;" 
950
|-
951
| 
952
{| style="text-align: left; margin:auto;width: 100%;" 
953
|-
954
| style="text-align: center;" | <math>\tau = \mu |\sigma |  </math>
955
|}
956
| style="width: 5px;text-align: right;white-space: nowrap;" | (58)
957
|}
958
959
where <math display="inline">\mu </math>  is the Coulomb friction coefficient. The frictional dissipation rate <math display="inline">\dot{D}</math> is calculated as
960
961
<span id="eq-59"></span>
962
{| class="formulaSCP" style="width: 100%; text-align: left;" 
963
|-
964
| 
965
{| style="text-align: left; margin:auto;width: 100%;" 
966
|-
967
| style="text-align: center;" | <math>\dot{D} = \tau v_t  </math>
968
|}
969
| style="width: 5px;text-align: right;white-space: nowrap;" | (59)
970
|}
971
972
where <math display="inline">v_t</math>  is the relative tangential velocity. The frictional dissipation is used to calculate the wear rate <math display="inline">\dot{w}</math> and the heat generated by friction <math display="inline">Q_{gen}</math>. This is computed by
973
974
<span id="eq-60"></span>
975
{| class="formulaSCP" style="width: 100%; text-align: left;" 
976
|-
977
| 
978
{| style="text-align: left; margin:auto;width: 100%;" 
979
|-
980
| style="text-align: center;" | <math>Q_{gen} = \chi \dot{D}  </math>
981
|}
982
| style="width: 5px;text-align: right;white-space: nowrap;" | (60)
983
|}
984
985
where <math display="inline">\chi </math> is the part of the friction work converted to heat (<math display="inline">\chi \le 1</math>). Heat generation through frictional dissipation is assumed to be absorbed equally by the two particles in contact.
986
987
Wear rate is calculated using the Archard law  <span id='citeF-15'></span>[[#cite-15|[15]]], which  assumes that the wear rate <math display="inline">\dot{w}</math> is proportional to the contact pressure <math display="inline">\sigma </math>   and to the slip velocity&nbsp;<math display="inline">v_T</math>
988
989
<span id="eq-61"></span>
990
{| class="formulaSCP" style="width: 100%; text-align: left;" 
991
|-
992
| 
993
{| style="text-align: left; margin:auto;width: 100%;" 
994
|-
995
| style="text-align: center;" | <math>\dot{w}=k \frac{\sigma v_T}{H}  </math>
996
|}
997
| style="width: 5px;text-align: right;white-space: nowrap;" | (61)
998
|}
999
1000
with <math display="inline">H</math> being the hardness of worn surface and <math display="inline">k</math> being a dimensionless parameter. The Archard law was derived originally for adhesive wear, the same form of law, however, can be obtained for abrasive wear, cf. <span id='citeF-16'></span>[[#cite-16|[16]]]. Values of adhesive and abrasive wear constants <math display="inline">k</math>  for different combinations of materials can be found in  <span id='citeF-16'></span>[[#cite-16|[16]]]. It is commonly accepted that wear is related to friction and thus friction coefficients are often introduced into Eq. ([[#eq-61|61]]). If the Coulomb friction law ([[#eq-58|58]]) holds,  the Archard wear law can be written in the following equivalent form:
1001
1002
<span id="eq-62"></span>
1003
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1004
|-
1005
| 
1006
{| style="text-align: left; margin:auto;width: 100%;" 
1007
|-
1008
| style="text-align: center;" | <math>\dot{w}=\bar{k} \frac{\tau v_T}{H}= \bar{k} \frac{\dot{D}}{H}  </math>
1009
|}
1010
| style="width: 5px;text-align: right;white-space: nowrap;" | (62)
1011
|}
1012
1013
in which the wear rate is related to frictional dissipation rate <math display="inline">\dot{D}</math>   and <math display="inline">\bar{k}=k/\mu </math>.
1014
1015
Influence of temperature on wear can be taken into account  by adaptation of the law of Archard given by Eq. ([[#eq-61|61]]). Thermal effects can be captured  approximately by taking the hardness <math display="inline">H</math> as a function of the temperature <math display="inline">T</math>
1016
1017
<span id="eq-63"></span>
1018
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1019
|-
1020
| 
1021
{| style="text-align: left; margin:auto;width: 100%;" 
1022
|-
1023
| style="text-align: center;" | <math>H=H(T)  </math>
1024
|}
1025
| style="width: 5px;text-align: right;white-space: nowrap;" | (63)
1026
|}
1027
1028
Then the Archard law ([[#eq-61|61]]) takes the following form:
1029
1030
<span id="eq-64"></span>
1031
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1032
|-
1033
| 
1034
{| style="text-align: left; margin:auto;width: 100%;" 
1035
|-
1036
| style="text-align: center;" | <math>\dot{w}=k \frac{\sigma v_T}{H(T)}  </math>
1037
|}
1038
| style="width: 5px;text-align: right;white-space: nowrap;" | (64)
1039
|}
1040
1041
Computation of temperature can be performed by solving thermal conduction problem on the discrete domain <span id='citeF-17'></span>[[#cite-17|[17]]].
1042
1043
==4 FORMULATION OF QUASI-INCOMPRESSIBLE SOLID FINITE ELEMENTS==
1044
1045
Many continuum finite elements exhibit so called “volumetric locking” in the analysis of incompressible or quasi-incompressible problems. Situations of this type are usual in the structural analysis of rubber materials, some geomechanical problems and most bulk metal forming processes. Volumetric locking is an undesirable effect leading to incorrect numerical results <span id='citeF-3'></span>[[#cite-3|[3]]].
1046
1047
Volumetric locking  is present in all low order solid elements based on the standard displacement formulation. The use of a mixed formulation or a selective integration technique eliminates the volumetric locking in many elements. These methods however, fail in some elements such as linear triangles and tetrahedra, due to lack of satisfaction of the Babuska-Brezzi conditions <span id='citeF-3'></span><span id='citeF-18'></span><span id='citeF-19'></span>[[#cite-3|[3,18,19]]] or alternatively the mixed patch test <span id='citeF-3'></span><span id='citeF-20'></span><span id='citeF-21'></span>[[#cite-3|[3,20,21]]] not being passed.
1048
1049
Considerable efforts have been made in recent years to develop linear triangles and tetrahedra producing correct (stable) results under incompressible situations. Brezzi and Pitkäranta <span id='citeF-22'></span>[[#cite-22|[22]]] proposed to extend the equation for the volumetric strain rate constraint for Stokes flows by adding a laplacian of pressure term. A similar method was derived for quasi-incompressible solids by Zienkiewicz and Taylor <span id='citeF-3'></span>[[#cite-3|[3]]]. Zienkiewicz ''et al.'' <span id='citeF-23'></span>[[#cite-23|[23]]] have proposed a stabilization technique which eliminates volumetric locking in incompressible solids based on a mixed formulation and a Characteristic Based Split (CBS) algorithm initially developed for fluids <span id='citeF-24'></span><span id='citeF-25'></span><span id='citeF-26'></span>[[#cite-24|[24,25,26]]] where a split of the pressure is introduced when solving the transient dynamic equations in time. Extensions of the CBS algorithm to solve bulk metal forming problems have been recently reported by Rojek ''et al.'' <span id='citeF-24'></span>[[#cite-24|[24]]]. Other methods to overcome volumetric locking are based on mixed displacement (or velocity)-pressure formulations using the Galerkin-Least-Square (GLS) method <span id='citeF-25'></span>[[#cite-25|[25]]], average nodal pressure  and average nodal deformation <span id='citeF-26'></span><span id='citeF-27'></span>[[#cite-26|[26,27]]] techniques, and Sub-Grid Scale (SGS) methods <span id='citeF-28'></span><span id='citeF-29'></span><span id='citeF-30'></span><span id='citeF-31'></span>[[#cite-28|[28,29,30,31]]].
1050
1051
In this paper volumetric locking is avoided by using a finite calculus (FIC) formulation. The basis of the FIC method is the satisfaction of the equations of balance of momentum and the constitutive equation relating the pressure with the volumetric strain in a domain of finite size.  The modified governing equations contain additional terms from standard infinitesimal theory. These terms introduce the necessary stability in the discretized equations to overcome the volumetric locking problem.
1052
1053
The FIC approach has been successfully used to derive stabilized finite element  and meshless methods for a wide range of advective-diffusive and fluid flow problems <span id='citeF-32'></span><span id='citeF-33'></span><span id='citeF-34'></span><span id='citeF-35'></span><span id='citeF-36'></span><span id='citeF-37'></span><span id='citeF-38'></span><span id='citeF-39'></span>[[#cite-32|[32,33,34,35,36,37,38,39]]].  The same ideas were applied in <span id='citeF-4'></span><span id='citeF-5'></span>[[#cite-4|[4,5]]]  to derive a stabilized formulation for quasi-incompressible and incompressible solids allowing the use of linear triangles and tetrahedra.
1054
1055
In the next section, the basis of the FIC method are given for static quasi-incompressible solid mechanics problems. The stabilized dynamic formulation for linear triangles and tetrahedra is presented and both semi-implicit and explicit monolithic solution schemes are described.
1056
1057
===4.1 Equilibrium equations===
1058
1059
The equilibrium equations for an elastic solid are written using the FIC technique as <span id='citeF-32'></span>[[#cite-32|[32]]]
1060
1061
<span id="eq-65"></span>
1062
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1063
|-
1064
| 
1065
{| style="text-align: left; margin:auto;width: 100%;" 
1066
|-
1067
| style="text-align: center;" | <math>{r}_{i}- \underline{\frac{{h}_k}{2}{\partial {r}_{i} \over \partial {x}_k}}\frac{{h}_k}{2}{\partial {r}_{i} \over \partial {x}_k}=0 \quad \mbox{in} \;\,\Omega \qquad k=1,n_d </math>
1068
|}
1069
| style="width: 5px;text-align: right;white-space: nowrap;" | (65)
1070
|}
1071
1072
where <math display="inline">n_d</math> is the number of space dimensions of the problems (i.e. <math display="inline">n_d =3</math> for 3D)
1073
1074
<span id="eq-66"></span>
1075
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1076
|-
1077
| 
1078
{| style="text-align: left; margin:auto;width: 100%;" 
1079
|-
1080
| style="text-align: center;" | <math>{r}_{i}: ={\partial {\sigma }_{ij} \over \partial {x}_j}+{b}_i </math>
1081
|}
1082
| style="width: 5px;text-align: right;white-space: nowrap;" | (66)
1083
|}
1084
1085
In  ([[#eq-65|65]]) and ([[#eq-66|66]]) <math>{\sigma }_{ij}</math>and <math>{b}_i</math>are the stresses and the body forces, respectively and <math>{h}_k</math>are characteristic length distances  of an arbitrary prismatic domain where equilibrium of forces is considered.
1086
1087
Equations ([[#eq-65|65]]) and ([[#eq-66|66]]) are completed with the boundary conditions on the displacements <math display="inline">u_i</math>
1088
1089
<span id="eq-67"></span>
1090
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1091
|-
1092
| 
1093
{| style="text-align: left; margin:auto;width: 100%;" 
1094
|-
1095
| style="text-align: center;" | <math>{u}_i-\bar{u}_i=0\quad \mbox{on} \;\,\Gamma _u </math>
1096
|}
1097
| style="width: 5px;text-align: right;white-space: nowrap;" | (67)
1098
|}
1099
1100
and the equilibrium of surface tractions
1101
1102
<span id="eq-68"></span>
1103
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1104
|-
1105
| 
1106
{| style="text-align: left; margin:auto;width: 100%;" 
1107
|-
1108
| style="text-align: center;" | <math>{\sigma }_{ij}{n}_j- \bar{t}_i- \underline{\frac{1}{2}{h}_k{n}_k{r}_{i}}\frac{1}{2}{h}_k{n}_k{r}_{i} =0\quad \mbox{on} \;\,\Gamma _t </math>
1109
|}
1110
| style="width: 5px;text-align: right;white-space: nowrap;" | (68)
1111
|}
1112
1113
In the above <math>\bar{u}_i</math>and <math>\bar{t}_i</math>are prescribed displacements and tractions over the boundaries <math display="inline">\Gamma _u</math>and <math display="inline">\Gamma _t</math>, respectively, <math display="inline">n_i</math> are the components of the unit normal vector and <math display="inline">h_k</math> are again the characteristic lengths.
1114
1115
The underlined terms in Eqs. ([[#eq-65|65]]) and ([[#eq-68|68]]) are a consequence of expressing the equilibrium of body forces and surface tractions in domains of finite size and retaining higher order terms than those usually accepted in the infinitesimal theory <span id='citeF-32'></span>[[#cite-32|[32]]].
1116
1117
These terms are essential to derive a stabilized finite element formulation as described below.
1118
1119
====4.1.1 Constitutive equations====
1120
1121
The stresses are split into deviatoric and volumetric (pressure) parts
1122
1123
<span id="eq-69"></span>
1124
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1125
|-
1126
| 
1127
{| style="text-align: left; margin:auto;width: 100%;" 
1128
|-
1129
| style="text-align: center;" | <math>{\sigma }_{ij}= s_{ij} +p{\delta }_{ij} </math>
1130
|}
1131
| style="width: 5px;text-align: right;white-space: nowrap;" | (69)
1132
|}
1133
1134
where <math>{\delta }_{ij}</math>is the Kronnecker delta function. The linear elastic constitutive equations for the deviatoric stresses <math display="inline">s_{ij}</math> are written as
1135
1136
<span id="eq-70"></span>
1137
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1138
|-
1139
| 
1140
{| style="text-align: left; margin:auto;width: 100%;" 
1141
|-
1142
| style="text-align: center;" | <math>s_{ij}\!\!=\!\! 2G \left({\varepsilon }_{ij}-\frac{1}{3}{\varepsilon }_{v}{\delta }_{ij}\right) </math>
1143
|}
1144
| style="width: 5px;text-align: right;white-space: nowrap;" | (70)
1145
|}
1146
1147
where <math display="inline">G</math>is the shear modulus,
1148
1149
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1150
|-
1151
| 
1152
{| style="text-align: left; margin:auto;width: 100%;" 
1153
|-
1154
| style="text-align: center;" | <math>{\varepsilon }_{ij}= \frac{1}{2} \left( {\partial {u}_i \over \partial {x}_j}+ {\partial {u}_j \over \partial x_i} \right) \quad  \mbox{and} \quad  {\varepsilon }_{v}={\varepsilon }_{ii}\,. </math>
1155
|}
1156
| style="width: 5px;text-align: right;white-space: nowrap;" | (71)
1157
|}
1158
1159
The constitutive equation for the pressure <math display="inline">p</math> can be written  for an arbitrary domain of finite size using the FIC formulation as [30]
1160
1161
<span id="eq-72"></span>
1162
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1163
|-
1164
| 
1165
{| style="text-align: left; margin:auto;width: 100%;" 
1166
|-
1167
| style="text-align: center;" | <math>\left({1\over K} p - {\varepsilon }_{v}\right)- \frac{{h}_k}{2}{\partial  \over \partial {x}_k} \left({1\over K} p - {\varepsilon }_{v}\right)=0\, ,\quad k=1,2,3\,\,\, \hbox{for }3D </math>
1168
|}
1169
| style="width: 5px;text-align: right;white-space: nowrap;" | (72)
1170
|}
1171
1172
where <math display="inline">K</math> is the bulk modulus of the material.
1173
1174
Note that for <math display="inline">h_k \to 0</math> the standard relationship between the pressure and the volumetric strain of the infinitesimal theory <math display="inline">(p=K\varepsilon _v )</math> is found.
1175
1176
For an incompressible material <math display="inline">K\to \infty </math> and Eq. ([[#eq-72|72]]) yields
1177
1178
<span id="eq-73"></span>
1179
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1180
|-
1181
| 
1182
{| style="text-align: left; margin:auto;width: 100%;" 
1183
|-
1184
| style="text-align: center;" | <math>\varepsilon _v - \frac{{h}_k}{2}{\partial {\varepsilon }_{v} \over \partial {x}_k}=0 </math>
1185
|}
1186
| style="width: 5px;text-align: right;white-space: nowrap;" | (73)
1187
|}
1188
1189
Eq. ([[#eq-73|73]]) expresses the limit incompressible behaviour of the solid.  This equation is typical in incompressible fluid dynamic problems and there arises from the mass continuity conditions <span id='citeF-32'></span><span id='citeF-33'></span>[[#cite-32|[32,33]]].
1190
1191
By combining Eqs. ([[#eq-65|65]]),  ([[#eq-66|66]]),  ([[#eq-69|69]]), ([[#eq-70|70]]) and ([[#eq-73|73]]) a mixed displacement&#8211;pressure formulation can be written as
1192
1193
<span id="eq-74"></span>
1194
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1195
|-
1196
| 
1197
{| style="text-align: left; margin:auto;width: 100%;" 
1198
|-
1199
| style="text-align: center;" | <math>{\partial s_{ij} \over \partial {x}_j}+{\partial p \over \partial x_i}+{b}_i-\frac{{h}_k}{2} {\partial r_i \over \partial {x}_k}=0 </math>
1200
|}
1201
| style="width: 5px;text-align: right;white-space: nowrap;" | (74)
1202
|}
1203
1204
<span id="eq-75"></span>
1205
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1206
|-
1207
| 
1208
{| style="text-align: left; margin:auto;width: 100%;" 
1209
|-
1210
| style="text-align: center;" | <math>\left(\frac{p}{K}-{\varepsilon }_{v}\right)-\frac{{h}_k}{2}{\partial  \over \partial {x}_k}\left(\frac{p}{K}-{\varepsilon }_{v}\right)=0 </math>
1211
|}
1212
| style="width: 5px;text-align: right;white-space: nowrap;" | (75)
1213
|}
1214
1215
Substituting Eq. ([[#eq-70|70]]) into ([[#eq-74|74]]) leads after some algebra to
1216
1217
<span id="eq-76"></span>
1218
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1219
|-
1220
| 
1221
{| style="text-align: left; margin:auto;width: 100%;" 
1222
|-
1223
| style="text-align: center;" | <math>{\partial {\varepsilon }_{v} \over \partial {x}_i}= \frac{3}{2G}\left[\hat{r}_{i}- \frac{{h}_k}{2}{\partial {r}_{i} \over \partial {x}_k} \right] </math>
1224
|}
1225
| style="width: 5px;text-align: right;white-space: nowrap;" | (76)
1226
|}
1227
1228
where <math>{r}_{i}</math>is defined by Eq. ([[#eq-66|66]]) and
1229
1230
<span id="eq-77"></span>
1231
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1232
|-
1233
| 
1234
{| style="text-align: left; margin:auto;width: 100%;" 
1235
|-
1236
| style="text-align: center;" | <math>\hat{r}_{i}=\frac{\partial }{\partial x_j}(2G \varepsilon _{ij})+{\partial p \over \partial x_i} + b_i </math>
1237
|}
1238
| style="width: 5px;text-align: right;white-space: nowrap;" | (77)
1239
|}
1240
1241
Substituting Eq. ([[#eq-76|76]]) into ([[#eq-75|75]]) gives after some simplification
1242
1243
<span id="eq-78"></span>
1244
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1245
|-
1246
| 
1247
{| style="text-align: left; margin:auto;width: 100%;" 
1248
|-
1249
| style="text-align: center;" | <math>\left(\frac{p}{K}-{\varepsilon }_{v}\right)- {h_i\over 2} \left({1\over K} {\partial p \over \partial x_i} - {3\over 2G}\hat r_i \right)- \left({3h_i\over 8}{{h}_k\over G}{\partial r_i \over \partial {x}_k}\right)=0 </math>
1250
|}
1251
| style="width: 5px;text-align: right;white-space: nowrap;" | (78)
1252
|}
1253
1254
<span id="eq-79"></span>
1255
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1256
|-
1257
| 
1258
{| style="text-align: left; margin:auto;width: 100%;" 
1259
|-
1260
| style="text-align: center;" | <math>{p\over K} - {\varepsilon }_{v}- \sum \limits _{i=1}^{n_d} \tau _i {\partial r_i \over \partial x_i}=0 </math>
1261
|}
1262
| style="width: 5px;text-align: right;white-space: nowrap;" | (79)
1263
|}
1264
1265
where
1266
1267
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1268
|-
1269
| 
1270
{| style="text-align: left; margin:auto;width: 100%;" 
1271
|-
1272
| style="text-align: center;" | <math>\tau _i = {3h_i^2\over 8G} </math>
1273
|}
1274
| style="width: 5px;text-align: right;white-space: nowrap;" | (80)
1275
|}
1276
1277
The coefficients <math display="inline">\tau _{i} </math> in Eq. ([[#eq-79|79]]) are also referred to as ''intrinsic time'' parameters per unit mass  (their dimensions are <math display="inline">{t^2 m^3\over Kg}</math> where <math display="inline">t</math> is the time).
1278
1279
===4.2 Non linear transient dynamic analysis===
1280
1281
The static formulation can be readily extended for the transient dynamic case accounting for geometrical and material non linear effects. Indeed in many situations of this kind,  material quasi-incompressibility develops in specific zones of the solid due to the accumulation of plastic strains. It is well known that in these cases the use of equal order interpolations for displacements and pressure leads to locking solutions unless some precautions are taken.
1282
1283
The transient FIC equilibrium equations  can be written in an identical form to eq. ([[#eq-65|65]]) (neglecting time stabilization terms <span id='citeF-32'></span><span id='citeF-37'></span><span id='citeF-38'></span>[[#cite-32|[32,37,38]]]) with
1284
1285
<span id="eq-81"></span>
1286
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1287
|-
1288
| 
1289
{| style="text-align: left; margin:auto;width: 100%;" 
1290
|-
1291
| style="text-align: center;" | <math>r_i:= - \rho \frac{\partial ^2{u_i}}{\partial{t}^2} +  {\partial \sigma _{ij} \over \partial x_j}+ b_i </math>
1292
|}
1293
| style="width: 5px;text-align: right;white-space: nowrap;" | (81)
1294
|}
1295
1296
where <math display="inline">\rho </math> is the  density.
1297
1298
Eq. ([[#eq-81|81]]) is completed with the constitutive equations for the deviatoric stresses (eq. ([[#eq-70|70]]))  and the  pressure (eq. ([[#eq-72|72]])), as well with the boundary conditions ([[#eq-67|67]]) and ([[#eq-68|68]]) and the initial conditions for <math display="inline">t=0</math>.
1299
1300
Following the arguments of the static case, the stabilized constitutive equation for the pressure can be expressed in terms of the residuals of the momentum equations by an expression identical to eq. ([[#eq-79|79]]). This  equation is now written in an incremental form more suitable for non linear transient analysis.
1301
1302
We note that the stabilization terms are not larger needed in the numerical (FEM) solution equilibrium equations which can be asted in the standard form. This simply implies neglecting the terms involving <math display="inline">h_i</math> in Eq. ([[#eq-74|74]]). The stabilized FIC form of the pressure constitutive equation (Eq. ([[#eq-79|79]])) is however mandatory in order to avoid volumetric locking.
1303
1304
The set of stabilized equations to be solved are now:
1305
1306
===Equilibrium===
1307
1308
<span id="eq-82"></span>
1309
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1310
|-
1311
| 
1312
{| style="text-align: left; margin:auto;width: 100%;" 
1313
|-
1314
| style="text-align: center;" | <math>r_i =0 </math>
1315
|}
1316
| style="width: 5px;text-align: right;white-space: nowrap;" | (82)
1317
|}
1318
1319
===Pressure~constitutive ~equation===
1320
1321
<span id="eq-83"></span>
1322
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1323
|-
1324
| 
1325
{| style="text-align: left; margin:auto;width: 100%;" 
1326
|-
1327
| style="text-align: center;" | <math>{\Delta p\over K} -{\partial (\Delta u_i) \over \partial x_i} - \sum \limits _{i=1}^{n_d} \tau _i {\partial r_i \over \partial x_i}=0 </math>
1328
|}
1329
| style="width: 5px;text-align: right;white-space: nowrap;" | (83)
1330
|}
1331
1332
where <math display="inline">\Delta p = p^{n+1}-p^n</math> and <math display="inline">\Delta u_i = u_i^{n+1}-u_i^n</math> are the increments of pressure and displacements, respectively. As usual <math display="inline">(\cdot )^n</math> denotes values at time <math display="inline">t_n</math>.
1333
1334
In the derivation of Eq. ([[#eq-83|83]]) we have accepted that <math display="inline">\Delta r_i =r_i^{n+1} \equiv r_i</math> as the infinitesimal equilibrium equations are assumed to be satisfied at time <math display="inline">t_n</math> (and hence <math display="inline">r_i^n =0</math>).
1335
1336
The weighted residual form of the FIC governing equations ([[#eq-82|82]]), ([[#eq-68|68]]) and ([[#eq-83|83]]) is
1337
1338
===Equilibrium===
1339
1340
<span id="eq-84"></span>
1341
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1342
|-
1343
| 
1344
{| style="text-align: left; margin:auto;width: 100%;" 
1345
|-
1346
| style="text-align: center;" | <math>\int _\Omega \!\delta u_i \left[-\rho \frac{\partial ^2{u_i}}{\partial{t}^2}+ {\partial \sigma _{ij} \over \partial x_j}+b_i\right] d\Omega + \!\!\int _{\Gamma _t}\! \delta u_i \left[\sigma _{ij}n_j - \bar t_i \right]d\Gamma =0 </math>
1347
|}
1348
| style="width: 5px;text-align: right;white-space: nowrap;" | (84)
1349
|}
1350
1351
===Pressure~constitutive ~equation===
1352
1353
<span id="eq-85"></span>
1354
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1355
|-
1356
| 
1357
{| style="text-align: left; margin:auto;width: 100%;" 
1358
|-
1359
| style="text-align: center;" | <math>\int _\Omega q \left({p\over K} - {\varepsilon }_{v}\right)d\Omega - \int _\Omega q \left(\sum \limits _{i=1}^{n_d} \tau _i {\partial {r}_{i} \over \partial x_i}\right)\, d\Omega =0 </math>
1360
|}
1361
| style="width: 5px;text-align: right;white-space: nowrap;" | (85)
1362
|}
1363
1364
where <math display="inline">\delta u_i </math> and <math display="inline">q</math> are arbitrary test functions representing virtual displacements and virtual pressure fields, respectively.
1365
1366
Integrating by parts, the terms involving <math display="inline">s_{ij}</math> and <math display="inline">p</math> in Eq. ([[#eq-84|84]]) and the term involving <math display="inline">r_i</math> in Eq. ([[#eq-85|85]]) and neglecting the space derivatives of the intrinsic time parameters leads to
1367
1368
<span id="eq-86"></span>
1369
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1370
|-
1371
| 
1372
{| style="text-align: left; margin:auto;width: 100%;" 
1373
|-
1374
| style="text-align: center;" | <math>\int _{\Omega }\delta u_i \rho \frac{\partial ^2{u_i}}{\partial{t}^2}d\Omega +  \int _{\Omega }\!\delta {\varepsilon }_{ij}\sigma _{ij}\,d\Omega{-} \int _{\Omega }\delta {u}_i{b}_i\, d\Omega - \int _{\Gamma _t}\delta {u}_i\bar{t}_i\,d\Omega =0 </math>
1375
|}
1376
| style="width: 5px;text-align: right;white-space: nowrap;" | (86)
1377
|}
1378
1379
<span id="eq-87"></span>
1380
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1381
|-
1382
| 
1383
{| style="text-align: left; margin:auto;width: 100%;" 
1384
|-
1385
| style="text-align: center;" | <math>\int _\Omega q \left({p\over K} - {\varepsilon }_{v}\right)d\Omega + \int _\Omega \left(\sum \limits _{i=1}^{n_d} {\partial q \over \partial x_i} \tau _{i} r_i \right)\,d\Omega =0 </math>
1386
|}
1387
| style="width: 5px;text-align: right;white-space: nowrap;" | (87)
1388
|}
1389
1390
The residual <math display="inline">r_i</math> is split now as
1391
1392
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1393
|-
1394
| 
1395
{| style="text-align: left; margin:auto;width: 100%;" 
1396
|-
1397
| style="text-align: center;" | <math>r_i = \pi _i + {\partial p \over \partial x_i} </math>
1398
|}
1399
| style="width: 5px;text-align: right;white-space: nowrap;" | (88)
1400
|}
1401
1402
where
1403
1404
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1405
|-
1406
| 
1407
{| style="text-align: left; margin:auto;width: 100%;" 
1408
|-
1409
| style="text-align: center;" | <math>\pi _i = - \frac{\partial ^2{u_i}}{\partial{t}^2} + {\partial s_{ij} \over \partial x_i}+b_i </math>
1410
|}
1411
| style="width: 5px;text-align: right;white-space: nowrap;" | (89)
1412
|}
1413
1414
Note that <math display="inline">\pi _i</math> is the part of <math display="inline">r_i</math> not containing the pressure gradient and may be interpreted as the negative of a projection of the pressure gradient. In a discrete setting the terms <math display="inline">\pi _i</math> can be considered belonging to a ''sub-scale'' space orthogonal to that of the pressure gradient terms.
1415
1416
In the infinitesimal limit <math display="inline">r_i=0</math> and <math display="inline">{\partial p \over \partial x_i} + \pi _i=0</math>. This limit relationship between <math display="inline">{\partial p \over \partial x_i}</math> and <math display="inline">\pi _i</math> can be weakly enforced by means of a weighted residual form.
1417
1418
The final set of integral equations is therefore
1419
1420
<span id="eq-90"></span>
1421
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1422
|-
1423
| 
1424
{| style="text-align: left; margin:auto;width: 100%;" 
1425
|-
1426
| style="text-align: center;" | <math>\int _{\Omega }\delta u_i \rho \frac{\partial ^2{u_i}}{\partial{t}^2}\,d\Omega + \int _{\Omega }\delta {\varepsilon }_{ij}{\sigma }_{ij}\,d\Omega{-} \int _{\Omega }\delta u_i b_i \,d\Omega  - \int _{\Gamma _t}\delta {u}_i\bar{t}_i \,d\Gamma _t=0 </math>
1427
|}
1428
| style="width: 5px;text-align: right;white-space: nowrap;" | (90)
1429
|}
1430
1431
<span id="eq-91"></span>
1432
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1433
|-
1434
| 
1435
{| style="text-align: left; margin:auto;width: 100%;" 
1436
|-
1437
| style="text-align: center;" | <math>\int _{\Omega }q \left({\Delta p\over K} - {\partial (\Delta u_i) \over \partial x_i} \right)\,d\Omega + \int _\Omega \left[ \sum \limits _{i=1}^{n_d}{\partial q \over \partial x_i} \tau _{i}\left({\partial p \over \partial x_i}+\pi _i \right)\right]d\Omega =0 </math>
1438
|}
1439
| style="width: 5px;text-align: right;white-space: nowrap;" | (91)
1440
|}
1441
1442
<span id="eq-92"></span>
1443
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1444
|-
1445
| 
1446
{| style="text-align: left; margin:auto;width: 100%;" 
1447
|-
1448
| style="text-align: center;" | <math>\int _{\Omega }\left[\sum \limits _{i=1}^{n_d} w_i \tau _{i} \left({\partial p \over \partial x_i}+\pi _i \right)\right]d\Omega =0 </math>
1449
|}
1450
| style="width: 5px;text-align: right;white-space: nowrap;" | (92)
1451
|}
1452
1453
where the <math display="inline">\tau _i</math> coefficients are introduced in Eq. ([[#eq-92|92]]) for convenience.
1454
1455
We will choose <math display="inline">C^0</math>continuous linear interpolations of the displacements, the pressure and the pressure gradient projection <math display="inline">\pi _i</math> over three-node triangles (2D) and four-node tetrahedra (3D) <span id='citeF-3'></span>[[#cite-3|[3]]]. The linear interpolations are  written as
1456
1457
<span id="eq-93"></span>
1458
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1459
|-
1460
| 
1461
{| style="text-align: left; margin:auto;width: 100%;" 
1462
|-
1463
| style="text-align: center;" | <math>{u}_i=\sum _{j=1}^{n} N_j \bar{u}_i^j \quad , \quad p=\sum _{j=1}^{n} N_j \bar{p}^j \quad , \quad  \pi _i =\sum _{j=1}^{n} N_j \bar \pi ^j_i </math>
1464
|}
1465
| style="width: 5px;text-align: right;white-space: nowrap;" | (93)
1466
|}
1467
1468
where <math display="inline">n=3(4)</math> for 2D(3D) problems and <math display="inline">(\bar{\cdot })</math> denotes nodal variables. As usual <math display="inline">N_j</math> are the linear shape functions <span id='citeF-3'></span>[[#cite-3|[3]]]. The nodal variables are a function of the time <math display="inline">t</math>. Substituting the approximations ([[#eq-93|93]])  into eqs. ([[#eq-90|90]])&#8211;([[#eq-92|92]]) gives the following system of discretized equations
1469
1470
<span id="eq-94"></span>
1471
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1472
|-
1473
| 
1474
{| style="text-align: left; margin:auto;width: 100%;" 
1475
|-
1476
| style="text-align: center;" | <math>{M} \ddot{\bar {u}} + {g} -  {f}={0} </math>
1477
|}
1478
| style="width: 5px;text-align: right;white-space: nowrap;" | (94)
1479
|}
1480
1481
<span id="eq-95"></span>
1482
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1483
|-
1484
| 
1485
{| style="text-align: left; margin:auto;width: 100%;" 
1486
|-
1487
| style="text-align: center;" | <math>{G}^T \Delta \bar {u} -{C} \Delta \bar {p}- {L} \bar {p} - {Q}\bar {\boldsymbol \pi }= {0}  </math>
1488
|}
1489
| style="width: 5px;text-align: right;white-space: nowrap;" | (95)
1490
|}
1491
1492
<span id="eq-96"></span>
1493
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1494
|-
1495
| 
1496
{| style="text-align: left; margin:auto;width: 100%;" 
1497
|-
1498
| style="text-align: center;" | <math>{Q}^T \bar {p}+ \bar {C} \bar {\boldsymbol \pi }=  {0}  </math>
1499
|}
1500
| style="width: 5px;text-align: right;white-space: nowrap;" | (96)
1501
|}
1502
1503
where <math display="inline">\ddot{\bar {u}}</math> is the nodal acceleration vector. The expression of the different matrices in vectors appearing in Eqs. ([[#eq-94|94]])&#8211;([[#eq-96|96]]) is presented in the Appendix.
1504
1505
A four steps semi-implicit time integration algorithm can be derived from eqs. ([[#eq-94|94]])&#8211;([[#eq-96|96]]) as follows
1506
1507
step
1508
1509
'''Step step.'''step
1510
* Compute the nodal velocities <math display="inline">\dot{\bar {u}}^{n+1/2}</math>
1511
1512
<span id="eq-97"></span>
1513
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1514
|-
1515
| 
1516
{| style="text-align: left; margin:auto;width: 100%;" 
1517
|-
1518
| style="text-align: center;" | <math>
1519
1520
\dot{\bar {u}}^{n+1/2} = \dot{\bar {u}}^{n-1/2}+\Delta t {M}^{-1}_d ({f}^n - {g}^n)  </math>
1521
|}
1522
| style="width: 5px;text-align: right;white-space: nowrap;" | (97)
1523
|}
1524
1525
* Compute the nodal displacements <math display="inline">\bar {u}^{n+1}</math>
1526
1527
<span id="eq-98"></span>
1528
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1529
|-
1530
| 
1531
{| style="text-align: left; margin:auto;width: 100%;" 
1532
|-
1533
| style="text-align: center;" | <math>
1534
1535
\bar {u}^{n+1} = \bar{u}^n + \Delta t \dot{\bar {u}}^{n+1/2}  </math>
1536
|}
1537
| style="width: 5px;text-align: right;white-space: nowrap;" | (98)
1538
|}
1539
* Compute the nodal pressures <math display="inline">\bar {p}^{n+1}</math>
1540
1541
<span id="eq-99"></span>
1542
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1543
|-
1544
| 
1545
{| style="text-align: left; margin:auto;width: 100%;" 
1546
|-
1547
| style="text-align: center;" | <math>
1548
1549
\bar {p}^{n+1}=[ {C} + {L}]^{-1} [\Delta t {G}^T \dot{\bar {u}}^{n+1/2}+  {C} \bar {p}^n- {Q} \bar {\boldsymbol \pi }^n]  </math>
1550
|}
1551
| style="width: 5px;text-align: right;white-space: nowrap;" | (99)
1552
|}
1553
* Compute the nodal projected pressure gradients <math display="inline">\bar {\boldsymbol \pi }^{n+1}</math>
1554
1555
<span id="eq-100"></span>
1556
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1557
|-
1558
| 
1559
{| style="text-align: left; margin:auto;width: 100%;" 
1560
|-
1561
| style="text-align: center;" | <math>
1562
1563
\bar {\boldsymbol \pi }^{n+1}=- \bar {C}_d^{-1} {Q}^T \bar {p}^{n+1}  </math>
1564
|}
1565
| style="width: 5px;text-align: right;white-space: nowrap;" | (100)
1566
|}
1567
1568
In above, all matrices are evaluated at <math display="inline">t^{n+1}</math>, <math display="inline">(\cdot )_d =\hbox{diag}\, (\cdot )</math> and
1569
1570
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1571
|-
1572
| 
1573
{| style="text-align: left; margin:auto;width: 100%;" 
1574
|-
1575
| style="text-align: center;" | <math>{g}^n =\int _{\Omega ^e} [{B}^T {\boldsymbol \sigma }]^n\,d\Omega  </math>
1576
|}
1577
| style="width: 5px;text-align: right;white-space: nowrap;" | (101)
1578
|}
1579
1580
where the stresses <math display="inline">{\boldsymbol \sigma }^n</math> are obtained by consistent integration of the adequate (non linear)  constitutive law [33].
1581
1582
Note that steps 1, 2 and 4 are fully explicit as a diagonal form of matrices <math display="inline">{C}</math> and <math display="inline">\bar{C}</math> has been chosen. The solution of step 3 with  a diagonal form for <math display="inline">{C}</math> still requires the inverse of a Laplacian matrix. This can be an inexpensive process using an iterative equation solution method (e.g. a preconditioned conjugate gradient method).
1583
1584
A ''three steps'' approach can be obtained by evaluating the projected pressure gradient variables <math display="inline">\bar {\boldsymbol \pi }^{n+1}</math>  at <math display="inline">t_{n+1}</math> in a fully implicit form in Eq. ([[#eq-99|99]]). Eliminating then <math display="inline">\bar {\boldsymbol \pi }^{n+1}</math>  from the fourth step using Eq. ([[#eq-100|100]]) and substituting this expression into Eq. ([[#eq-99|99]]) leads to
1585
1586
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1587
|-
1588
| 
1589
{| style="text-align: left; margin:auto;width: 100%;" 
1590
|-
1591
| style="text-align: center;" | <math>\bar{p}^{n+1}= [{C} + {L}-{S}]^{-1} [\Delta t {G}^T \dot{\bar {u}}^{n+1/2} + {C} \bar{p}^n ] </math>
1592
|}
1593
| style="width: 5px;text-align: right;white-space: nowrap;" | (102)
1594
|}
1595
1596
where
1597
1598
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1599
|-
1600
| 
1601
{| style="text-align: left; margin:auto;width: 100%;" 
1602
|-
1603
| style="text-align: center;" | <math>{S} = {Q} \bar {C}_d^{-1} {Q}^T </math>
1604
|}
1605
| style="width: 5px;text-align: right;white-space: nowrap;" | (103)
1606
|}
1607
1608
Recall that for the full incompressible case <math display="inline">K=\infty </math> and <math display="inline">{C}=0 </math> in all above equations.
1609
1610
The critical time step <math display="inline">\Delta t</math> is taken as that of the standard explicit dynamic scheme [28].
1611
1612
===Explicit algorithm===
1613
1614
A fully explicit algorithm can be obtained by computing <math display="inline">\bar {p}^{n+1}</math> from step 3 in Eq. ([[#eq-99|99]])  as follows
1615
1616
<span id="eq-104"></span>
1617
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1618
|-
1619
| 
1620
{| style="text-align: left; margin:auto;width: 100%;" 
1621
|-
1622
| style="text-align: center;" | <math>\bar{p}^{n+1}={C}_d^{-1}[\Delta t {G}^T \dot{\bar {u}}^{n+1/2} + ({C}_d+{L}) \bar{p}^n - {Q} {\boldsymbol \sigma }^n] </math>
1623
|}
1624
| style="width: 5px;text-align: right;white-space: nowrap;" | (104)
1625
|}
1626
1627
Obviously, solution of eq. ([[#eq-104|104]]) breaks down for <math display="inline">K=\infty </math> as <math display="inline">{C}=0 </math> in this case. Therefore, the explicit algorithm is not applicable in the full incompressible limit. The explicit form can however be used with success in problems where quasi-incompressible regions exist adjacent to standard “compressible” zones.
1628
1629
==5 DEM&#8211;FEM CONTACT ALGORITHM==
1630
1631
Use of combined DEM/FEM models involves treatment of contact between spherical discrete elements and boundary of continuum subdomains discretized with finite elements as shown in Fig.  [[#img-13|13]]. Similarly as in case of contact between two spheres (Fig. [[#img-2|2]]) the contact force between the sphere and  external edge <math display="inline">{F}</math>  of a finite element is decomposed into normal and tangential components, <math display="inline">{F}_n</math>and <math display="inline">{ F}_T</math>. In a general case the contact model between discrete elements and finite element edges can include cohesion, damping, friction, wear, heat generation and exchange.
1632
1633
<div id='img-13'></div>
1634
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1635
|-
1636
|[[Image:Draft_Samper_102318540-fesphere.png|270px|Contact between a sphere and finite element edge]]
1637
|- style="text-align: center; font-size: 75%;"
1638
| colspan="1" | '''Figure 13:''' Contact between a sphere and finite element edge
1639
|}
1640
1641
Similarly to the relationship ([[#eq-10|10]])  the normal contact force <math display="inline">F_n</math> is decomposed into the elastic part <math display="inline">F_{ne}</math> and the damping contact force <math display="inline">F_{nd}</math>
1642
1643
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1644
|-
1645
| 
1646
{| style="text-align: left; margin:auto;width: 100%;" 
1647
|-
1648
| style="text-align: center;" | <math>F_n = F_{ne} + F_{nd}\,. </math>
1649
|}
1650
| style="width: 5px;text-align: right;white-space: nowrap;" | (105)
1651
|}
1652
1653
The elastic part of the normal contact force <math display="inline">F_{ne}</math> is proportional to the normal stiffness <math display="inline">k_n</math> and to the penetration of the two particle surfaces <math display="inline">u_{rn}</math>
1654
1655
<span id="eq-106"></span>
1656
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1657
|-
1658
| 
1659
{| style="text-align: left; margin:auto;width: 100%;" 
1660
|-
1661
| style="text-align: center;" | <math>F_{ne}=k_n u_{rn}\,. </math>
1662
|}
1663
| style="width: 5px;text-align: right;white-space: nowrap;" | (106)
1664
|}
1665
1666
The penetration is calculated here as
1667
1668
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1669
|-
1670
| 
1671
{| style="text-align: left; margin:auto;width: 100%;" 
1672
|-
1673
| style="text-align: center;" | <math>u_{rn} = d - r\,, </math>
1674
|}
1675
| style="width: 5px;text-align: right;white-space: nowrap;" | (107)
1676
|}
1677
1678
where <math display="inline">d</math> is the distance of the particle centre <math display="inline">C</math> and <math display="inline">r</math> is  the particle radius.
1679
1680
The contact damping force is assumed to be of viscous type
1681
1682
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1683
|-
1684
| 
1685
{| style="text-align: left; margin:auto;width: 100%;" 
1686
|-
1687
| style="text-align: center;" | <math>F_{nd}=c_n v_{rn} </math>
1688
|}
1689
| style="width: 5px;text-align: right;white-space: nowrap;" | (108)
1690
|}
1691
1692
where <math display="inline">c_n</math> is the damping coefficient and <math display="inline">v_{rn}</math> is the normal component of the relative velocity at the contact point. The relative velocity at the contact point is calculated as
1693
1694
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1695
|-
1696
| 
1697
{| style="text-align: left; margin:auto;width: 100%;" 
1698
|-
1699
| style="text-align: center;" | <math>{v}_{r}= ({{v}}_{C} +{\omega }_{C} \times {r})- ({{v}}_{1}H_{1} + {{v}}_{2}H_{2}) \,, </math>
1700
|}
1701
| style="width: 5px;text-align: right;white-space: nowrap;" | (109)
1702
|}
1703
1704
where <math display="inline">{{v}}_{C} +{\omega }_{C} \times {r}</math> is the discrete element velocity at the contact point and <math display="inline">{{v}}_{1}H_{1} + {{v}}_{2}H_{2}</math> is the velocity of the finite element edge at the contact point expressed in terms of the nodal velocities, <math display="inline">{{v}}_{1}</math> and <math display="inline">{{v}}_{2}</math>, and respective shape functions <math display="inline">H_{1}</math> and <math display="inline">H_{2}</math>. The relative velocity is decomposed into the normal and tangential velocity using the formulae ([[#eq-15|15]]) and ([[#eq-17|17]]).
1705
1706
The cohesive contact force is calculated using the  elastic-perfectly brittle model described in Section [[#3.1 Elastic perfectly brittle model|3.1]]. In the absence of cohesion the frictional contact reaction is calculated using the regularized Coulomb friction model described in section [[#lb-2.3.3|2.3.3]].
1707
1708
==6 GENERAL DISCRETE ELEMENT-FINITE ELEMENT DYNAMIC FORMULATION==
1709
1710
The general algorithm for the transient dynamic problem involving discrete elements and (stabilized) finite elements includes the following steps.
1711
1712
===Step1. Compute the nodal velocities===
1713
1714
====Discrete elements====
1715
1716
<span id="eq-110"></span>
1717
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1718
|-
1719
| 
1720
{| style="text-align: left; margin:auto;width: 100%;" 
1721
|-
1722
| style="text-align: center;" | <math>\dot{{u}}_i^{n+1/2}=\dot{{u}}_i^{n-1/2}+\ddot{{u}}_i^{n}{\Delta t}\,, </math>
1723
|}
1724
| style="width: 5px;text-align: right;white-space: nowrap;" | (110)
1725
|}
1726
1727
with
1728
1729
<span id="eq-111"></span>
1730
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1731
|-
1732
| 
1733
{| style="text-align: left; margin:auto;width: 100%;" 
1734
|-
1735
| style="text-align: center;" | <math>\ddot{{u}}_i^{n}=\frac{{{F}}_i^n}{m_i}\,, </math>
1736
|}
1737
| style="width: 5px;text-align: right;white-space: nowrap;" | (111)
1738
|}
1739
1740
<span id="eq-112"></span>
1741
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1742
|-
1743
| 
1744
{| style="text-align: left; margin:auto;width: 100%;" 
1745
|-
1746
| style="text-align: center;" | <math>{{\omega }}_i^{n+1/2}={{\omega }}_i^{n-1/2}+\dot{{\omega }}_i^{n}{\Delta t}\,. </math>
1747
|}
1748
| style="width: 5px;text-align: right;white-space: nowrap;" | (112)
1749
|}
1750
1751
with
1752
1753
<span id="eq-113"></span>
1754
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1755
|-
1756
| 
1757
{| style="text-align: left; margin:auto;width: 100%;" 
1758
|-
1759
| style="text-align: center;" | <math>\dot{{\omega }}_i^{n}=\frac{{{T}}_i^n}{I_i}\,, </math>
1760
|}
1761
| style="width: 5px;text-align: right;white-space: nowrap;" | (113)
1762
|}
1763
1764
====Finite elements====
1765
1766
<span id="eq-114"></span>
1767
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1768
|-
1769
| 
1770
{| style="text-align: left; margin:auto;width: 100%;" 
1771
|-
1772
| style="text-align: center;" | <math>\dot{\bar {u}}^{n+1/2} = \dot{\bar {u}}^{n-1/2}+\Delta t {M}^{-1}_d ({f}^n - {g}^n)  </math>
1773
|}
1774
| style="width: 5px;text-align: right;white-space: nowrap;" | (114)
1775
|}
1776
1777
===Step2. Compute the nodal displacements===
1778
1779
====Discrete elements====
1780
1781
<span id="eq-115"></span>
1782
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1783
|-
1784
| 
1785
{| style="text-align: left; margin:auto;width: 100%;" 
1786
|-
1787
| style="text-align: center;" | <math>{{u}}_i^{n+1}={{u}}_i^{n}+\dot{{u}}_i^{n+1/2}{\Delta t}\,. </math>
1788
|}
1789
| style="width: 5px;text-align: right;white-space: nowrap;" | (115)
1790
|}
1791
1792
<span id="eq-116"></span>
1793
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1794
|-
1795
| 
1796
{| style="text-align: left; margin:auto;width: 100%;" 
1797
|-
1798
| style="text-align: center;" | <math>{{\Delta \theta }}_i={{\omega }}_i^{n+1/2}{\Delta t}\,, </math>
1799
|}
1800
| style="width: 5px;text-align: right;white-space: nowrap;" | (116)
1801
|}
1802
1803
====Finite elements====
1804
1805
<span id="eq-117"></span>
1806
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1807
|-
1808
| 
1809
{| style="text-align: left; margin:auto;width: 100%;" 
1810
|-
1811
| style="text-align: center;" | <math>\bar {u}^{n+1} = \bar{u}^n + \Delta t \dot{\bar {u}}^{n+1/2}  </math>
1812
|}
1813
| style="width: 5px;text-align: right;white-space: nowrap;" | (117)
1814
|}
1815
1816
===Step3. Compute the nodal pressures===
1817
1818
====Finite elements====
1819
1820
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1821
|-
1822
| 
1823
{| style="text-align: left; margin:auto;width: 100%;" 
1824
|-
1825
| style="text-align: center;" | <math>\bar{p}^{n+1}= [{C} + {L}-{S}]^{-1} [\Delta t {G}^T \dot{\bar {u}}^{n+1/2} + {C} \bar{p}^n ] </math>
1826
|}
1827
| style="width: 5px;text-align: right;white-space: nowrap;" | (118)
1828
|}
1829
1830
===Step4. Update the nodal coordinates===
1831
1832
====Discrete and finite elements====
1833
1834
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1835
|-
1836
| 
1837
{| style="text-align: left; margin:auto;width: 100%;" 
1838
|-
1839
| style="text-align: center;" | <math>{X}_i^{n+1} = {X}_i^{n}+\Delta \bar {u}_i </math>
1840
|}
1841
| style="width: 5px;text-align: right;white-space: nowrap;" | (119)
1842
|}
1843
1844
{| class="formulaSCP" style="width: 100%; text-align: left;" 
1845
|-
1846
| 
1847
{| style="text-align: left; margin:auto;width: 100%;" 
1848
|-
1849
| style="text-align: center;" | <math>\Delta \bar {u}_i =\bar {u}_i^{n+1} - \bar {u}_i^n </math>
1850
|}
1851
| style="width: 5px;text-align: right;white-space: nowrap;" | (120)
1852
|}
1853
1854
===Step 5. Check frictional contact forces===
1855
1856
===Step 6. Update residual force vector  Go to step 1===
1857
1858
==7 EXAMPLES==
1859
1860
===7.1 Uniaxial compression and tension tests of a rock material===
1861
1862
<div id='img-14'></div>
1863
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1864
|-
1865
|[[Image:Draft_Samper_102318540-ucs-lab.png|462px|Rock specimen after laboratory unconfined compression test]]
1866
|- style="text-align: center; font-size: 75%;"
1867
| colspan="1" | '''Figure 14:''' Rock specimen after laboratory unconfined compression test
1868
|}
1869
1870
Mechanical properties of rock materials are determined from laboratory tests such as compression, triaxial and tension tests. Figure [[#img-14|14]] shows a cube specimen of a sandstone rock after unconfined compression test. The parameters of numerical discrete element model can be obtained carrying out simulations of basic laboratory tests. Such analyses allow us to determine microscopic constitutive parameters for a material sample modelled with discrete elements yielding required macroscopic properties. The most important macroscopic properties are Young modulus, Poisson ratio and compressive and tensile strengths. Microscopic parameters are in turn all the   constitutive model parameters governing the frictional contact interaction between a pair of particles. Here elastic perfectly brittle constitutive model will be used with the contact stiffness in normal and tangential direction, interface tensile and shear strength and Coulomb friction coefficient being the microscopic parameters.
1871
1872
<div id='img-15'></div>
1873
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1874
|-
1875
|[[Image:Draft_Samper_102318540-uni-comp-geo.png|360px|Set-up  of the numerical model of uncofined compression test]]
1876
|- style="text-align: center; font-size: 75%;"
1877
| colspan="1" | '''Figure 15:''' Set-up  of the numerical model of uncofined compression test
1878
|}
1879
1880
Figure [[#img-15|15]] presents the material sample prepared for testing. A material sample of <math display="inline">109\times 109</math> mm is represented by an assembly of randomly compacted 2100 discs of radii 1&#8211;1.5 mm. It is shown in <span id='citeF-40'></span>[[#cite-40|[40]]] that preparing a well-connected densely packed irregular assembly of particles is the key to successful simulation with discrete elements. Compaction of the particle assembly shown in  Fig. [[#img-15|15]] is characterized by a porosity of <math display="inline">13</math>%
1881
1882
The macroscopic response of a square material sample subjected to uniaxial compression has been studied. The loading has been introduced under kinematic control by prescribing the motion of rigid plates pressing on the top and bottom  of the sample. The deformation in the <math display="inline">x</math> direction was free. The velocity of the wall displacement   was 1 mm/s which was  found to be sufficiently low to obtain quasi-static loading.
1883
1884
The purpose of the present study was to determine the microscopic parameters for sandstone with the following macroscopic properties: Young modulus <math display="inline">E=14</math> GPa, Poisson ratio <math display="inline">\nu=0.3</math>  and unconfined compressive strength <math display="inline">\sigma _{c}=60</math> MPa. The stress-strain relationship shown in Fig. [[#img-16|16]] satisfying these requirements has been obtained with the following set of micromechanical parameters: contact stiffness in the normal and tangential directions <math display="inline">k_n=k_T=20</math>  GPa, Coulomb friction coefficient <math display="inline">\mu=0.839</math> and cohesive bond strengths in the normal and tangential direction, <math display="inline">R_n=0.1</math> MN/m and <math display="inline">R_T=1</math> MN/m, respectively. The failure mode of the specimen obtained in the simulation is shown in Fig. [[#img-17|17]]. Particles with broken cohesive bonds are marked with different colour. Comparison of Figs. [[#img-14|14]] and [[#img-17|17]] show that numerical analysis yields a failure mode  similar to that observed in experiments.
1885
1886
<div id='img-16'></div>
1887
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1888
|-
1889
|[[Image:Draft_Samper_102318540-com-brittle.png|450px|Stress&#8211;strain relationship for uniaxial unconfined compression test]]
1890
|- style="text-align: center; font-size: 75%;"
1891
| colspan="1" | '''Figure 16:''' Stress&#8211;strain relationship for uniaxial unconfined compression test
1892
|}
1893
1894
<div id='img-17'></div>
1895
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1896
|-
1897
|[[Image:Draft_Samper_102318540-com-fail.png|390px|Failure mode obtained in the simulation of unconfined compression test]]
1898
|- style="text-align: center; font-size: 75%;"
1899
| colspan="1" | '''Figure 17:''' Failure mode obtained in the simulation of unconfined compression test
1900
|}
1901
1902
Tensile strength of rocks is obtained experimentally from indirect tension (Brasilian) test. In the numerical procedure we checked the tensile strength of the particle assembly carrying out simulations of direct tension test. The stress&#8211;strain relationship obtained in this simulation is plotted in Fig. [[#img-18|18]].  The elastic modulus for tension and tensile strength obtained from the plots are the following: <math display="inline">E_t=11.5</math> GPa, <math display="inline">\sigma _{t}=12.7</math> MPa. It can be noted that the curves both for compression and tension are typical for materials with brittle failure. The failure mode obtained in the analysis of a direct tension test is shown in Fig. [[#img-19|19]].
1903
1904
<div id='img-18'></div>
1905
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1906
|-
1907
|[[Image:Draft_Samper_102318540-ten-brittle.png|450px|Stress&#8211;strain relationship for uniaxial direct tension test]]
1908
|- style="text-align: center; font-size: 75%;"
1909
| colspan="1" | '''Figure 18:''' Stress&#8211;strain relationship for uniaxial direct tension test
1910
|}
1911
1912
<div id='img-19'></div>
1913
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1914
|-
1915
|[[Image:Draft_Samper_102318540-ten-bb.png|354px|Failure mode obtained in the simulation of direct tension test]]
1916
|- style="text-align: center; font-size: 75%;"
1917
| colspan="1" | '''Figure 19:''' Failure mode obtained in the simulation of direct tension test
1918
|}
1919
1920
===7.2 Simulation of rock cutting===
1921
1922
<div id='img-20'></div>
1923
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1924
|-
1925
|[[Image:Draft_Samper_102318540-cut-fem12-ini.png|600px|Rock cutting problem. 2100 discrete elements are used to model the rock whereas 6800 linear triangles are used to discretize the tool.]]
1926
|- style="text-align: center; font-size: 75%;"
1927
| colspan="1" | '''Figure 20:''' Rock cutting problem. 2100 discrete elements are used to model the rock whereas 6800 linear triangles are used to discretize the tool.
1928
|}
1929
1930
A material sample representing sandstone with parameters determined in the previous section has been used in the simulation of a rock cutting process. Two different solutions with the tool discretized with discrete and finite elements have been studied. The combined FEM/DEM model is shown in Fig. [[#img-20|20]]. The material sample <math display="inline">109\times 109</math> mm is represented by an assembly of randomly compacted 2100 discs of radii 1&#8211;1.5 mm. Contact interface between particles representing the rock is characterized with the stiffness in the normal and tangential directions <math display="inline">k_n=k_T=20</math>  GPa, the cohesive bond strengths <math display="inline">R_n=0.1</math> MN/m, <math display="inline">R_T=1</math> MN/m, and the friction coefficient <math display="inline">\mu=0.839</math>. The  tool has been modelled with 6800 linear triangles. The following elasto-plastic parameters have been assumed for the tool material: Young modulus <math display="inline">E=2\cdot 10^5</math> MPa, Poisson's ratio <math display="inline">\nu=0.3</math>, yield stress <math display="inline">\sigma _Y=600</math> MPa and hardening modulus <math display="inline">300</math> MPa. The following parameters have been assumed for the tool-rock interface: contact stiffness modulus <math display="inline">k_n=20</math> GPa, Coulomb friction coefficient <math display="inline">\mu = 0.839</math>.
1931
1932
Cutting process has been carried out with prescribed horizontal velocity of the tool of 4 m/s. The cutting is shown in Fig. [[#img-21|21]]. In this figure we can also see the failure mode. Particles with broken bonds are coloured in blue. It can be clearly seen the formation of a chip during cutting which is typical for brittle materials.
1933
1934
Variation of the cutting force in the initial phase of cutting is shown in Fig. [[#img-22|22]]. The Mises stress and efffective plastic strain distributions in the tool in the initial phase of cutting are shown in Figs. [[#img-23|23]] and [[#img-24|24]], respectively.
1935
1936
<div id='img-21'></div>
1937
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1938
|-
1939
|[[Image:Draft_Samper_102318540-cut-fem12-bb-0015.png|468px|]]
1940
|[[Image:Draft_Samper_102318540-cut-fem12-bb-0025.png|450px|]]
1941
|-
1942
|[[Image:Draft_Samper_102318540-cut-fem12-bb-0065.png|426px|]]
1943
|[[Image:Draft_Samper_102318540-cut-fem12-bb-012.png|378px|Rock cutting process analysed using the FEM/DEM model. The failure mode]]
1944
|- style="text-align: center; font-size: 75%;"
1945
| colspan="2" | '''Figure 21:''' Rock cutting process analysed using the FEM/DEM model. The failure mode
1946
|}
1947
1948
<div id='img-22'></div>
1949
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1950
|-
1951
|[[Image:Draft_Samper_102318540-cut-fem-tcx.png|450px|Variation of the cutting force in the initial phase of cutting]]
1952
|- style="text-align: center; font-size: 75%;"
1953
| colspan="1" | '''Figure 22:''' Variation of the cutting force in the initial phase of cutting
1954
|}
1955
1956
<div id='img-23'></div>
1957
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1958
|-
1959
|[[Image:Draft_Samper_102318540-cut-fem12ns-vms-0009.png|582px|Rock cutting process analysed using the FEM/DEM model. The Mises stresses in the tool in the initial phase of cutting]]
1960
|- style="text-align: center; font-size: 75%;"
1961
| colspan="1" | '''Figure 23:''' Rock cutting process analysed using the FEM/DEM model. The Mises stresses in the tool in the initial phase of cutting
1962
|}
1963
1964
<div id='img-24'></div>
1965
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1966
|-
1967
|[[Image:Draft_Samper_102318540-cut-fem12ns-eps-0009.png|582px|Rock cutting process analysed using the FEM/DEM model. The effective plastic strain in the tool in the initial phase of cutting]]
1968
|- style="text-align: center; font-size: 75%;"
1969
| colspan="1" | '''Figure 24:''' Rock cutting process analysed using the FEM/DEM model. The effective plastic strain in the tool in the initial phase of cutting
1970
|}
1971
1972
The model of rock cutting process presented above has been used to demonstrate the capabilities of the developed algorithm for wear estimation. Next the tool has been discretized with 4500 discrete elements (Fig. [[#img-25|25]]). This allows us to easily modify the shape of the tool. The tool shape is changed by eliminating particles if the accumulated wear exceeds the particle diameter. Particles are eliminated automatically and the analysis is continued with the modified shape of the tool.
1973
1974
<div id='img-25'></div>
1975
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1976
|-
1977
|[[Image:Draft_Samper_102318540-w2d-ini.png|408px|Rock cutting process analyzed with the tool discretized with 4500 discrete elements]]
1978
|- style="text-align: center; font-size: 75%;"
1979
| colspan="1" | '''Figure 25:''' Rock cutting process analyzed with the tool discretized with 4500 discrete elements
1980
|}
1981
1982
The required material constants <math display="inline">k</math> and <math display="inline">H</math> necessary for  wear evaluation have been obtained in  laboratory tests. The average Brinell hardness of 400 has been obtained for the tool steel and a value of the Archard wear constant <math display="inline">k=0.0025</math> has been taken in the wear simulation.
1983
1984
Evolution of the tool shape is obtained in a number of consecutive simulations of work cycles. The tool shape is modified continuously during the analysis. After finishing one work cycle the tool with a new shape is again used in the simulation of next work cycle. The whole procedure can be repeated as needed to reproduce the real working conditions of the tool. In this manner the evolution of change of tool shape in wear process during its whole service life can be  obtained.
1985
1986
Tool shapes at different stages of wear are shown in Fig. [[#img-26|26]]. As much as 40% of the tool mass has been lost due to abrasive wear. Figure [[#img-27|27]] presents the tool mass loss as a fraction of the initial tool mass versus analysis time. The wear curve obtained in the analysis  agrees quite well with experimental curves obtained in field tests.
1987
1988
<div id='img-26'></div>
1989
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
1990
|-
1991
|[[Image:Draft_Samper_102318540-tshape-a.png|264px|]]
1992
|[[Image:Draft_Samper_102318540-tshape-b.png|390px|]]
1993
|-
1994
|[[Image:Draft_Samper_102318540-tshape-c.png|378px|]]
1995
|[[Image:Draft_Samper_102318540-tshape-d.png|378px|Evolution of the cutting tool shape due to wear during the rock cutting process]]
1996
|- style="text-align: center; font-size: 75%;"
1997
| colspan="2" | '''Figure 26:''' Evolution of the cutting tool shape due to wear during the rock cutting process
1998
|}
1999
2000
<div id='img-27'></div>
2001
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
2002
|-
2003
|[[Image:Draft_Samper_102318540-tmass-exp.png|600px|Mass loss fraction in the tool versus analysis time &#8211;- comparison of experimental and simulation wear curves]]
2004
|- style="text-align: center; font-size: 75%;"
2005
| colspan="1" | '''Figure 27:''' Mass loss fraction in the tool versus analysis time &#8211;- comparison of experimental and simulation wear curves
2006
|}
2007
2008
===7.3 Strip punch test===
2009
2010
The strip punch test investigated experimentally by Dede <span id='citeF-41'></span>[[#cite-41|[41]]] and studied numerically using a FE model by Klerck <span id='citeF-42'></span>[[#cite-42|[42]]] is analysed here using the DEM model. <div id='img-28'></div>
2011
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
2012
|-
2013
|[[Image:Draft_Samper_102318540-dede.png|528px|Strip punch geometry and loading configuration]]
2014
|- style="text-align: center; font-size: 75%;"
2015
| colspan="1" | '''Figure 28:''' Strip punch geometry and loading configuration
2016
|}
2017
Prismatic quartzite specimens of dimensions <math display="inline">80\times 80\times 80</math> mm are crushed by a punch load applied in the strip 10 mm wide on the top side along the symmetry line as shown in Fig. [[#img-28|28]]. A confining pressure of  4.5 MPa was applied to the specimen sides.
2018
2019
The experimental fracture pattern of the specimen is shown in Fig. [[#img-29|29]] and the experimental axial stress versus the specimen axial strain is plotted in Fig. [[#img-30|30]]. <div id='img-29'></div>
2020
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
2021
|-
2022
|[[Image:Draft_Samper_102318540-strip-exp.png|450px|Experimental fracture pattern, after <span id='citeF-41'></span>[[#cite-41|[41]]]]]
2023
|- style="text-align: center; font-size: 75%;"
2024
| colspan="1" | '''Figure 29:''' Experimental fracture pattern, after <span id='citeF-41'></span>[[#cite-41|[41]]]
2025
|}
2026
<div id='img-30'></div>
2027
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
2028
|-
2029
|[[Image:Draft_Samper_102318540-strip-forcex.png|600px|Experimental axial stress vs. axial strain plot, after <span id='citeF-41'></span>[[#cite-41|[41]]]]]
2030
|- style="text-align: center; font-size: 75%;"
2031
| colspan="1" | '''Figure 30:''' Experimental axial stress vs. axial strain plot, after <span id='citeF-41'></span>[[#cite-41|[41]]]
2032
|}
2033
2034
<div id='img-31'></div>
2035
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
2036
|-
2037
|[[Image:Draft_Samper_102318540-strip-ini.png|432px|2D discrete element model]]
2038
|- style="text-align: center; font-size: 75%;"
2039
| colspan="1" | '''Figure 31:''' 2D discrete element model
2040
|}
2041
2042
The 2D discrete element model used in our analysis is shown in Fig. [[#img-31|31]]. The square rock specimen <math display="inline">80\times 80</math> mm is represented by an assembly of randomly compacted 2100 discs of radii 0.8&#8211;1.2 mm. Sliding conditions are prescribed at the bottom specimen side and confining pressure is introduced by means of rigid plates transmitting a confining load. Axial load is introduced by a rigid punch under a prescribed displacement.
2043
2044
Required  microscopic parameters for the DEM model corresponding to the macroscopic properties of quartzite rock given in <span id='citeF-41'></span>[[#cite-41|[41]]] are: Young's modulus <math display="inline">E=68</math> GPA, compressive strength <math display="inline">200</math> MPa and tensile strength <math display="inline">27</math> MPa. Those values have been obtained using the procedure presented in example [[#7.1 Uniaxial compression and tension tests of a rock material|7.1]] for the following microscopic parameters: the stiffness in the normal and tangential directions <math display="inline">k_n=k_T=100</math>  GPa, the cohesive bond strengths <math display="inline">R_n=0.25</math> MN/m, <math display="inline">R_T=2.5</math> MN/m, and the friction coefficient <math display="inline">\mu=0.839</math>.
2045
2046
The numerical failure pattern is shown in Fig. [[#img-32|32]]. The cracking pattern obtained is similar to experimental results shown in Fig. [[#img-29|29]]. The triangular crushed zone below the punch and the axial cleavage crack are reproduced in the numerical simulation.
2047
2048
Quite a good agreement can also be observed between experimental (Fig. [[#img-30|30]]) and numerical (Fig. [[#img-33|33]]) critical values of axial stresses.
2049
2050
<div id='img-32'></div>
2051
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
2052
|-
2053
|[[Image:Draft_Samper_102318540-strip-punch2-bb-035.png|420px|]]
2054
|[[Image:Draft_Samper_102318540-strip-punch2-bb-040.png|426px|]]
2055
|-
2056
|[[Image:Draft_Samper_102318540-strip-punch2-bb-050.png|414px|]]
2057
|[[Image:Draft_Samper_102318540-strip-punch2-bb-1.png|432px|Numerical failure pattern, after <span id='citeF-41'></span>[[#cite-41|[41]]]]]
2058
|- style="text-align: center; font-size: 75%;"
2059
| colspan="2" | '''Figure 32:''' Numerical failure pattern, after <span id='citeF-41'></span>[[#cite-41|[41]]]
2060
|}
2061
2062
<div id='img-33'></div>
2063
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
2064
|-
2065
|[[Image:Draft_Samper_102318540-strip-punch2-f.png|450px|Numerical axial stress vs. axial strain plot]]
2066
|- style="text-align: center; font-size: 75%;"
2067
| colspan="1" | '''Figure 33:''' Numerical axial stress vs. axial strain plot
2068
|}
2069
2070
===7.4 Impact of projectile against rock plate===
2071
2072
The impact of the projectile against a thick rock plate as defined in Fig. [[#img-34|34]] has been studied in order to illustrate possibilities of combined FEM/DEM modelling for this kind of problems.
2073
2074
<div id='img-34'></div>
2075
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
2076
|-
2077
|[[Image:Draft_Samper_102318540-rock-impact.png|276px|Impact of projectile &#8211; definition of geometry]]
2078
|- style="text-align: center; font-size: 75%;"
2079
| colspan="1" | '''Figure 34:''' Impact of projectile &#8211; definition of geometry
2080
|}
2081
2082
The projectile has been discretized with either equal order linear triangular elements using  the FIC formulation previously described, and also   with quadrilateral Q1/P0 elements based on a mixed formulation <span id='citeF-3'></span>[[#cite-3|[3]]]. Projectile has been assumed to be made of steel with the following properties: Young's modulus <math display="inline">E=2\cdot 10^5</math> MPa,  Poisson's ratio <math display="inline">\nu=0.3</math>, yield stress <math display="inline">\sigma _Y=400</math> MPA and isotropic linear hardening modulus <math display="inline">H=300</math> MPA. The macroscopic parameters of rock has been assumed the same as in the previous example (section [[#7.3 Strip punch test|7.3]]). Again  the same collection of discrete elements has been used.
2083
2084
Fracturing of the plate under the impact as well as projectile deformation is shown in Figs. [[#img-35|35]] and [[#img-36|36]] for triangular and quadrilateral projectile discretizations, respectively. Results obtained with quadrilateral elements are considered as the reference ones for the solution obtained with the FIC formulation. The agreement observed between Figs. [[#img-35|35]] and [[#img-36|36]] confirms correctness of the FIC algorithm developed by the authors. Advantages of this algorithm can be observed especially in   3D models with more complicated geometry where  meshing with tetrahedra is  easier than with hexahedra. Effective plastic strains obtained with triangles and quadrilaterals also coincide as it can be seen in Fig. [[#img-37|37]].
2085
2086
Non-symmetry observed in the deformation of the projectile can be attributed in the non-symmetry of the discrete elements in the zone of impact which initiates global buckling of the projectile.
2087
2088
<div id='img-35'></div>
2089
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
2090
|-
2091
|[[Image:Draft_Samper_102318540-rock-bara-ps3t-c2-ini.png|258px|]]
2092
|[[Image:Draft_Samper_102318540-rock-bara-ps3t-c2-bb-024.png|258px|]]
2093
|-
2094
|[[Image:Draft_Samper_102318540-rock-bara-ps3t-c2-bb-040.png|264px|]]
2095
|[[Image:Draft_Samper_102318540-rock-bara-ps3t-c2-bb-056.png|270px|]]
2096
|-
2097
|[[Image:Draft_Samper_102318540-rock-bara-ps3t-c2-bb-080.png|282px|]]
2098
|[[Image:Draft_Samper_102318540-rock-bara-ps3t-c2-bb-112a.png|288px|Impact of projectile &#8211; numerical results obtained using linear triangular elements based on FIC formulation]]
2099
|- style="text-align: center; font-size: 75%;"
2100
| colspan="2" | '''Figure 35:''' Impact of projectile &#8211; numerical results obtained using linear triangular elements based on FIC formulation
2101
|}
2102
2103
<div id='img-36'></div>
2104
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
2105
|-
2106
|[[Image:Draft_Samper_102318540-rock-bara-ps3-c2-ini.png|252px|]]
2107
|[[Image:Draft_Samper_102318540-rock-bara-ps3-c2-bb-024.png|258px|]]
2108
|-
2109
|[[Image:Draft_Samper_102318540-rock-bara-ps3-c2-bb-040.png|264px|]]
2110
|[[Image:Draft_Samper_102318540-rock-bara-ps3-c2-bb-056.png|276px|]]
2111
|-
2112
|[[Image:Draft_Samper_102318540-rock-bara-ps3-c2-bb-080.png|282px|]]
2113
|[[Image:Draft_Samper_102318540-rock-bara-ps3-c2-bb-112.png|288px|Impact of projectile &#8211; numerical results obtained using mixed bilinear Q1/P0 quadrilateral elements]]
2114
|- style="text-align: center; font-size: 75%;"
2115
| colspan="2" | '''Figure 36:''' Impact of projectile &#8211; numerical results obtained using mixed bilinear Q1/P0 quadrilateral elements
2116
|}
2117
2118
<div id='img-37'></div>
2119
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
2120
|-
2121
|[[Image:Draft_Samper_102318540-rock-bara-ps3t-c2-eps-112.png|360px|]]
2122
|[[Image:Draft_Samper_102318540-rock-bara-ps3-c2-eps-112.png|384px|Impact of projectile &#8211; effective plastic strain distribution; comparison of numerical results obtained with FIC-based linear triangles with those obtained using mixed bilinear Q1/P0 quadrilateral elements]]
2123
|- style="text-align: center; font-size: 75%;"
2124
| colspan="2" | '''Figure 37:''' Impact of projectile &#8211; effective plastic strain distribution; comparison of numerical results obtained with FIC-based linear triangles with those obtained using mixed bilinear Q1/P0 quadrilateral elements
2125
|}
2126
2127
===7.5 Study of pipe ovalization===
2128
2129
Interaction of soil and pipe leading to pipe ovalization has been studied. This example demonstrates another possibility of combined FEM/DEM modelling in geomechanics.
2130
2131
The definition of the problem geometry is shown in Fig. [[#img-38|38]]. A pipe of 1m diameter and 1 cm thickness is buried 1 m under the surface of a soil. Due to symmetry only half of the geometry is modelled. The surface load <math display="inline">Q</math> increases linearly as
2132
2133
{| class="formulaSCP" style="width: 100%; text-align: left;" 
2134
|-
2135
| 
2136
{| style="text-align: left; margin:auto;width: 100%;" 
2137
|-
2138
| style="text-align: center;" | <math> Q =k \bar{Q} </math>
2139
|}
2140
|}
2141
2142
where <math display="inline">k</math> is load factor and <math display="inline">\bar{Q}</math> is the constant reference load defined by a set of uniformly distributed concentrated forces <math display="inline">F_i</math>
2143
2144
{| class="formulaSCP" style="width: 100%; text-align: left;" 
2145
|-
2146
| 
2147
{| style="text-align: left; margin:auto;width: 100%;" 
2148
|-
2149
| style="text-align: center;" | <math> \bar{Q}=\{ F_i\} , i=1,12 </math>
2150
|}
2151
|}
2152
2153
The following values of reference forces have been assumed:
2154
2155
{| class="formulaSCP" style="width: 100%; text-align: left;" 
2156
|-
2157
| 
2158
{| style="text-align: left; margin:auto;width: 100%;" 
2159
|-
2160
| style="text-align: center;" | <math> F_1=200 \mbox{N}, \quad  F_2=400 \mbox{N}, \quad  F_3=400 \mbox{N}, \quad  F_4=300 \mbox{N}, \quad  F_5=200 \mbox{N}, \quad  </math>
2161
|}
2162
|}
2163
2164
{| class="formulaSCP" style="width: 100%; text-align: left;" 
2165
|-
2166
| 
2167
{| style="text-align: left; margin:auto;width: 100%;" 
2168
|-
2169
| style="text-align: center;" | <math> F_6=100 \mbox{N}, \quad  F_7=50 \mbox{N}, \quad  F_8=5 \mbox{N}, \quad  F_9=2 \mbox{N}, \quad  F_{10}=1 \mbox{N}, \quad  F_{11}=1 \mbox{N}. </math>
2170
|}
2171
|}
2172
2173
The load has been transmitted to soil by a kind of diaphragm with hinges in the places of force application.
2174
2175
The following material properties for the pipe have been assumed: Young's modulus <math display="inline">E=2\cdot 10^5</math> MPa, Poisson's coefficient <math display="inline">\nu=0.3</math>, density <math display="inline">\rho=7800</math> kg/m<math display="inline">^3</math>,  yield stress <math display="inline">\sigma _Y = 300</math> MPa, isotropic hardening modulus <math display="inline">H=300</math> MPa. The pipe has been discretized with 42 plane strain triangular elements based on the FIC formulation.
2176
2177
The soil specimen has been modelled by 15 000 cylindrical discrete elements of diameters in the range from 2 cm to 3 cm. The following microscopic parameters in the soil model have been used: density <math display="inline">\rho=2800</math> kg/m<math display="inline">^3</math>, the stiffness in the normal and tangential directions <math display="inline">k_n=k_T=200</math>  GPa, Coulomb friction coefficient <math display="inline">\mu=0.8</math>.
2178
2179
The pipe&#8211;soil interface has been defined by the  stiffness in the normal and tangential directions <math display="inline">k_n=k_T=20</math>  GPa,  Coulomb friction coefficient <math display="inline">\mu=0.3</math>.
2180
2181
Numerical results for the deformed pipe shape are shown in Fig. [[#img-39|39]]. The increasing soil pressure leads to ovalization of the pipe. Due to large soil pressure the pipe instability is observed. Pipe material undergoes elasto-plastic deformation. Pipe cross-sections with effective plastic strain distribution at different load level are shown in Fig. [[#img-40|40]].
2182
2183
Pipe cross-sections deformation can be described using an ovalization factor parameter <math display="inline">f_{ov}</math>, defined in terms of the maximum and minimum pipe diameters, <math display="inline">D_{max}</math> and <math display="inline">D_{min}</math>, respectively as
2184
2185
{| class="formulaSCP" style="width: 100%; text-align: left;" 
2186
|-
2187
| 
2188
{| style="text-align: left; margin:auto;width: 100%;" 
2189
|-
2190
| style="text-align: center;" | <math> f_{ov}= \frac{D_{max}-D_{min}}{D_{min}+D_{max}}\cdot 100% </math>
2191
|}
2192
|}
2193
2194
Variation of the ovalization factor for different load levels is shown in Fig. [[#img-41|41]]. Numerical results agree with those given in <span id='citeF-43'></span>[[#cite-43|[43]]].
2195
2196
<div id='img-38'></div>
2197
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
2198
|-
2199
|[[Image:Draft_Samper_102318540-pipeoval.png|492px|Pipe ovalization &#8211; geometry definition]]
2200
|- style="text-align: center; font-size: 75%;"
2201
| colspan="1" | '''Figure 38:''' Pipe ovalization &#8211; geometry definition
2202
|}
2203
2204
<div id='img-39'></div>
2205
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
2206
|-
2207
|[[Image:Draft_Samper_102318540-pipeovaln3s-wd-geo-0.png|480px|]]
2208
|[[Image:Draft_Samper_102318540-pipeovaln3s-wd-geo-02.png|480px|]]
2209
|-
2210
|[[Image:Draft_Samper_102318540-pipeovaln3s-wd-geo-04.png|480px|]]
2211
|[[Image:Draft_Samper_102318540-pipeovaln3s-wd-geo-06.png|480px|Pipe ovalization &#8211; deformed configuration at different load levels]]
2212
|- style="text-align: center; font-size: 75%;"
2213
| colspan="2" | '''Figure 39:''' Pipe ovalization &#8211; deformed configuration at different load levels
2214
|}
2215
2216
<div id='img-40'></div>
2217
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
2218
|-
2219
|[[Image:Draft_Samper_102318540-pipeovaln3s-wd-eps-0.png|246px|]]
2220
|[[Image:Draft_Samper_102318540-pipeovaln3s-wd-eps-02.png|330px|]]
2221
|-
2222
|[[Image:Draft_Samper_102318540-pipeovaln3s-wd-eps-04.png|348px|]]
2223
|[[Image:Draft_Samper_102318540-pipeovaln3s-wd-eps-06.png|366px|Pipe ovalization &#8211; deformed configuration]]
2224
|- style="text-align: center; font-size: 75%;"
2225
| colspan="2" | '''Figure 40:''' Pipe ovalization &#8211; deformed configuration
2226
|}
2227
2228
<div id='img-41'></div>
2229
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
2230
|-
2231
|[[Image:Draft_Samper_102318540-p3s-ovf.png|450px|Pipe ovalization factor]]
2232
|- style="text-align: center; font-size: 75%;"
2233
| colspan="1" | '''Figure 41:''' Pipe ovalization factor
2234
|}
2235
2236
==Conclusions==
2237
2238
The combination of spherical rigid elements and  finite elements is a promising approach for simulation of geomechanical problems involving fracture, penetration and wear. The transient dynamic formulation presented allows to model the motion of discrete and finite elements using a unified algorithm. Frictional contact effects between the rigid spheres and between the spheres and finite elements can be simply accounted for. FIC stabilization procedure allows to use low order triangles and tetrahedra elements with equal order interpolation of the displacement and pressure variables free of the volumetric locking problem. This allows us to apply the presented formulation to problems where large plastic deformation occur in the continuum part of the domain. The discrete formulation can be used to model tools in rock cutting operations. This allows to reproduce tool wear simply by removing the worn particles from the tool surface, thus reproducing the material loss in a realistic manner. Further applications of this approach to the prediction of wear of rock cutting process can be found in <span id='citeF-17'></span><span id='citeF-44'></span>[[#cite-17|[17,44]]].
2239
2240
===BIBLIOGRAPHY===
2241
2242
<div id="cite-1"></div>
2243
'''[[#citeF-1|[1]]]'''  H.&nbsp;Konietzky, editor. ''Numerical Modeling in Micromechanics via Particle   Methods''. Balkema Publishers, 2002.
2244
2245
<div id="cite-2"></div>
2246
'''[[#citeF-2|[2]]]'''  B.K. Cook and R.P. Jensen, editors. ''Discrete Element Methods. Numerical Modeling of Discontinua''. ASCE, 2002.
2247
2248
<div id="cite-3"></div>
2249
'''[[#citeF-3|[3]]]'''  O.C. Zienkiewicz and R.C. Taylor. ''The Finite Element Method''. Butterworth-Heinemann, Oxford, V edition, 2000.
2250
2251
<div id="cite-4"></div>
2252
'''[[#citeF-4|[4]]]'''  E.&nbsp;Oñate, J.&nbsp;Rojek, R.L. Taylor, and O.C. Zienkiewicz. Linear triangles and tetrahedra for incompressible problem using a   finite calculus formulation. In ''Proc. 2nd European Conference on Computational Mechanics   ECCM-2001'', Cracow, 26-29 June, 2001 (on CD ROM), 2001.
2253
2254
<div id="cite-5"></div>
2255
'''[[#citeF-5|[5]]]'''  E.&nbsp;Oñate, R.L. Taylor, O.C. Zienkiewicz, and J.&nbsp;Rojek. A residual correction method based on finite calculus. In ''FEMClass of 42 Meeting'', 30&#8211;31 May, Ibiza   (www.femclass42.com), 2002.
2256
2257
<div id="cite-6"></div>
2258
'''[[#citeF-6|[6]]]'''  P.A. Cundall and O.D.L. Strack. A discrete numerical method for granular assemblies. ''Geotechnique'', 29:47&#8211;65, 1979.
2259
2260
<div id="cite-7"></div>
2261
'''[[#citeF-7|[7]]]'''  C.S. Campbell. Rapid granular flows. ''Annual Review of Fluid Mechanics'', 2:57&#8211;92, 1990.
2262
2263
<div id="cite-8"></div>
2264
'''[[#citeF-8|[8]]]'''  ''Eng. Comput.'', 9(2), 1992. Special issue on Discrete Element Methods, Editor: G. Mustoe.
2265
2266
<div id="cite-9"></div>
2267
'''[[#citeF-9|[9]]]'''  J.R. Williams and R.&nbsp;O'Connor. Discrete Element Simulation and the Contact Problem. ''Archives Comp. Meth. Engng'', 6(4):279&#8211;304, 1999.
2268
2269
<div id="cite-10"></div>
2270
'''[[#citeF-10|[10]]]'''  P.A. Cundall. Formulation of a Three Dimensional Distinct Element Model   &#8211;- Part I. A Scheme to Detect and Represent Contacts in a   System of Many Polyhedral Blocks. ''Int. J. Rock Mech., Min. Sci. & Geomech. Abstr.'',   25(3):107&#8211;116, 1988.
2271
2272
<div id="cite-11"></div>
2273
'''[[#citeF-11|[11]]]'''  J.&nbsp;Rojek, E.&nbsp;Oñate, F.&nbsp;Zarate, and J.&nbsp;Miquel. Modelling of rock, soil and granular materials using spherical   elements. In ''2nd European Conference on Computational Mechanics   ECCM-2001'', Cracow, 26-29 June, 2001.
2274
2275
<div id="cite-12"></div>
2276
'''[[#citeF-12|[12]]]'''  J.&nbsp;Argyris. An excursion into large rotations. ''Comput. Meth. Appl. Mech. Eng.'', 32:85&#8211;155, 1982.
2277
2278
<div id="cite-13"></div>
2279
'''[[#citeF-13|[13]]]'''  L.M. Taylor and D.S. Preece. Simulation of blasting induced rock motion. ''Eng. Comput.'', 9(2):243&#8211;252, 1992.
2280
2281
<div id="cite-14"></div>
2282
'''[[#citeF-14|[14]]]'''  T.J.R. Hughes. ''The Finite Element Method. Linear Static and Dynamic   Analysis''. Prentice-Hall, 1987.
2283
2284
<div id="cite-15"></div>
2285
'''[[#citeF-15|[15]]]'''  J.F. Archard. Contact and rubbing of flat surfaces. ''J. Appl. Phys.'', 24(8):981&#8211;988, 1953.
2286
2287
<div id="cite-16"></div>
2288
'''[[#citeF-16|[16]]]'''  E.&nbsp;Rabinowicz. ''Friction and wear of materials''. 1995.
2289
2290
<div id="cite-17"></div>
2291
'''[[#citeF-17|[17]]]'''  J.&nbsp;Rojek, E.&nbsp;Oñate, F.&nbsp;Zarate, and J.&nbsp;Miquel. Thermomechanical discrete element formulation for wear analysis of   rock cutting tools. (to be published).
2292
2293
<div id="cite-18"></div>
2294
'''[[#citeF-18|[18]]]'''  I.&nbsp;Babuska. The finite element method with lagrangian multipliers. ''Numer. Math.'', 20:179&#8211;192, 1973.
2295
2296
<div id="cite-19"></div>
2297
'''[[#citeF-19|[19]]]'''  F.&nbsp;Brezzi. On the existence, uniqueness and approximation of saddle-point   problems arising from Lagrange multipliers. ''Rev. Francaise d'Automatique Inform. Rech. Opér.,   Ser. Rouge Anal. Numér.'', 8R-2:129&#8211;151, 1974.
2298
2299
<div id="cite-20"></div>
2300
'''[[#citeF-20|[20]]]'''  O.C. Zienkiewicz, S.&nbsp;Qu, R.L. Taylor, and S.&nbsp;Nakazawa. The patch test for mixed formulations. ''Int. J. Num. Meth. Eng.'', 23:1873&#8211;1883, 1986.
2301
2302
<div id="cite-21"></div>
2303
'''[[#citeF-21|[21]]]'''  O.C. Zienkiewicz and R.L. Taylor. The finite element patch test revisited: A computer test for   convergence, validation and error estimates. ''Comput. Meth. Appl. Mech. Eng.'', 149:523&#8211;544, 1997.
2304
2305
<div id="cite-22"></div>
2306
'''[[#citeF-22|[22]]]'''  F.&nbsp;Brezzi and J.&nbsp;Pitkäranta. On the stabilization of finite element approximations of the Stokes   problem. In W.&nbsp;Hackbusch, editor, ''Efficient Solution of Elliptic   Problems, Notes on Numerical Fluid Mechanics''. Vieweg, Wiesbaden, 1984.
2307
2308
<div id="cite-23"></div>
2309
'''[[#citeF-23|[23]]]'''  O.C. Zienkiewicz, J.&nbsp;Rojek, R.L. Taylor, and M.&nbsp;Pastor. Triangles and tetrahedra in explicit dynamic codes for solids. ''Int. J. Num. Meth. Eng.'', 43:565&#8211;583, 1998.
2310
2311
<div id="cite-24"></div>
2312
'''[[#citeF-24|[24]]]'''  J.&nbsp;Rojek, O.C. Zienkiewicz, E.&nbsp;Oñate, and R.L. Taylor. Simulation of metal forming using new formulation of triangular and   tetrahedral elements. In ''<math>8^th</math> Int. Conf. on Metal Forming 2000'', Kraków,   Poland, 2000. Balkema.
2313
2314
<div id="cite-25"></div>
2315
'''[[#citeF-25|[25]]]'''  T.J.R. Hughes, L.P. Franca, and G.M. Hulbert. A new finite element formulation for computational fluid dynamics:   VIII. The galerkin/least-squares method for advective-diffusive   equations. ''Comput. Meth. Appl. Mech. Eng.'', 73:173&#8211;189, 1989.
2316
2317
<div id="cite-26"></div>
2318
'''[[#citeF-26|[26]]]'''  J.&nbsp;Bonet, H.&nbsp;Marriot, and O.&nbsp;Hassan. Stability and comparison of different linear tetrahedral formulations   for nearly incompressible explicit dynamic applications. ''Int. J. Num. Meth. Eng.'', 50:119&#8211;133, 2001.
2319
2320
<div id="cite-27"></div>
2321
'''[[#citeF-27|[27]]]'''  J.&nbsp;Bonet, H.&nbsp;Marriot, and O.&nbsp;Hassan. An average nodal deformation tetrahedron for large strain explicity   dynamic applications. ''Commun. Num. Meth. in Engrg.'', 17:551&#8211;561, 2001.
2322
2323
<div id="cite-28"></div>
2324
'''[[#citeF-28|[28]]]'''  R.&nbsp;Codina and J.&nbsp;Blasco. Stabilized finite element method for transient Navier-Stokes   equations based on pressure gradient projection. ''Comput. Meth. Appl. Mech. Eng.'', 182:287&#8211;300, 2000.
2325
2326
<div id="cite-29"></div>
2327
'''[[#citeF-29|[29]]]'''  R.&nbsp;Codina. Stabilization of incompressibility and convection through orthogonal   sub-scales in finite element methods. ''Comput. Meth. Appl. Mech. Eng.'', 190:1579&#8211;1599, 2000.
2328
2329
<div id="cite-30"></div>
2330
'''[[#citeF-30|[30]]]'''  M.&nbsp;Chiumenti, Q.&nbsp;Valverde, C.&nbsp;Agelet de&nbsp;Saracibar, and M.&nbsp;Cervera. A stabilized formulation for incompressible elasticity using linear   displacement and pressure interpolations. ''Comput. Meth. Appl. Mech. Eng.'', 2002.
2331
2332
<div id="cite-31"></div>
2333
'''[[#citeF-31|[31]]]'''  M.&nbsp;Chiumenti, Q.&nbsp;Valverde, C.&nbsp;Agelet de&nbsp;Saracibar, and M.&nbsp;Cervera. A stabilized formulation for incompressible plasticity using linear   triangles and tetrahedra. ''Int. Journal of Plasticity'', 2002.
2334
2335
<div id="cite-32"></div>
2336
'''[[#citeF-32|[32]]]'''  E.&nbsp;O<math display="inline">\tilde{\mbox{n}}</math>ate. Derivation of stabilized equations for advective-diffusive transport   and fluid flow problems. ''Comput. Meth. Appl. Mech. Eng.'', 151:233&#8211;267, 1998.
2337
2338
<div id="cite-33"></div>
2339
'''[[#citeF-33|[33]]]'''  E.&nbsp;O<math display="inline">\tilde{\mbox{n}}</math>ate. A stabilized finite element method for incompressible flows using a   finite increment calculus formulation. ''Comput. Meth. Appl. Mech. Eng.'', 182:355&#8211;370, 2000.
2340
2341
<div id="cite-34"></div>
2342
'''[[#citeF-34|[34]]]'''  E.&nbsp;O<math display="inline">\tilde{\mbox{n}}</math>ate, J.&nbsp;García, and S.&nbsp;Idelhson. Computation of the stabilization parameter for the finite element   solution of advective-diffusive problems. ''Int. J. Num. Meth. Fluids.'', 25:1385&#8211;1407, 1997.
2343
2344
<div id="cite-35"></div>
2345
'''[[#citeF-35|[35]]]'''  E.&nbsp;O<math display="inline">\tilde{\mbox{n}}</math>ate, J.&nbsp;García, and S.&nbsp;Idelhson. An alpha-adaptive approach for stabilized finite element solution of   advective-diffusive problems. In P.&nbsp;Ladeveze and J.T. Oden, editors, ''New Adv. in Adaptive   Comp. Meth. in Mech.'' Elsevier, 1998.
2346
2347
<div id="cite-36"></div>
2348
'''[[#citeF-36|[36]]]'''  E.&nbsp;O<math display="inline">\tilde{\mbox{n}}</math>ate and M.&nbsp;Manzan. A general procedure for deriving stabilized space-time finite element   methods for advective-diffusive problems. ''Int. J. Num. Meth. Fluids.'', 31:203 &#8211; 221, 1999.
2349
2350
<div id="cite-37"></div>
2351
'''[[#citeF-37|[37]]]'''  E.&nbsp;O<math display="inline">\tilde{\mbox{n}}</math>ate and J.&nbsp;García. A finite element method for fluid&#8211;structure interaction with surface   waves using a finite calculus formulation. ''Comput. Meth. Appl. Mech. Eng.'', 191:635&#8211;660, 2001.
2352
2353
<div id="cite-38"></div>
2354
'''[[#citeF-38|[38]]]'''  J.&nbsp;García and E.&nbsp;O<math display="inline">\tilde{\mbox{n}}</math>ate. An unstructured finite element solver for ship hydrodynamic problems. ''J. Appl. Mech.'', 2002.
2355
2356
<div id="cite-39"></div>
2357
'''[[#citeF-39|[39]]]'''  E.&nbsp;O<math display="inline">\tilde{\mbox{n}}</math>ate, C.&nbsp;Sacco, and S.&nbsp;Idelhson. A finite point method for incompressible flow problems. ''Computing and Visualization in Science'', 3:67&#8211;75, 2000.
2358
2359
<div id="cite-40"></div>
2360
'''[[#citeF-40|[40]]]'''  H.&nbsp;Huang. ''Discrete Element Modeling of Tool-Rock Interaction''. PhD thesis, University of Minnesota, 1999.
2361
2362
<div id="cite-41"></div>
2363
'''[[#citeF-41|[41]]]'''  T.&nbsp;Dede. ''Fracture onset and propagation in layered media''. M.Sc. Dissertation, University of the Witwaterstrand, Johannesburg,   South Africa, 1997.
2364
2365
<div id="cite-42"></div>
2366
'''[[#citeF-42|[42]]]'''  P.A. Klerck. ''The Finite Element Modelling of Discrete Fracture in   Quasi-brittle Materials''. PhD thesis, University of Wales, Swansea, 2000.
2367
2368
<div id="cite-43"></div>
2369
'''[[#citeF-43|[43]]]'''  http://geosim.engr.mun.ca/pipesoil.htm.
2370
2371
<div id="cite-44"></div>
2372
'''[[#citeF-44|[44]]]'''  E.&nbsp;Oñate, C.&nbsp;Recarey, F.&nbsp;Zarate, J.&nbsp;Miquel, and J.&nbsp;Rojek. Characterization of micro-macro parameters in discrete element   formulation. (to be published).
2373
2374
==8 APPENDIX==
2375
2376
===8.1 Matrices and vectors in FEM equations of motion using the FIC formulation===
2377
2378
In FEM equations of motion using the FIC formulation ([[#eq-94|94]])&#8211;([[#eq-96|96]]) the following vectors and matrices have been used:
2379
2380
{| class="formulaSCP" style="width: 100%; text-align: left;" 
2381
|-
2382
| 
2383
{| style="text-align: left; margin:auto;width: 100%;" 
2384
|-
2385
| style="text-align: center;" | <math>{\bf M}_{ij}=\int_{\Omega^e} \rho {\bf N}_i {\bf N}_j\,d\Omega\quad ,\quad {\bf N}_i = N_I {\bf I}_3 </math>
2386
|}
2387
| style="width: 5px;text-align: right;white-space: nowrap;" | (A.1)
2388
|}
2389
2390
is the mass matrix
2391
2392
<span id="eq-2"></span>
2393
{| class="formulaSCP" style="width: 100%; text-align: left;" 
2394
|-
2395
| 
2396
{| style="text-align: left; margin:auto;width: 100%;" 
2397
|-
2398
| style="text-align: center;" | <math>{\bf g}= \int_\Omega {\bf B}^T {\boldsymbol\sigma} d\Omega </math>
2399
|}
2400
| style="width: 5px;text-align: right;white-space: nowrap;" | (A.2)
2401
|}
2402
2403
is the internal nodal force vector and the rest of matrices and vectors are
2404
2405
{| class="formulaSCP" style="width: 100%; text-align: left;" 
2406
|-
2407
| 
2408
{| style="text-align: left; margin:auto;width: 100%;" 
2409
|-
2410
| style="text-align: center;" | <math>{\bf G}_{ij} =  \int_{\Omega^e} ({\boldsymbol\nabla} N_i) N_j \, d\Omega a </math>
2411
|-
2412
| style="text-align: center;" | <math> L_{ij}= \int _{\Omega ^e} {\boldsymbol \nabla }^T N_i [\tau ] {\boldsymbol \nabla } N_j \, d\Omega \qquad ,\quad C_{ij} =  \int _{\Omega ^e} {1\over K} N_i N_j \, d\Omega </math>
2413
|-
2414
| style="text-align: center;" | <math> {\bf C}  =   \left[\begin{matrix} \bar {\bf C}^1 & {\bf 0}&{\bf 0}\\
2415
 {\bf 0}&\bar {\bf C}^2 & {\bf 0}\\ {\bf 0}&{\bf 0}&\bar {\bf C}^3\end{matrix}\right]\qquad ,\quad \bar C_{ij}^k =  \int _{\Omega ^e} \tau _k N_i N_j \, d\Omega </math>
2416
| style="width: 5px;text-align: right;white-space: nowrap;" | (A.3)
2417
|-
2418
| style="text-align: center;" | <math> {\bf Q} = [{\bf Q}^1, {\bf Q}^2, {\bf Q}^3]\quad , \quad Q_{ij}^k = \int _{\Omega ^e} \tau _k {\partial N_i \over \partial x_k} N_j d\Omega </math>
2419
|-
2420
| style="text-align: center;" | <math> {\bf f}_i = \int_{\Omega^e} N_i {\bf b} \, d\Omega +
2421
\int_\Gamma N_i \bar {\bf t} \, d\Gamma \quad , \quad i,j=1, n_d </math>
2422
|}
2423
|}
2424
2425
In above <math display="inline">{\bf b}=[b_1,b_2,b_3]^T</math> and <math display="inline">\bar{\bf t} = [\bar t_1, \bar t_2,\bar t_3]^T</math>,
2426
2427
{| class="formulaSCP" style="width: 100%; text-align: left;" 
2428
|-
2429
| 
2430
{| style="text-align: left; margin:auto;width: 100%;" 
2431
|-
2432
| style="text-align: center;" | <math>{\boldsymbol \nabla } = \left\{\begin{matrix} {\partial  \over \partial x_1}\\ {\partial  \over \partial x_2}\\ {\partial  \over \partial x_3}\\\end{matrix}\right\}\quad , \quad [\tau ]  =\left[\begin{matrix} \tau _{1} &  &  {\bf 0}\\ &  \tau _{2} &\\{\bf 0}&&\tau _{3}\\\end{matrix}\right]\quad , \quad {I}_3=\left[\begin{matrix}1&0&0\\ 0&1&0\\ 0&0&1\\\end{matrix}  \right] </math>
2433
|}
2434
| style="width: 5px;text-align: right;white-space: nowrap;" | (A.4)
2435
|}
2436
2437
<math display="inline">{\bf B}</math> is the standard infinitesimal strain matrix and <math display="inline">{\bf D}_d</math> is the deviatoric constitutive matrix <span id='citeF-3'></span>[[#cite-3|[3]]]. Note that the expression of <math display="inline">{\bf g}</math> of eq. ([[#eq-2|A.2]]) is adequate for non linear structural analysis.
2438

Return to Onate Rojek 2004a.

Back to Top

Document information

Published on 01/01/2004

DOI: 10.1016/j.cma.2003.12.056
Licence: CC BY-NC-SA license

Document Score

0

Times cited: 141
Views 192
Recommendations 0

Share this document

claim authorship

Are you one of the authors of this document?