m (Move page script moved page Draft Samper 102318540 to Onate Rojek 2004a)
 
(37 intermediate revisions by one other user not shown)
Line 1: Line 1:
−
<!-- metadata commented in wiki content
+
Published in ''Comput. Methods Appl. Mech. Engrg.'', Vol 193, 3087-3128, 2004<br />
−
==COMBINATION OF DISCRETE ELEMENT AND FINITE ELEMENT METHODS FOR DYNAMIC ANALYSIS OF GEOMECHANICS PROBLEMS==
+
doi: 10.1016/j.cma.2003.12.056
  
−
'''E. Oñate<math>^\dagger </math> and J. Rojek<math>^{\dagger \ddagger }</math>'''
 
−
 
−
{|
 
−
|-
 
−
|<math>^\dagger </math>International Centre for Numerical Methods in Engineering
 
−
|-
 
−
| Universidad Politécnica de Cataluña
 
−
|-
 
−
| Gran Capitán s/n, 08034 Barcelona, Spain
 
−
|-
 
−
| E&#8211;mail: [mailto:onate@cimne.upc.es onate@cimne.upc.es]
 
−
|-
 
−
| Web page: [http://www.cimne.upc.es www.cimne.upc.es]
 
−
|-
 
−
| <math>^{\ddagger }</math>Institute of Fundamental Technological Research
 
−
|-
 
−
| Polish Academy of Sciences
 
−
|-
 
−
| Swietokrzyska 21, 00-049 Warsaw, Poland
 
−
|-
 
−
| E&#8211;mail: [mailto:jrojek@ippt.gov.pl jrojek@ippt.gov.pl]
 
−
|-
 
−
|
 
−
|}
 
−
-->
 
  
 
==Abstract==
 
==Abstract==
Line 62: Line 37:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| style="text-align: center;" | <math>m_i \ddot{{u}}_i= {{F}}_i\,, </math>
+
| style="text-align: center;" | <math>m_i \ddot{\textbf{u}}_i= {\mathbf{F}}_i\,, </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (1)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (1)
Line 78: Line 53:
 
|}
 
|}
  
−
where <math>{u}</math>is the element centroid displacement in a fixed (inertial) coordinate frame <math>{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>{F}</math>&#8211; the resultant force, and <math>{T}</math>&#8211; the resultant moment about the central axes. Vectors <math>{F}</math>and <math>{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>{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>{X}</math>, and hence  the rotational motion can be easily described in this system.
+
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.
  
 
<div id='img-1'></div>
 
<div id='img-1'></div>
Line 96: Line 71:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| style="text-align: center;" | <math>\ddot{{u}}_i^{n}=\frac{{{F}}_i^n}{m_i}\,,  </math>
+
| style="text-align: center;" | <math>\ddot{{\textbf u}}_i^{n}=\frac{{{\textbf F}}_i^n}{m_i}\,,  </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
Line 107: Line 82:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| style="text-align: center;" | <math>\dot{{u}}_i^{n+1/2}=\dot{{u}}_i^{n-1/2}+\ddot{{u}}_i^{n}{\Delta t}\,,  </math>
+
| 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>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
Line 118: Line 93:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| style="text-align: center;" | <math>{{u}}_i^{n+1}={{u}}_i^{n}+\dot{{u}}_i^{n+1/2}{\Delta t}\,. </math>
+
| style="text-align: center;" | <math>{{\textbf u}}_i^{n+1}={{\textbf u}}_i^{n}+\dot{{\textbf u}}_i^{n+1/2}{\Delta t}\,. </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (5)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (5)
Line 168: Line 143:
 
===2.3 Evaluation of contact forces===
 
===2.3 Evaluation of contact forces===
  
−
====2.3.1 <span id='lb-2.3.1'></span>Decomposition ~of ~the ~contact ~force====
+
====2.3.1 <span id='lb-2.3.1'></span>Decomposition of the contact force====
  
−
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:
+
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">{\textbf F}_1</math> and <math display="inline">{\textbf F}_2</math>, which by the Newton's third law satisfy the following relation:
  
 
{| class="formulaSCP" style="width: 100%; text-align: left;"  
 
{| class="formulaSCP" style="width: 100%; text-align: left;"  
Line 177: Line 152:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| style="text-align: center;" | <math>{F}_{1}=-{F}_2\,. </math>
+
| style="text-align: center;" | <math>{\textbf F}_{1}=-{\textbf F}_2\,. </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (9)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (9)
 
|}
 
|}
  
−
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]])
+
We take <math display="inline">{\textbf F}={\textbf F}_1</math> and decompose <math>{\textbf F}</math>into the normal and tangential components, <math display="inline">{F}_n</math>and <math display="inline">{\textbf F}_T</math>, respectively (Fig. [[#img-2|2]])
  
 
{| class="formulaSCP" style="width: 100%; text-align: left;"  
 
{| class="formulaSCP" style="width: 100%; text-align: left;"  
Line 189: Line 164:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| style="text-align: center;" | <math> {F}={F}_n+{F}_T =F_n{n}+{F}_T \,, </math>
+
| style="text-align: center;" | <math> {\textbf F}={\textbf F}_n+{\textbf F}_T =F_n{\textbf n}+{\textbf F}_T \,, </math>
 
|}
 
|}
 
|}
 
|}
  
−
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.
+
where <math>{\textbf 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.
  
−
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>.
+
The contact forces <math display="inline">F_n</math> and <math display="inline">{\textbf 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>.
  
 
<div id='img-2'></div>
 
<div id='img-2'></div>
Line 213: Line 188:
 
|}
 
|}
  
−
====2.3.2 <span id='lb-2.3.2'></span>Normal ~contact ~force====
+
====2.3.2 <span id='lb-2.3.2'></span>Normal contact force====
  
 
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>
 
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>
Line 290: Line 265:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| style="text-align: center;" | <math>{v}_{rn}=(\dot{{u}}_2-\dot{{u}}_1)\cdot { n}\,. </math>
+
| style="text-align: center;" | <math>{v}_{rn}=(\dot{{\textbf u}}_2-\dot{{\textbf u}}_1)\cdot {\textbf  n}\,. </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (15)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (15)
Line 308: Line 283:
 
|}
 
|}
  
−
====2.3.3 <span id='lb-2.3.3'></span>Tangential ~frictional ~contact====
+
====2.3.3 <span id='lb-2.3.3'></span>Tangential frictional contact====
  
−
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
+
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">{\textbf 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">{\textbf v}_{rT}</math> is calculated from the following relationship
  
 
<span id="eq-17"></span>
 
<span id="eq-17"></span>
Line 318: Line 293:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| style="text-align: center;" | <math>{v}_{rT}= {v}_{r}-{v}_{r}\cdot {n}\,, </math>
+
| style="text-align: center;" | <math>{\textbf v}_{rT}= {\textbf v}_{r}-{\textbf v}_{r}\cdot {\textbf n}\,, </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (17)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (17)
 
|}
 
|}
−
 
+
with
 
{| class="formulaSCP" style="width: 100%; text-align: left;"  
 
{| class="formulaSCP" style="width: 100%; text-align: left;"  
 
|-
 
|-
Line 328: Line 303:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| 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>
+
| style="text-align: center;" | <math>{\textbf v}_r=(\dot{{\textbf u}}_{2} +{\omega }_{2} \times {\textbf r}_{c2})- (\dot{{\textbf u}}_{1} +{\omega }_1 \times {\textbf r}_{c1})\,, </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (18)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (18)
 
|}
 
|}
  
−
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.
+
where <math display="inline">\dot{{\textbf u}}_{1}</math>, <math display="inline">\dot{{\textbf 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">{\textbf  r}_{c1}</math> and <math display="inline">{\textbf r}_{c2}</math> are the vectors connecting particle centres with contact points.
  
 
<div id='img-4'></div>
 
<div id='img-4'></div>
Line 344: Line 319:
 
|}
 
|}
  
−
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:
+
The relationship between the friction force <math display="inline">\| {\textbf 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">{\textbf v}_{rT}^{\mathrm{r}}</math> and <math display="inline">{\textbf v}_{rT}^{\mathrm{ir}}</math>, respectively as:
  
 
{| class="formulaSCP" style="width: 100%; text-align: left;"  
 
{| class="formulaSCP" style="width: 100%; text-align: left;"  
Line 351: Line 326:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| style="text-align: center;" | <math>{v}_{rT}= {v}_{rT}^{\mathrm{r}}+{v}_{r}^{\mathrm{ir}}\,. </math>
+
| style="text-align: center;" | <math>{\textbf v}_{rT}= {\textbf v}_{rT}^{\mathrm{r}}+{\textbf v}_{r}^{\mathrm{ir}}\,. </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (19)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (19)
Line 364: Line 339:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| style="text-align: center;" | <math>{F}_T^{\mathrm{trial}}={F}_T^{\mathrm{old}}-k_T {v}_{rT} \Delta t\,,  </math>
+
| style="text-align: center;" | <math>{\textbf F}_T^{\mathrm{trial}}={\textbf F}_T^{\mathrm{old}}-k_T {\textbf v}_{rT} \Delta t\,,  </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (20)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (20)
Line 376: Line 351:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| style="text-align: center;" | <math>\phi ^{\mathrm{trial}}=\| {F}_T^{\mathrm{trial}}\|{-\mu}|F_n|\,. </math>
+
| style="text-align: center;" | <math>\phi ^{\mathrm{trial}}=\| {\textbf F}_T^{\mathrm{trial}}\|{-\mu}|F_n|\,. </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (21)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (21)
Line 388: Line 363:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| style="text-align: center;" | <math>{F}_T^{\mathrm{new}}={F}_T^{\mathrm{trial}}\,, </math>
+
| style="text-align: center;" | <math>{\textbf F}_T^{\mathrm{new}}={\textbf F}_T^{\mathrm{trial}}\,, </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (22)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (22)
Line 400: Line 375:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| style="text-align: center;" | <math>{F}_T^{\mathrm{new}}=\mu |F_n|\frac{{F}_T^{\mathrm{trial}}}{\| {F}_T^{\mathrm{trial}} \| }\,. </math>
+
| style="text-align: center;" | <math>{\textbf F}_T^{\mathrm{new}}=\mu |F_n|\frac{{\textbf F}_T^{\mathrm{trial}}}{\| {\textbf F}_T^{\mathrm{trial}} \| }\,. </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (23)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (23)
Line 415: Line 390:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| style="text-align: center;" | <math>m_i \ddot{{u}}_i= {{F}}_i+{{F}}_i^{\mathrm{damp}}\,, </math>
+
| style="text-align: center;" | <math>m_i \ddot{\mathbf{u}}_i= {\mathbf{F}}_i+{\mathbf{F}}_i^{\mathrm{damp}}\,, </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (24)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (24)
Line 443: Line 418:
 
| style="text-align: center;" | <math>
 
| style="text-align: center;" | <math>
  
−
{{F}}_{i}^{\mathrm{damp}}=-\alpha ^{vt}m_{i}\dot{{ u}}_{i}\,,  </math>
+
{\mathbf{F}}_{i}^{\mathrm{damp}}=-\alpha ^{vt}m_{i}\dot{\mathbf{ u}}_{i}\,,  </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (26)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (26)
Line 471: Line 446:
 
| style="text-align: center;" | <math>
 
| style="text-align: center;" | <math>
  
−
{{F}}_{i}^{\mathrm{damp}}=-\alpha ^{nvt}\Vert {{F}}_{i}\Vert \frac{\dot{{u}}_{i}}{\Vert \dot{{u}}_{i}\Vert }\,,  </math>
+
{\mathbf{F}}_{i}^{\mathrm{damp}}=-\alpha ^{nvt}\Vert {\mathbf{F}}_{i}\Vert \frac{\dot{\mathbf{u}}_{i}}{\Vert \dot{\mathbf{u}}_{i}\Vert }\,,  </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (28)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (28)
Line 778: Line 753:
 
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
 
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
 
|-
 
|-
−
|[[Image:Draft_Samper_102318540-brittlesurf.png|273px|Failure surface for the elastic perfectly brittle model]]
+
|[[Image:Draft_Samper_102318540-brittlesurf.png|473px|Failure surface for the elastic perfectly brittle model]]
 
|- style="text-align: center; font-size: 75%;"
 
|- style="text-align: center; font-size: 75%;"
 
| colspan="1" | '''Figure 7:''' Failure surface for the elastic perfectly brittle model
 
| colspan="1" | '''Figure 7:''' Failure surface for the elastic perfectly brittle model
Line 817: Line 792:
 
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>.
 
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>.
  
−
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.
+
An elasto-plastic contact law with linear softening is assumed for the shear and normal tensile contact forces (Figs. [[#img-8|8]] and [[#img-9|9]]). An elastic linear law is assumed for compression.
  
 
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.
 
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.
Line 857: Line 832:
 
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
 
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
 
|-
 
|-
−
|[[Image:Draft_Samper_102318540-yielsurf.png|231px|Yield surface for elasto-plastic model]]
+
|[[File:Draft_Samper_102318540_6786_Fig10.png]]
 
|- style="text-align: center; font-size: 75%;"
 
|- style="text-align: center; font-size: 75%;"
 
| colspan="1" | '''Figure 10:''' Yield surface for elasto-plastic model
 
| colspan="1" | '''Figure 10:''' Yield surface for elasto-plastic model
Line 908: Line 883:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| style="text-align: center;" | <math>\psi (u_n)=\left\{ \begin{array}{lll}1 & {for} & \dst u_n\le \frac{R_n}{k_n}\\[2ex] \dst \frac{k^2_n u_n}{HR_n +k_n R_n-Hk_n u_n} & {for}                  & \dst \frac{R_n}{k_n}\le u_n\le \frac{R_n}{k_n}+\frac{R_n}{H}\\[2ex] \dst \infty  & {for}                  & \dst u_n\ge \frac{R_n}{k_n}+\frac{R_n}{H} \end{array} \right.  </math>
+
| 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>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (56)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (56)
Line 998: Line 973:
 
|}
 
|}
  
−
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:
+
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 <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:
  
 
<span id="eq-62"></span>
 
<span id="eq-62"></span>
Line 1,065: Line 1,040:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| style="text-align: center;" | <math>\ri - \underline{\frac{\hk }{2}{\partial \ri  \over \partial \xk }}\frac{\hk }{2}{\partial \ri  \over \partial \xk }=0 \quad \mbox{in} \;\,\Omega \qquad k=1,n_d </math>
+
| 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>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (65)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (65)
Line 1,078: Line 1,053:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| style="text-align: center;" | <math>\ri : ={\partial \sij  \over \partial \xj }+\bi  </math>
+
| style="text-align: center;" | <math>{r}_{i}: ={\partial {\sigma }_{ij} \over \partial {x}_j}+{b}_i </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (66)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (66)
 
|}
 
|}
  
−
In  ([[#eq-65|65]]) and ([[#eq-66|66]]) and are the stresses and the body forces, respectively and are characteristic length distances  of an arbitrary prismatic domain where equilibrium of forces is considered.
+
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.
  
 
Equations ([[#eq-65|65]]) and ([[#eq-66|66]]) are completed with the boundary conditions on the displacements <math display="inline">u_i</math>
 
Equations ([[#eq-65|65]]) and ([[#eq-66|66]]) are completed with the boundary conditions on the displacements <math display="inline">u_i</math>
Line 1,093: Line 1,068:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| style="text-align: center;" | <math>\ui -\bui =0\quad \mbox{on} \;\,\Gu  </math>
+
| style="text-align: center;" | <math>{u}_i-\bar{u}_i=0\quad \mbox{on} \;\,\Gamma _u </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (67)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (67)
Line 1,106: Line 1,081:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| style="text-align: center;" | <math>\sij \nj - \bti - \underline{\frac{1}{2}\hk \nk \ri }\frac{1}{2}\hk \nk \ri  =0\quad \mbox{on} \;\,\Gt  </math>
+
| 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>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (68)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (68)
 
|}
 
|}
  
−
In the above and are prescribed displacements and tractions over the boundaries <math display="inline">\Gu </math>and <math display="inline">\Gt </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.
+
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.
  
 
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]]].
 
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]]].
Line 1,127: Line 1,102:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| style="text-align: center;" | <math>\sij = s_{ij} +p\dij  </math>
+
| style="text-align: center;" | <math>{\sigma }_{ij}= s_{ij} +p{\delta }_{ij} </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (69)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (69)
 
|}
 
|}
  
−
where is the Kronnecker delta function. The linear elastic constitutive equations for the deviatoric stresses <math display="inline">s_{ij}</math> are written as
+
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
  
 
<span id="eq-70"></span>
 
<span id="eq-70"></span>
Line 1,140: Line 1,115:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| style="text-align: center;" | <math>s_{ij}\!\!=\!\! 2G \left(\eij -\frac{1}{3}\ev \dij \right) </math>
+
| style="text-align: center;" | <math>s_{ij}\!\!=\!\! 2G \left({\varepsilon }_{ij}-\frac{1}{3}{\varepsilon }_{v}{\delta }_{ij}\right) </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (70)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (70)
Line 1,152: Line 1,127:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| style="text-align: center;" | <math>\eij = \frac{1}{2} \left( {\partial \ui  \over \partial \xj }+ {\partial \uj  \over \partial x_i} \right) \quad  \mbox{and} \quad  \ev =\eii \,. </math>
+
| 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>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (71)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (71)
Line 1,165: Line 1,140:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| style="text-align: center;" | <math>\left({1\over K} p - \ev \right)- \frac{\hk }{2}{\partial  \over \partial \xk } \left({1\over K} p - \ev \right)=0\, ,\quad k=1,2,3\,\,\, \hbox{for }3D </math>
+
| 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>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (72)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (72)
Line 1,182: Line 1,157:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| style="text-align: center;" | <math>\varepsilon _v - \frac{\hk }{2}{\partial \ev  \over \partial \xk }=0 </math>
+
| style="text-align: center;" | <math>\varepsilon _v - \frac{{h}_k}{2}{\partial {\varepsilon }_{v} \over \partial {x}_k}=0 </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (73)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (73)
Line 1,197: Line 1,172:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| style="text-align: center;" | <math>{\partial s_{ij} \over \partial \xj }+{\partial p \over \partial x_i}+\bi -\frac{\hk }{2} {\partial r_i \over \partial \xk }=0 </math>
+
| 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>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (74)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (74)
Line 1,208: Line 1,183:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| style="text-align: center;" | <math>\left(\frac{p}{K}-\ev \right)-\frac{\hk }{2}{\partial  \over \partial \xk }\left(\frac{p}{K}-\ev \right)=0 </math>
+
| 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>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (75)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (75)
Line 1,221: Line 1,196:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| style="text-align: center;" | <math>{\partial \ev  \over \partial \ixi }= \frac{3}{2G}\left[\hri - \frac{\hk }{2}{\partial \ri  \over \partial \xk } \right] </math>
+
| 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>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (76)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (76)
 
|}
 
|}
  
−
where is defined by Eq. ([[#eq-66|66]]) and
+
where <math>{r}_{i}</math>is defined by Eq. ([[#eq-66|66]]) and
  
 
<span id="eq-77"></span>
 
<span id="eq-77"></span>
Line 1,234: Line 1,209:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| style="text-align: center;" | <math>\hri =\frac{\partial }{\partial x_j}(2G \varepsilon _{ij})+{\partial p \over \partial x_i} + b_i </math>
+
| 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>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (77)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (77)
Line 1,247: Line 1,222:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| style="text-align: center;" | <math>\left(\frac{p}{K}-\ev \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}{\hk \over G}{\partial r_i \over \partial \xk }\right)=0 </math>
+
| 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>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (78)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (78)
Line 1,258: Line 1,233:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| style="text-align: center;" | <math>{p\over K} - \ev - \sum \limits _{i=1}^{n_d} \tau _i {\partial r_i \over \partial x_i}=0 </math>
+
| 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>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (79)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (79)
Line 1,289: Line 1,264:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| style="text-align: center;" | <math>r_i:= - \rho \ppder{u_i}{t} +  {\partial \sigma _{ij} \over \partial x_j}+ b_i </math>
+
| 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>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (81)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (81)
Line 1,304: Line 1,279:
 
The set of stabilized equations to be solved are now:
 
The set of stabilized equations to be solved are now:
  
−
===Equilibrium===
+
'''''Equilibrium'''''
  
 
<span id="eq-82"></span>
 
<span id="eq-82"></span>
Line 1,317: Line 1,292:
 
|}
 
|}
  
−
===Pressure~constitutive ~equation===
+
'''''Pressure constitutive equation'''''
  
 
<span id="eq-83"></span>
 
<span id="eq-83"></span>
Line 1,336: Line 1,311:
 
The weighted residual form of the FIC governing equations ([[#eq-82|82]]), ([[#eq-68|68]]) and ([[#eq-83|83]]) is
 
The weighted residual form of the FIC governing equations ([[#eq-82|82]]), ([[#eq-68|68]]) and ([[#eq-83|83]]) is
  
−
===Equilibrium===
+
'''''Equilibrium''''
  
 
<span id="eq-84"></span>
 
<span id="eq-84"></span>
Line 1,344: Line 1,319:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| style="text-align: center;" | <math>\int _\Omega \!\delta u_i \left[-\rho \ppder{u_i}{t}+ {\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>
+
| 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>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (84)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (84)
 
|}
 
|}
  
−
===Pressure~constitutive ~equation===
+
'''''Pressure constitutive equation'''''
  
 
<span id="eq-85"></span>
 
<span id="eq-85"></span>
Line 1,357: Line 1,332:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| style="text-align: center;" | <math>\int _\Omega q \left({p\over K} - \ev \right)d\Omega - \int _\Omega q \left(\sum \limits _{i=1}^{n_d} \tau _i {\partial \ri  \over \partial x_i}\right)\, d\Omega =0 </math>
+
| 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>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (85)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (85)
Line 1,372: Line 1,347:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| style="text-align: center;" | <math>\IV \delta u_i \rho \ppder{u_i}{t}d\Omega +  \IV \!\delta \eij \sigma _{ij}\,d\Omega{-} \IV \delta \ui \bi \, d\Omega - \IGt \delta \ui \bti \,d\Omega =0 </math>
+
| 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>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (86)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (86)
Line 1,383: Line 1,358:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| style="text-align: center;" | <math>\int _\Omega q \left({p\over K} - \ev \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>
+
| 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>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (87)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (87)
Line 1,407: Line 1,382:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| style="text-align: center;" | <math>\pi _i = - \ppder{u_i}{t} + {\partial s_{ij} \over \partial x_i}+b_i </math>
+
| 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>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (89)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (89)
Line 1,424: Line 1,399:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| style="text-align: center;" | <math>\IV \delta u_i \rho \ppder{u_i}{t}\,d\Omega + \IV \delta \eij \sij \,d\Omega{-} \IV \delta u_i b_i \,d\Omega  - \IGt \delta \ui \bti  \,d\Gamma _t=0 </math>
+
| 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>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (90)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (90)
Line 1,435: Line 1,410:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| style="text-align: center;" | <math>\IV 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>
+
| 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>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (91)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (91)
Line 1,446: Line 1,421:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| style="text-align: center;" | <math>\IV \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>
+
| 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>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (92)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (92)
Line 1,461: Line 1,436:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| style="text-align: center;" | <math>\ui =\sum _{j=1}^{n} N_j \bui ^j \quad , \quad p=\sum _{j=1}^{n} N_j \bp ^j \quad , \quad  \pi _i =\sum _{j=1}^{n} N_j \bar \pi ^j_i </math>
+
| 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>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (93)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (93)
Line 1,474: Line 1,449:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| style="text-align: center;" | <math>{M} \ddot{\bar {u}} + {g} -  {f}={0} </math>
+
| style="text-align: center;" | <math>\mathbf{M} \ddot{\bar \mathbf{u}} + \mathbf{g} -  \mathbf{f}=\mathbf{0} </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (94)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (94)
Line 1,485: Line 1,460:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| style="text-align: center;" | <math>{G}^T \Delta \bar {u} -{C} \Delta \bar {p}- {L} \bar {p} - {Q}\bar {\boldsymbol \pi }= {0}  </math>
+
| style="text-align: center;" | <math>\mathbf{G}^T \Delta \bar \mathbf{u} -\mathbf{C} \Delta \bar \mathbf{p}- \mathbf{L} \bar \mathbf{p} - \mathbf{Q}\bar {\boldsymbol \pi }= {0}  </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (95)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (95)
Line 1,496: Line 1,471:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| style="text-align: center;" | <math>{Q}^T \bar {p}+ \bar {C} \bar {\boldsymbol \pi }=  {0}  </math>
+
| style="text-align: center;" | <math>\mathbf{Q}^T \bar \mathbf{p}+ \bar \mathbf{C} \bar {\boldsymbol \pi }=  \mathbf{0}  </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (96)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (96)
 
|}
 
|}
  
−
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.
+
where <math display="inline">\ddot{\bar \mathbf{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.
  
 
A four steps semi-implicit time integration algorithm can be derived from eqs. ([[#eq-94|94]])&#8211;([[#eq-96|96]]) as follows
 
A four steps semi-implicit time integration algorithm can be derived from eqs. ([[#eq-94|94]])&#8211;([[#eq-96|96]]) as follows
Line 1,507: Line 1,482:
 
step
 
step
  
−
'''Step step.'''step
+
'''Step 1'''. Compute the nodal velocities <math display="inline">\dot{\bar \mathbf{u}}^{n+1/2}</math>
−
* Compute the nodal velocities <math display="inline">\dot{\bar {u}}^{n+1/2}</math>
+
  
 
<span id="eq-97"></span>
 
<span id="eq-97"></span>
Line 1,518: Line 1,492:
 
| style="text-align: center;" | <math>
 
| 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>
+
\dot{\bar \mathbf{u}}^{n+1/2} = \dot{\bar \mathbf{u}}^{n-1/2}+\Delta t \mathbf{M}^{-1}_d (\mathbf{f}^n - \mathbf{g}^n)  </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (97)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (97)
 
|}
 
|}
  
−
* Compute the nodal displacements <math display="inline">\bar {u}^{n+1}</math>
+
'''Step 2'''. Compute the nodal displacements <math display="inline">\bar \mathbf{u}^{n+1}</math>
  
 
<span id="eq-98"></span>
 
<span id="eq-98"></span>
Line 1,533: Line 1,507:
 
| style="text-align: center;" | <math>
 
| style="text-align: center;" | <math>
  
−
\bar {u}^{n+1} = \bar{u}^n + \Delta t \dot{\bar {u}}^{n+1/2}  </math>
+
\bar \mathbf{u}^{n+1} = \bar\mathbf{u}^n + \Delta t \dot{\bar \mathbf{u}}^{n+1/2}  </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (98)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (98)
 
|}
 
|}
−
* Compute the nodal pressures <math display="inline">\bar {p}^{n+1}</math>
+
 
 +
'''Step 3'''. Compute the nodal pressures <math display="inline">\bar \mathbf{p}^{n+1}</math>
  
 
<span id="eq-99"></span>
 
<span id="eq-99"></span>
Line 1,547: Line 1,522:
 
| style="text-align: center;" | <math>
 
| style="text-align: center;" | <math>
  
−
\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>
+
\bar \mathbf{p}^{n+1}=[ \mathbf{C} + \mathbf{L}]^{-1} [\Delta t \mathbf{G}^T \dot{\bar \mathbf{u}}^{n+1/2}+  \mathbf{C} \bar \mathbf{p}^n- \mathbf{Q} \bar {\boldsymbol \pi }^n]  </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (99)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (99)
 
|}
 
|}
−
* Compute the nodal projected pressure gradients <math display="inline">\bar {\boldsymbol \pi }^{n+1}</math>
+
 
 +
'''Step 4'''. Compute the nodal projected pressure gradients <math display="inline">\bar {\boldsymbol \pi }^{n+1}</math>
  
 
<span id="eq-100"></span>
 
<span id="eq-100"></span>
Line 1,561: Line 1,537:
 
| style="text-align: center;" | <math>
 
| style="text-align: center;" | <math>
  
−
\bar {\boldsymbol \pi }^{n+1}=- \bar {C}_d^{-1} {Q}^T \bar {p}^{n+1}  </math>
+
\bar {\boldsymbol \pi }^{n+1}=- \bar \mathbf{C}_d^{-1} \mathbf{Q}^T \bar \mathbf{p}^{n+1}  </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (100)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (100)
Line 1,573: Line 1,549:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| style="text-align: center;" | <math>{g}^n =\int _{\Omega ^e} [{B}^T {\boldsymbol \sigma }]^n\,d\Omega  </math>
+
| style="text-align: center;" | <math>\mathbf{g}^n =\int _{\Omega ^e} [\mathbf{B}^T {\boldsymbol \sigma }]^n\,d\Omega  </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (101)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (101)
Line 1,580: Line 1,556:
 
where the stresses <math display="inline">{\boldsymbol \sigma }^n</math> are obtained by consistent integration of the adequate (non linear)  constitutive law [33].
 
where the stresses <math display="inline">{\boldsymbol \sigma }^n</math> are obtained by consistent integration of the adequate (non linear)  constitutive law [33].
  
−
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).
+
Note that steps 1, 2 and 4 are fully explicit as a diagonal form of matrices <math display="inline">\mathbf{C}</math> and <math display="inline">\bar\mathbf{C}</math> has been chosen. The solution of step 3 with  a diagonal form for <math display="inline">\mathbf{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).
  
 
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
 
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
Line 1,589: Line 1,565:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| 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>
+
| style="text-align: center;" | <math>\bar\mathbf{p}^{n+1}= [\mathbf{C} + \mathbf{L}-\mathbf{S}]^{-1} [\Delta t \mathbf{G}^T \dot{\bar \mathbf{u}}^{n+1/2} + \mathbf{C} \bar\mathbf{p}^n ] </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (102)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (102)
Line 1,601: Line 1,577:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| style="text-align: center;" | <math>{S} = {Q} \bar {C}_d^{-1} {Q}^T </math>
+
| style="text-align: center;" | <math>\mathbf{S} = \mathbf{Q} \bar \mathbf{C}_d^{-1} \mathbf{Q}^T </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (103)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (103)
 
|}
 
|}
  
−
Recall that for the full incompressible case <math display="inline">K=\infty </math> and <math display="inline">{C}=0 </math> in all above equations.
+
Recall that for the full incompressible case <math display="inline">K=\infty </math> and <math display="inline">\mathbf{C}=0 </math> in all above equations.
  
 
The critical time step <math display="inline">\Delta t</math> is taken as that of the standard explicit dynamic scheme [28].
 
The critical time step <math display="inline">\Delta t</math> is taken as that of the standard explicit dynamic scheme [28].
  
−
===Explicit algorithm===
+
'''''Explicit algorithm'''''
  
−
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
+
A fully explicit algorithm can be obtained by computing <math display="inline">\bar \mathbf{p}^{n+1}</math> from step 3 in Eq. ([[#eq-99|99]])  as follows
  
 
<span id="eq-104"></span>
 
<span id="eq-104"></span>
Line 1,620: Line 1,596:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| 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>
+
| style="text-align: center;" | <math>\bar\mathbf{p}^{n+1}=\mathbf{C}_d^{-1}[\Delta t \mathbf{G}^T \dot{\bar \mathbf{u}}^{n+1/2} + (\mathbf{C}_d+\mathbf{L}) \bar\mathbf{p}^n - \mathbf{Q} {\boldsymbol \sigma }^n] </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (104)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (104)
Line 1,629: Line 1,605:
 
==5 DEM&#8211;FEM CONTACT ALGORITHM==
 
==5 DEM&#8211;FEM CONTACT ALGORITHM==
  
−
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.
+
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">\mathbf{F}</math>  of a finite element is decomposed into normal and tangential components, <math display="inline">\mathbf{F}_n</math> and <math display="inline">\mathbf{ 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.
  
 
<div id='img-13'></div>
 
<div id='img-13'></div>
Line 1,697: Line 1,673:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| style="text-align: center;" | <math>{v}_{r}= ({{v}}_{C} +{\omega }_{C} \times {r})- ({{v}}_{1}H_{1} + {{v}}_{2}H_{2}) \,, </math>
+
| style="text-align: center;" | <math>\mathbf{v}_{r}= ({\mathbf{v}}_{C} +{\omega }_{C} \times \mathbf{r})- ({\mathbf{v}}_{1}H_{1} + {\mathbf{v}}_{2}H_{2}) \,, </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (109)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (109)
 
|}
 
|}
  
−
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]]).
+
where <math display="inline">{\mathbf{v}}_{C} +{\omega }_{C} \times \mathbf{r}</math> is the discrete element velocity at the contact point and <math display="inline">{\mathbf{v}}_{1}H_{1} + {\mathbf{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">{\mathbf{v}}_{1}</math> and <math display="inline">{\mathbf{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]]).
  
 
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]].
 
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]].
Line 1,710: Line 1,686:
 
The general algorithm for the transient dynamic problem involving discrete elements and (stabilized) finite elements includes the following steps.
 
The general algorithm for the transient dynamic problem involving discrete elements and (stabilized) finite elements includes the following steps.
  
−
===Step1. Compute the nodal velocities===
+
====Step1. Compute the nodal velocities====
  
−
====Discrete elements====
+
'''Discrete elements'''
  
 
<span id="eq-110"></span>
 
<span id="eq-110"></span>
Line 1,720: Line 1,696:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| style="text-align: center;" | <math>\dot{{u}}_i^{n+1/2}=\dot{{u}}_i^{n-1/2}+\ddot{{u}}_i^{n}{\Delta t}\,, </math>
+
| style="text-align: center;" | <math>\dot{\mathbf{u}}_i^{n+1/2}=\dot{\mathbf{u}}_i^{n-1/2}+\ddot{\mathbf{u}}_i^{n}{\Delta t}\,, </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (110)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (110)
Line 1,733: Line 1,709:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| style="text-align: center;" | <math>\ddot{{u}}_i^{n}=\frac{{{F}}_i^n}{m_i}\,, </math>
+
| style="text-align: center;" | <math>\ddot{\mathbf{u}}_i^{n}=\frac{{\mathbf{F}}_i^n}{m_i}\,, </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (111)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (111)
Line 1,757: Line 1,733:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| style="text-align: center;" | <math>\dot{{\omega }}_i^{n}=\frac{{{T}}_i^n}{I_i}\,, </math>
+
| style="text-align: center;" | <math>\dot{{\omega }}_i^{n}=\frac{{\mathbf{T}}_i^n}{I_i}\,, </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (113)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (113)
 
|}
 
|}
  
−
====Finite elements====
+
'''Finite elements'''
  
 
<span id="eq-114"></span>
 
<span id="eq-114"></span>
Line 1,770: Line 1,746:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| 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>
+
| style="text-align: center;" | <math>\dot{\bar \mathbf{u}}^{n+1/2} = \dot{\bar \mathbf{u}}^{n-1/2}+\Delta t \mathbf{M}^{-1}_d (\mathbf{f}^n - \mathbf{g}^n)  </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (114)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (114)
 
|}
 
|}
  
−
===Step2. Compute the nodal displacements===
+
====Step2. Compute the nodal displacements====
  
−
====Discrete elements====
+
'''Discrete elements'''
  
 
<span id="eq-115"></span>
 
<span id="eq-115"></span>
Line 1,785: Line 1,761:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| style="text-align: center;" | <math>{{u}}_i^{n+1}={{u}}_i^{n}+\dot{{u}}_i^{n+1/2}{\Delta t}\,. </math>
+
| style="text-align: center;" | <math>{\mathbf{u}}_i^{n+1}={\mathbf{u}}_i^{n}+\dot{\mathbf{u}}_i^{n+1/2}{\Delta t}\,. </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (115)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (115)
Line 1,801: Line 1,777:
 
|}
 
|}
  
−
====Finite elements====
+
'''Finite elements'''
  
 
<span id="eq-117"></span>
 
<span id="eq-117"></span>
Line 1,809: Line 1,785:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| style="text-align: center;" | <math>\bar {u}^{n+1} = \bar{u}^n + \Delta t \dot{\bar {u}}^{n+1/2}  </math>
+
| style="text-align: center;" | <math>\bar \mathbf{u}^{n+1} = \bar\mathbf{u}^n + \Delta t \dot{\bar \mathbf{u}}^{n+1/2}  </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (117)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (117)
 
|}
 
|}
  
−
===Step3. Compute the nodal pressures===
+
====Step3. Compute the nodal pressures====
  
−
====Finite elements====
+
'''Finite elements'''
  
 
{| class="formulaSCP" style="width: 100%; text-align: left;"  
 
{| class="formulaSCP" style="width: 100%; text-align: left;"  
Line 1,823: Line 1,799:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| 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>
+
| style="text-align: center;" | <math>\bar\mathbf{p}^{n+1}= [\mathbf{C} + \mathbf{L}-\mathbf{S}]^{-1} [\Delta t \mathbf{G}^T \dot{\bar \mathbf{u}}^{n+1/2} + \mathbf{C} \bar\mathbf{p}^n ] </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (118)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (118)
 
|}
 
|}
  
−
===Step4. Update the nodal coordinates===
+
====Step4. Update the nodal coordinates====
  
−
====Discrete and finite elements====
+
'''Discrete and finite elements'''
  
 
{| class="formulaSCP" style="width: 100%; text-align: left;"  
 
{| class="formulaSCP" style="width: 100%; text-align: left;"  
Line 1,837: Line 1,813:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| style="text-align: center;" | <math>{X}_i^{n+1} = {X}_i^{n}+\Delta \bar {u}_i </math>
+
| style="text-align: center;" | <math>\mathbf{X}_i^{n+1} = \mathbf{X}_i^{n}+\Delta \bar \mathbf{u}_i </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (119)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (119)
Line 1,847: Line 1,823:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| style="text-align: center;" | <math>\Delta \bar {u}_i =\bar {u}_i^{n+1} - \bar {u}_i^n </math>
+
| style="text-align: center;" | <math>\Delta \bar \mathbf{u}_i =\bar \mathbf{u}_i^{n+1} - \bar \mathbf{u}_i^n </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (120)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (120)
 
|}
 
|}
  
−
===Step 5. Check frictional contact forces===
+
====Step 5. Check frictional contact forces====
  
−
===Step 6. Update residual force vector Go to step 1===
+
====Step 6. Update residual force vector====
 +
 
 +
====Go to step 1====
  
 
==7 EXAMPLES==
 
==7 EXAMPLES==
  
 
===7.1 Uniaxial compression and tension tests of a rock material===
 
===7.1 Uniaxial compression and tension tests of a rock material===
 +
 +
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.
 +
  
 
<div id='img-14'></div>
 
<div id='img-14'></div>
Line 1,868: Line 1,849:
 
|}
 
|}
  
−
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.
+
 
 +
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>%.
 +
 
  
 
<div id='img-15'></div>
 
<div id='img-15'></div>
Line 1,878: Line 1,861:
 
|}
 
|}
  
−
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>%
 
  
 
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.
 
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.
Line 1,919: Line 1,901:
  
 
===7.2 Simulation of rock cutting===
 
===7.2 Simulation of rock cutting===
 +
 +
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>.
  
 
<div id='img-20'></div>
 
<div id='img-20'></div>
Line 1,928: Line 1,912:
 
|}
 
|}
  
−
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>.
 
  
 
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.
 
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.
Line 2,008: Line 1,991:
 
===7.3 Strip punch test===
 
===7.3 Strip punch test===
  
−
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>
+
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. 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.
 +
 
 +
<div id='img-28'></div>
 
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
 
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
 
|-
 
|-
Line 2,015: Line 2,000:
 
| colspan="1" | '''Figure 28:''' Strip punch geometry and loading configuration
 
| colspan="1" | '''Figure 28:''' Strip punch geometry and loading configuration
 
|}
 
|}
−
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.
+
 
  
 
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>
 
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>
Line 2,031: Line 2,016:
 
| colspan="1" | '''Figure 30:''' Experimental axial stress vs. axial strain plot, after <span id='citeF-41'></span>[[#cite-41|[41]]]
 
| colspan="1" | '''Figure 30:''' Experimental axial stress vs. axial strain plot, after <span id='citeF-41'></span>[[#cite-41|[41]]]
 
|}
 
|}
 +
 +
 +
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.
  
 
<div id='img-31'></div>
 
<div id='img-31'></div>
Line 2,040: Line 2,028:
 
|}
 
|}
  
−
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.
 
  
 
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>.
 
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>.
Line 2,091: Line 2,078:
 
|[[Image:Draft_Samper_102318540-rock-bara-ps3t-c2-ini.png|258px|]]
 
|[[Image:Draft_Samper_102318540-rock-bara-ps3t-c2-ini.png|258px|]]
 
|[[Image:Draft_Samper_102318540-rock-bara-ps3t-c2-bb-024.png|258px|]]
 
|[[Image:Draft_Samper_102318540-rock-bara-ps3t-c2-bb-024.png|258px|]]
 +
|[[Image:Draft_Samper_102318540-rock-bara-ps3t-c2-bb-040.png|264px|]]
 
|-
 
|-
−
|[[Image:Draft_Samper_102318540-rock-bara-ps3t-c2-bb-040.png|264px|]]
 
 
|[[Image:Draft_Samper_102318540-rock-bara-ps3t-c2-bb-056.png|270px|]]
 
|[[Image:Draft_Samper_102318540-rock-bara-ps3t-c2-bb-056.png|270px|]]
−
|-
 
 
|[[Image:Draft_Samper_102318540-rock-bara-ps3t-c2-bb-080.png|282px|]]
 
|[[Image:Draft_Samper_102318540-rock-bara-ps3t-c2-bb-080.png|282px|]]
 
|[[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]]
 
|[[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]]
 
|- style="text-align: center; font-size: 75%;"
 
|- style="text-align: center; font-size: 75%;"
−
| colspan="2" | '''Figure 35:''' Impact of projectile &#8211; numerical results obtained using linear triangular elements based on FIC formulation
+
| colspan="3" | '''Figure 35:''' Impact of projectile &#8211; numerical results obtained using linear triangular elements based on FIC formulation
 
|}
 
|}
  
Line 2,106: Line 2,092:
 
|[[Image:Draft_Samper_102318540-rock-bara-ps3-c2-ini.png|252px|]]
 
|[[Image:Draft_Samper_102318540-rock-bara-ps3-c2-ini.png|252px|]]
 
|[[Image:Draft_Samper_102318540-rock-bara-ps3-c2-bb-024.png|258px|]]
 
|[[Image:Draft_Samper_102318540-rock-bara-ps3-c2-bb-024.png|258px|]]
 +
|[[Image:Draft_Samper_102318540-rock-bara-ps3-c2-bb-040.png|264px|]]
 
|-
 
|-
−
|[[Image:Draft_Samper_102318540-rock-bara-ps3-c2-bb-040.png|264px|]]
 
 
|[[Image:Draft_Samper_102318540-rock-bara-ps3-c2-bb-056.png|276px|]]
 
|[[Image:Draft_Samper_102318540-rock-bara-ps3-c2-bb-056.png|276px|]]
−
|-
 
 
|[[Image:Draft_Samper_102318540-rock-bara-ps3-c2-bb-080.png|282px|]]
 
|[[Image:Draft_Samper_102318540-rock-bara-ps3-c2-bb-080.png|282px|]]
 
|[[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]]
 
|[[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]]
 
|- style="text-align: center; font-size: 75%;"
 
|- style="text-align: center; font-size: 75%;"
−
| colspan="2" | '''Figure 36:''' Impact of projectile &#8211; numerical results obtained using mixed bilinear Q1/P0 quadrilateral elements
+
| colspan="3" | '''Figure 36:''' Impact of projectile &#8211; numerical results obtained using mixed bilinear Q1/P0 quadrilateral elements
 
|}
 
|}
  
Line 2,128: Line 2,113:
  
 
Interaction of soil and pipe leading to pipe ovalization has been studied. This example demonstrates another possibility of combined FEM/DEM modelling in geomechanics.
 
Interaction of soil and pipe leading to pipe ovalization has been studied. This example demonstrates another possibility of combined FEM/DEM modelling in geomechanics.
 +
 +
<div id='img-38'></div>
 +
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
 +
|-
 +
|[[Image:Draft_Samper_102318540-pipeoval.png|452px|Pipe ovalization &#8211; geometry definition]]
 +
|- style="text-align: center; font-size: 75%;"
 +
| colspan="1" | '''Figure 38:''' Pipe ovalization &#8211; geometry definition
 +
|}
 +
  
 
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
 
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
Line 2,194: Line 2,188:
 
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]]].
 
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]]].
  
−
<div id='img-38'></div>
 
−
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
 
−
|-
 
−
|[[Image:Draft_Samper_102318540-pipeoval.png|492px|Pipe ovalization &#8211; geometry definition]]
 
−
|- style="text-align: center; font-size: 75%;"
 
−
| colspan="1" | '''Figure 38:''' Pipe ovalization &#8211; geometry definition
 
−
|}
 
  
 
<div id='img-39'></div>
 
<div id='img-39'></div>
Line 2,217: Line 2,204:
 
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
 
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
 
|-
 
|-
−
|[[Image:Draft_Samper_102318540-pipeovaln3s-wd-eps-0.png|246px|]]
+
|[[Image:Draft_Samper_102318540-pipeovaln3s-wd-eps-0.png|180px|]]
−
|[[Image:Draft_Samper_102318540-pipeovaln3s-wd-eps-02.png|330px|]]
+
|[[Image:Draft_Samper_102318540-pipeovaln3s-wd-eps-02.png|230px|]]
 
|-
 
|-
−
|[[Image:Draft_Samper_102318540-pipeovaln3s-wd-eps-04.png|348px|]]
+
|[[Image:Draft_Samper_102318540-pipeovaln3s-wd-eps-04.png|270px|]]
−
|[[Image:Draft_Samper_102318540-pipeovaln3s-wd-eps-06.png|366px|Pipe ovalization &#8211; deformed configuration]]
+
|[[Image:Draft_Samper_102318540-pipeovaln3s-wd-eps-06.png|306px|Pipe ovalization &#8211; deformed configuration]]
 
|- style="text-align: center; font-size: 75%;"
 
|- style="text-align: center; font-size: 75%;"
 
| colspan="2" | '''Figure 40:''' Pipe ovalization &#8211; deformed configuration
 
| colspan="2" | '''Figure 40:''' Pipe ovalization &#8211; deformed configuration
Line 2,238: Line 2,225:
 
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]]].
 
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]]].
  
−
===BIBLIOGRAPHY===
+
==References==
  
 
<div id="cite-1"></div>
 
<div id="cite-1"></div>
Line 2,372: Line 2,359:
 
'''[[#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).
 
'''[[#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).
  
−
==8 APPENDIX==
+
==APPENDIX A==
  
−
===8.1 Matrices and vectors in FEM equations of motion using the FIC formulation===
+
===A.1 Matrices and vectors in FEM equations of motion using the FIC formulation===
  
 
In FEM equations of motion using the FIC formulation ([[#eq-94|94]])&#8211;([[#eq-96|96]]) the following vectors and matrices have been used:
 
In FEM equations of motion using the FIC formulation ([[#eq-94|94]])&#8211;([[#eq-96|96]]) the following vectors and matrices have been used:
Line 2,383: Line 2,370:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| style="text-align: center;" | <math>{M}_{ij}=\int _{\Omega ^e} \rho {N}_i {N}_j\,d\Omega \quad ,\quad {N}_i = N_I {I}_3 </math>
+
| 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>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (A.1)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (A.1)
Line 2,396: Line 2,383:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| style="text-align: center;" | <math>{g}= \int _\Omega {B}^T {\boldsymbol \sigma } d\Omega </math>
+
| style="text-align: center;" | <math>{\bf g}= \int_\Omega {\bf B}^T {\boldsymbol\sigma} d\Omega </math>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (A.2)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (A.2)
Line 2,402: Line 2,389:
  
 
is the internal nodal force vector and the rest of matrices and vectors are
 
is the internal nodal force vector and the rest of matrices and vectors are
−
 
  
 
{| class="formulaSCP" style="width: 100%; text-align: left;"  
 
{| class="formulaSCP" style="width: 100%; text-align: left;"  
Line 2,409: Line 2,395:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| style="text-align: center;" | <math>\begin{array}{ll}{\bf G}_{ij} =  \int_{\Omega^e} ({\boldsymbol\nabla} N_i) N_j \, d\Omega & \\
+
| style="text-align: center;" | <math>{\bf G}_{ij} =  \int_{\Omega^e} ({\boldsymbol\nabla} N_i) N_j \, d\Omega a </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\\
+
|-
−
\bar {\bf C}  =  \left[\matrix{ \bar {\bf C}^1 & {\bf 0}&{\bf 0}\cr
+
| 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>
−
  {\bf 0}&\bar {\bf C}^2 & {\bf 0}\cr {\bf 0}&{\bf 0}&\bar {\bf C}^3\cr}\right]\qquad ,\quad & \bar C_{ij}^k = \int_{\Omega^e} \tau_k N_i N_j \, d\Omega \\
+
|-
−
{\bf Q} = [{\bf Q}^1, {\bf Q}^2, {\bf Q}^3]\quad , \quad & Q_{ij}^k =
+
| style="text-align: center;" | <math> {\bf C}  =  \left[\begin{matrix} \bar {\bf C}^1 & {\bf 0}&{\bf 0}\\
−
\int_{\Omega^e} \tau_k \frac{\partial N_i}{\partial x_k} N_j d\Omega\nonumber\\
+
  {\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>
−
{\bf f}_i = \int_{\Omega^e} N_i {\bf b} \, d\Omega +
+
−
\int_\Gamma N_i \bar {\bf t} \, d\Gamma \quad , \quad i,j=1, n_d &\end{array}</math>
+
−
|}
+
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (A.3)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (A.3)
 +
|-
 +
| 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>
 +
|-
 +
| style="text-align: center;" | <math> {\bf f}_i = \int_{\Omega^e} N_i {\bf b} \, d\Omega +
 +
\int_\Gamma N_i \bar {\bf t} \, d\Gamma \quad , \quad i,j=1, n_d </math>
 +
|}
 
|}
 
|}
  
−
In above <math display="inline">{b}=[b_1,b_2,b_3]^T</math> and <math display="inline">\bar {t} = [\bar t_1, \bar t_2,\bar t_3]^T</math>,
+
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>,
  
 
{| class="formulaSCP" style="width: 100%; text-align: left;"  
 
{| class="formulaSCP" style="width: 100%; text-align: left;"  
Line 2,428: Line 2,417:
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
{| style="text-align: left; margin:auto;width: 100%;"  
 
|-
 
|-
−
| 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} &  &{\textbf 0}\\ &  \tau _{2} & \\ {\textbf 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>
+
| 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>
 
|}
 
|}
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (A.4)
 
| style="width: 5px;text-align: right;white-space: nowrap;" | (A.4)
 
|}
 
|}
  
−
'''B''' is the standard infinitesimal strain matrix and <math display="inline">{D}_d</math> is the deviatoric constitutive matrix <span id='citeF-3'></span>[[#cite-3|[3]]]. Note that the expression of <math display="inline">{g}</math> of eq. ([[#eq-2|2]]) is adequate for non linear structural analysis.
+
<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.

Latest revision as of 12:14, 12 November 2018

Published in Comput. Methods Appl. Mech. Engrg., Vol 193, 3087-3128, 2004
doi: 10.1016/j.cma.2003.12.056


Abstract

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.

1 Introduction

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 [1,2].

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

In this paper a stabilization method based on the finite calculus (FIC) approach [4,5] is used for deriving locking-free linear triangles and tetrahedra which can be effectively used in conjunction with the DEM.

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.

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.

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.

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.

2 DISCRETE ELEMENT METHOD FORMULATION

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 [6,7,8,9]. Formulation of spherical discrete elements (called also distinct elements) following the main assumptions of Cundall [6,10] has been developed by Rojek et al. in [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.

2.1 Equations of motion

The translational and rotational motion of rigid spherical or cylindrical elements (particles) (Fig. 1) is described by means of the standard equations of rigid body dynamics. For the Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle i} -th element we have

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): m_i \ddot{\textbf{u}}_i= {\mathbf{F}}_i\,,
(1)

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): I_i \dot{{\omega }}_i= {\mathbf{T}}_i\,,
(2)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {u} is the element centroid displacement in a fixed (inertial) coordinate frame Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \mathbf{X} , Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \omega – the angular velocity, Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle m} 

– the element mass, Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle I}
– the moment of inertia, Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \mathbf{F}

– the resultant force, and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \mathbf{T} – the resultant moment about the central axes. Vectors Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \mathbf{F} 

and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \mathbf{T}
are sums of all forces and moments applied to the Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle i}

-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 (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 Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \mathbf{x} 

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 Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \mathbf{X}

, and hence the rotational motion can be easily described in this system.

Motion of a rigid particle
Figure 1: Motion of a rigid particle

Equations of motion (1) and (2) are integrated in time using a central difference scheme. The time integration operator for the translational motion at the Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle n} -th time step is as follows:

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \ddot{{\textbf u}}_i^{n}=\frac{{{\textbf F}}_i^n}{m_i}\,,
(3)

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \dot{{\textbf u}}_i^{n+1/2}=\dot{{\textbf u}}_i^{n-1/2}+\ddot{{\textbf u}}_i^{n}{\Delta t}\,,
(4)

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {{\textbf u}}_i^{n+1}={{\textbf u}}_i^{n}+\dot{{\textbf u}}_i^{n+1/2}{\Delta t}\,.
(5)

The first two steps in the integration scheme for the rotational motion are identical to those given by Eqs.(3) and (4):

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \dot{{\omega }}_i^{n}=\frac{{{T}}_i^n}{I_i}\,,
(6)

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {{\omega }}_i^{n+1/2}={{\omega }}_i^{n-1/2}+\dot{{\omega }}_i^{n}{\Delta t}\,.
(7)

The vector of incremental rotation Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\Delta \theta }=\{ \Delta \theta _x\,\Delta \theta _y\,\Delta \theta _z\} ^{\mathrm{T}}} 

is calculated as

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {{\Delta \theta }}_i={{\omega }}_i^{n+1/2}{\Delta t}\,,
(8)

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

2.2 Contact search algorithm

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 Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle n^2} , where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle n} 

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 Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle n\ln{n}}

, which allows to solve large frictional contact systems.

2.3 Evaluation of contact forces

2.3.1 Decomposition of the contact force

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 Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\textbf F}_1} 

and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\textbf F}_2}

, which by the Newton's third law satisfy the following relation:

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textbf F}_{1}=-{\textbf F}_2\,.
(9)

We take Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\textbf F}={\textbf F}_1} 

and decompose Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textbf F}

into the normal and tangential components, Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {F}_n} and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\textbf F}_T} , respectively (Fig. 2)

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textbf F}={\textbf F}_n+{\textbf F}_T =F_n{\textbf n}+{\textbf F}_T \,,

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textbf n} 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.

The contact forces Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle F_n} 

and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\textbf F}_T}
are obtained using a constitutive model formulated for the contact between two rigid spheres (or discs in 2D) (Fig. 3). The contact interface in our formulation  is characterized by the normal and tangential stiffness Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle k_n}
and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle k_T}

, respectively, the Coulomb friction coefficient Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \mu } , and the contact damping coefficient Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle c_n} .

Decomposition of the contact force into the normal and tangential components
Figure 2: Decomposition of the contact force into the normal and tangential components
Model of the contact interface
Figure 3: Model of the contact interface

2.3.2 Normal contact force

The normal contact force Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle F_n} 

is decomposed into the elastic part Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle F_{ne}}
and the damping contact force Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle F_{nd}}


Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): F_n = F_{ne} + F_{nd}\,.
(10)

The damping part is used to decrease oscillations of the contact forces and to dissipate kinetic energy.

The elastic part of the normal contact force Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle F_{ne}} 

is proportional to the normal stiffness Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle k_n}
and to the penetration of the two particle surfaces Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle u_{rn}}

, i.e.

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): F_{ne}=k_n u_{rn}\,.
(11)

The penetration is calculated as

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): u_{rn} = d - r_1 - r_2\,,
(12)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle d} 

is the distance between the particle centres, and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle r_1}

, Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle r_2} 

their radiae. In the  present study no cohesion is allowed, so  no tensile normal contact forces are allowed and hence

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): F_{ne}\le 0\,.
(13)

If Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle u_{rn}<0} , Eq. (11) holds, otherwise Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle F_{ne}= 0} .

The contact damping force is assumed to be of viscous type and given by

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): F_{nd}=c_n v_{rn}
(14)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle v_{rn}} 

is the normal relative velocity of the centres of the two particles in contact defined by

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {v}_{rn}=(\dot{{\textbf u}}_2-\dot{{\textbf u}}_1)\cdot {\textbf n}\,.
(15)

The damping coefficient Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle c_n} 

can be taken as a fraction of the critical damping Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle C_{cr}}
for the system of two rigid bodies with masses Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle m_1}
and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle m_2}

, connected with a spring of stiffness Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle k_n} 

([13]) with

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): C_{cr}=2\sqrt{\frac{m_1m_2k_n}{m_1+m_2} }\,.
(16)

2.3.3 Tangential frictional contact

In the absence of cohesion (if the particles were not bonded at all or the initial cohesive bond has been broken) the tangential reaction Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\textbf F}_T} appears by friction opposing the relative motion at the contact point. The relative tangential velocity at the contact point Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\textbf v}_{rT}} 

is calculated from the following relationship

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textbf v}_{rT}= {\textbf v}_{r}-{\textbf v}_{r}\cdot {\textbf n}\,,
(17)

with

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textbf v}_r=(\dot{{\textbf u}}_{2} +{\omega }_{2} \times {\textbf r}_{c2})- (\dot{{\textbf u}}_{1} +{\omega }_1 \times {\textbf r}_{c1})\,,
(18)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \dot{{\textbf u}}_{1}} , Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \dot{{\textbf u}}_{2}} , and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\omega }_1} , Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\omega }_2} 

are the translational and rotational velocities of the particles, and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\textbf  r}_{c1}}
and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\textbf r}_{c2}}
are the vectors connecting particle centres with contact points.
Draft Samper 102318540-coulomb.png Friction force vs. relative tangential displacement a) Coulomb law, b) regularized Coulomb law
Figure 4: Friction force vs. relative tangential displacement a) Coulomb law, b) regularized Coulomb law

The relationship between the friction force Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \| {\textbf F}_T\| } and the relative tangential displacement Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle u_{rT}} 

for the classical Coulomb model (for a constant normal force Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle F_n}

) is shown in Fig. 4a. 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, Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\textbf v}_{rT}^{\mathrm{r}}} 

and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\textbf v}_{rT}^{\mathrm{ir}}}

, respectively as:

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textbf v}_{rT}= {\textbf v}_{rT}^{\mathrm{r}}+{\textbf v}_{r}^{\mathrm{ir}}\,.
(19)

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. 4b. This analogy allows us to calculate the friction force employing the standard radial return algorithm []. First a trial state is calculated

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textbf F}_T^{\mathrm{trial}}={\textbf F}_T^{\mathrm{old}}-k_T {\textbf v}_{rT} \Delta t\,,
(20)

and then the slip condition is checked

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \phi ^{\mathrm{trial}}=\| {\textbf F}_T^{\mathrm{trial}}\|{-\mu}|F_n|\,.
(21)

If Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \phi ^{\mathrm{trial}}\le{0}} , we have stick contact and the friction force is assigned the trial value

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textbf F}_T^{\mathrm{new}}={\textbf F}_T^{\mathrm{trial}}\,,
(22)

otherwise (slip contact) a return mapping is performed giving

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textbf F}_T^{\mathrm{new}}=\mu |F_n|\frac{{\textbf F}_T^{\mathrm{trial}}}{\| {\textbf F}_T^{\mathrm{trial}} \| }\,.
(23)

2.4 Background damping

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 Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {{F}}_i^{\mathrm{damp}}} 

and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\mathbf{T}}_{i}^{\mathrm{damp}}}
are  added to equations of motion (1) and (2)

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): m_i \ddot{\mathbf{u}}_i= {\mathbf{F}}_i+{\mathbf{F}}_i^{\mathrm{damp}}\,,
(24)

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): I_i \dot{{\omega }}_i= {\mathbf{T}}_i+{\mathbf{T}}_i^{\mathrm{damp}}\,.
(25)

where

  • for viscous damping

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\mathbf{F}}_{i}^{\mathrm{damp}}=-\alpha ^{vt}m_{i}\dot{\mathbf{ u}}_{i}\,,
(26)

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\mathbf{T}}_{i}^{\mathrm{damp}}=-\alpha ^{vr}I_{i}{\omega }_{i}\,,
(27)
  • for non-viscous damping

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\mathbf{F}}_{i}^{\mathrm{damp}}=-\alpha ^{nvt}\Vert {\mathbf{F}}_{i}\Vert \frac{\dot{\mathbf{u}}_{i}}{\Vert \dot{\mathbf{u}}_{i}\Vert }\,,
(28)

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\mathbf{T}}_{i}^{\mathrm{damp}}=-\alpha ^{nvr}\Vert {\mathbf{T}}_{i}\Vert \frac{{\omega }_{i}}{\Vert {\omega }_{i}\Vert }\,.
(29)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \alpha ^{vt}} , Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \alpha ^{vr}} , Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \alpha ^{nvt}} 

and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \alpha ^{nvr}}
are damping constants. It can be seen from Eqs. (26)–(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.

2.5 Numerical stability

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 Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\Delta t}} , i.e.

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\Delta t} \le {\Delta t}_{cr}
(30)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\Delta t}_{cr} 

is a critical time step determined by the highest natural frequency of the system Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \omega _{ max}


Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\Delta t}_{cr}=\frac{2}{\omega _{max}}\,.
(31)

If damping exists, the critical time increment is given by

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\Delta t}_{cr}=\frac{2}{\omega _{max}}\left(\sqrt{1+\xi ^2}-\xi \right)\,,
(32)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \xi } 

is the fraction of the critical damping corresponding to the highest frequency Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \omega _{max}

. Exact calculation of the highest frequency Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \omega _{max} 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

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\mathbf{m}}_i \ddot{\mathbf{r}}_i+ {\mathbf{k}}_i {\mathbf{r}}_i= \mathbf{0}\,,
(33)

where

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\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}}\,,
(34)

and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\mathbf{k}}_i} 

is the stiffness matrix accounting for the contributions from all the penalty constraints active for the Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle i}

-th particle. Equation (34) defines vectors Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\mathbf{m}}_i} 

and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\mathbf{r}}_i}
for a spherical particle in a 3D space. For a cylindrical particle in a 2D  model the respective vectors are defined as

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\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}}\,.
(35)

Equation (33) leads to the eigenproblem

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\mathbf{k}}_i {\mathbf{r}}_i= \lambda _j {\mathbf{m}}_i {\mathbf{r}}_i\,,
(36)

whose eigenvalues Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \lambda _j} 

(Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle j=1,\ldots , 6}
in 3D case, and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle j=1,2,3}
for 2D case) are the squared frequencies of the free vibration frequencies:

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \lambda _j=\omega _j^2\,.
(37)

In a 3D problem, three of six frequencies Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \omega _j} 

are translational, and the other three – rotational.

In the algorithm implemented, a further simplification is assumed. The maximum frequency is estimated as the maximum of natural frequencies of mass–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:

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): m_i \ddot{u}_n + k_n u_n=0\,,
(38)

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): I_i \ddot{\theta } + k_\theta \theta=0\,,
(39)

where it is assumed that the translational motion Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle u_n} is due to the contact interaction along the normal direction (the spring stiffness Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle k_n} 

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 Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle k_T}

, the rotational stiffness Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle k_\theta } 

can be obtained as

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): k_{\theta }=k_T r^2\,,
(40)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle r} 

is the length of the vector connecting the mass centre to the contact point.

The natural frequency of the translational vibrations is given by:

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \omega _n=\sqrt{\frac{k_n}{m_i}}\,,
(41)

while the rotational frequency Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \omega _{\theta }} 

can be obtained from

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \omega _{\theta }=\sqrt{\frac{k_{\theta }}{I_i}}\,.
(42)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle I_i} 

is the rotational inertia of a sphere
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): I= \frac{2}{5} mr^2
(43)

and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle k_\theta } 

is given by Eq. (40). Substituting these expressions into Eq. (42) gives  the rotational frequency  as
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \omega _{\theta }=\sqrt{\frac{5 k_{T}}{2 m_i}}\,.
(44)

If Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle k_T=k_n} , the rotational frequency Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \omega _{\theta }} 

is considerably higher than the translational frequency Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \omega _n}
obtained from Eq. (41), which results in a smaller critical time increment (see Eq. (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 [14].

3 MICROMECHANICAL CONSTITUTIVE MODELS

3.1 Elastic perfectly brittle model

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:

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \sigma \!=\!k_n u_n\,, (45)
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \tau \!=\!k_t \, u_t\,, (46)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \sigma } 

and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \tau }
are the  normal and tangential contact force, respectively, Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle k_n}
and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle k_t}
are the interface stiffness in the normal and tangential directions and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle u_n}
and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle u_t}
the normal and tangential relative displacements, respectively.

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:

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \sigma \!\le \!R_n\,, (47)
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): |\tau | \!\le \!R_t\,, (48)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle R_n} 

and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle R_t}
are the interface strengths in the normal and tangential directions, respectively.

In the absence of cohesion the normal contact force can be compressive only, i.e.

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \sigma \le 0
(49)

and the (positive) tangential contact force is given by

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \tau = \mu |\sigma |
(50)

if Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \sigma{<0}} 

 or zero otherwise. The friction force is given by Eq. (50) expressing the Coulomb friction law, with Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \mu }
 being the Coulomb friction coefficient.
Normal contact force in the elastic perfectly brittle model
Figure 5: Normal contact force in the elastic perfectly brittle model
Tangential contact force in the elastic perfectly brittle model (tensile normal contact force)
Figure 6: Tangential contact force in the elastic perfectly brittle model (tensile normal contact force)
Failure surface for the elastic perfectly brittle model
Figure 7: Failure surface for the elastic perfectly brittle model

Contact laws for the normal and tangential directions for the elastic perfectly brittle model are shown in Figs. 5 and 6, respectively. The failure surface for the elastic perfectly brittle model defined by conditions (47) and (48) is shown in Fig. 7.

3.2 Elasto-plastic contact model with linear softening

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:

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \sigma _i=k_i(u_i-u_i^p)\,, \quad i=n,t
(51)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \sigma _i} 

is the normal or tangential contact force, Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle u_i}
the total relative displacement, Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle u_i^p}
the plastic component of the relative displacement, Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle k_i}
the elastic parameters which characterize the contact interface stiffness and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle n, t}
the  normal and tangential directions.

A general yield criterion has been assumed in the following form:

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): f(\sigma ,\alpha )=|\sigma | -[\sigma _Y -H \alpha ]\,,
(52)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle H} 

 is  the plastic softening modulus (assumed to be positive) and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \alpha }
 is  the softening  parameter represented by the plastic part of the relative displacement Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle u^p}

.

An elasto-plastic contact law with linear softening is assumed for the shear and normal tensile contact forces (Figs. 8 and 9). An elastic linear law is assumed for compression.

Dependence of the tangential force on the normal force is expressed by the application of the 2D failure criterion shown in Fig. 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.

Flow rule

The plastic strain increments for both formulations (i.e. the normal and tangential motions) are calculated from the following flow rule:

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \dot{u}^p_i=\gamma \, \mbox{sign}(\sigma _i) \,,
(53)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \gamma } 

is the absolute value of the slip rate (Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \gamma{>0}}

), Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \dot{u}^p_i} 

the derivative of the plastic relative displacement (normal or tangential) and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \sigma _i}
the contact force in the normal or tangential direction.
Elasto-plastic contact law with softening for the normal direction
Figure 8: Elasto-plastic contact law with softening for the normal direction
Elasto-plastic contact law with softening for the tangential direction
Figure 9: Elasto-plastic contact law with softening for the tangential direction
Draft Samper 102318540 6786 Fig10.png
Figure 10: Yield surface for elasto-plastic model

3.3 Contact model with elastic damage

The elastic perfectly brittle model can be easily extended to account for elastic damage by assuming a softening behaviour defined by a softening modulus Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle H}

introduced into the force-displacement relationship (Fig. 11).
Normal contact force in the contact model with elastic damage
Figure 11: Normal contact force in the contact model with elastic damage

The constitutive relationship for 1D elastic damage is given by:

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \sigma =k_n^D \,u_n =(1-\omega )k_n u_n
(54)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle k_n^D } 

is the elastic damaged secant modulus and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \omega }
the scalar  damage variable.

The scalar damage variable Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \omega } 

is a measure of material damage. For the undamaged state Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \omega=0}
and for a damaged state Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle 0<\omega \le 1}

. The scalar damage variable Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \omega } 

can be written in the following form:

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \omega = \frac{\psi (u_n)-1}{\psi (u_n)}
(55)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \psi (u_n)} 

 is a function of total relative displacement. For linear strain softening Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \psi (u_n)}
  is defined by

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \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.
(56)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle R_n} 

is the initial tensile strength and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle H}
is the softening modulus (taken as positive).

A similar contact force–displacement law with damage (Fig. 12) can be introduced for the tangential direction. Then the bond decohesion can occur either due to tension or shear.

Tangential contact force in the contact model with elastic damage
Figure 12: Tangential contact force in the contact model with elastic damage

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.

After the contact bonds are broken due to damage (Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \omega=0} ) standard frictional contact can occur between the spherical elements.

3.4 Contact model with friction, wear and heat generation

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 Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \sigma } 

 is calculated from the linear relationships:

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \sigma = k_n u_n
(57)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle u_n\le 0} 

 is the normal relative displacement (penetration of one sphere against another). The frictional contact force is evaluated using the regularized Coulomb law

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \tau = \mu |\sigma |
(58)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \mu } 

 is the Coulomb friction coefficient. The frictional dissipation rate Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \dot{D}}
is calculated as

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \dot{D} = \tau v_t
(59)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle v_t} 

 is the relative tangential velocity. The frictional dissipation is used to calculate the wear rate Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \dot{w}}
and the heat generated by friction Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle Q_{gen}}

. This is computed by

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): Q_{gen} = \chi \dot{D}
(60)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \chi } 

is the part of the friction work converted to heat (Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \chi \le 1}

). Heat generation through frictional dissipation is assumed to be absorbed equally by the two particles in contact.

Wear rate is calculated using the Archard law [15], which assumes that the wear rate Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \dot{w}} 

is proportional to the contact pressure Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \sigma }
  and to the slip velocity Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle v_T}


Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \dot{w}=k \frac{\sigma v_T}{H}
(61)

with Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle H} 

being the hardness of worn surface and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle k}
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 [16]. Values of adhesive and abrasive wear constants Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle k}
 for different combinations of materials can be found in  [16]. It is commonly accepted that wear is related to friction and thus friction coefficients are often introduced into Eq. (61). If the Coulomb friction law (58) holds,  the Archard wear law can be written in the following equivalent form:

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \dot{w}=\bar{k} \frac{\tau v_T}{H}= \bar{k} \frac{\dot{D}}{H}
(62)

in which the wear rate is related to frictional dissipation rate Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \dot{D}} 

  and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \bar{k}=k/\mu }

.

Influence of temperature on wear can be taken into account by adaptation of the law of Archard given by Eq. (61). Thermal effects can be captured approximately by taking the hardness Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle H} 

as a function of the temperature Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle T}


Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): H=H(T)
(63)

Then the Archard law (61) takes the following form:

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \dot{w}=k \frac{\sigma v_T}{H(T)}
(64)

Computation of temperature can be performed by solving thermal conduction problem on the discrete domain [17].

4 FORMULATION OF QUASI-INCOMPRESSIBLE SOLID FINITE ELEMENTS

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

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 [3,18,19] or alternatively the mixed patch test [3,20,21] not being passed.

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 [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 [3]. Zienkiewicz et al. [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 [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. [24]. Other methods to overcome volumetric locking are based on mixed displacement (or velocity)-pressure formulations using the Galerkin-Least-Square (GLS) method [25], average nodal pressure and average nodal deformation [26,27] techniques, and Sub-Grid Scale (SGS) methods [28,29,30,31].

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.

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 [32,33,34,35,36,37,38,39]. The same ideas were applied in [4,5] to derive a stabilized formulation for quasi-incompressible and incompressible solids allowing the use of linear triangles and tetrahedra.

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.

4.1 Equilibrium equations

The equilibrium equations for an elastic solid are written using the FIC technique as [32]

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {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
(65)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle n_d} 

is the number of space dimensions of the problems (i.e. Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle n_d =3}
for 3D)

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {r}_{i}: ={\partial {\sigma }_{ij} \over \partial {x}_j}+{b}_i
(66)

In (65) and (66) Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\sigma }_{ij} and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {b}_i are the stresses and the body forces, respectively and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {h}_k are characteristic length distances of an arbitrary prismatic domain where equilibrium of forces is considered.

Equations (65) and (66) are completed with the boundary conditions on the displacements Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle u_i}


Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {u}_i-\bar{u}_i=0\quad \mbox{on} \;\,\Gamma _u
(67)

and the equilibrium of surface tractions

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\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
(68)

In the above Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \bar{u}_i and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \bar{t}_i are prescribed displacements and tractions over the boundaries Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Gamma _u} and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Gamma _t} , respectively, Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle n_i} 

are the components of the unit normal vector and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle h_k}
are again the characteristic lengths.

The underlined terms in Eqs. (65) and (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 [32].

These terms are essential to derive a stabilized finite element formulation as described below.

4.1.1 Constitutive equations

The stresses are split into deviatoric and volumetric (pressure) parts

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\sigma }_{ij}= s_{ij} +p{\delta }_{ij}
(69)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\delta }_{ij} is the Kronnecker delta function. The linear elastic constitutive equations for the deviatoric stresses Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle s_{ij}} 

are written as

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): s_{ij}\!\!=\!\! 2G \left({\varepsilon }_{ij}-\frac{1}{3}{\varepsilon }_{v}{\delta }_{ij}\right)
(70)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle G} is the shear modulus,

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\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}\,.
(71)

The constitutive equation for the pressure Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle p} 

can be written  for an arbitrary domain of finite size using the FIC formulation as [30]

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \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
(72)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle K} 

is the bulk modulus of the material.

Note that for Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle h_k \to 0} 

the standard relationship between the pressure and the volumetric strain of the infinitesimal theory Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle (p=K\varepsilon _v )}
is found.

For an incompressible material Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle K\to \infty } 

and Eq. (72) yields

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \varepsilon _v - \frac{{h}_k}{2}{\partial {\varepsilon }_{v} \over \partial {x}_k}=0
(73)

Eq. (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 [32,33].

By combining Eqs. (65), (66), (69), (70) and (73) a mixed displacement–pressure formulation can be written as

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\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
(74)

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \left(\frac{p}{K}-{\varepsilon }_{v}\right)-\frac{{h}_k}{2}{\partial \over \partial {x}_k}\left(\frac{p}{K}-{\varepsilon }_{v}\right)=0
(75)

Substituting Eq. (70) into (74) leads after some algebra to

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\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]
(76)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {r}_{i} is defined by Eq. (66) and

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \hat{r}_{i}=\frac{\partial }{\partial x_j}(2G \varepsilon _{ij})+{\partial p \over \partial x_i} + b_i
(77)

Substituting Eq. (76) into (75) gives after some simplification

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \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
(78)

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {p\over K} - {\varepsilon }_{v}- \sum \limits _{i=1}^{n_d} \tau _i {\partial r_i \over \partial x_i}=0
(79)

where

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \tau _i = {3h_i^2\over 8G}
(80)

The coefficients Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \tau _{i} } 

in Eq. (79) are also referred to as intrinsic time parameters per unit mass  (their dimensions are Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {t^2 m^3\over Kg}}
where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle t}
is the time).

4.2 Non linear transient dynamic analysis

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.

The transient FIC equilibrium equations can be written in an identical form to eq. (65) (neglecting time stabilization terms [32,37,38]) with

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): r_i:= - \rho \frac{\partial ^2{u_i}}{\partial{t}^2} + {\partial \sigma _{ij} \over \partial x_j}+ b_i
(81)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \rho } 

is the  density.

Eq. (81) is completed with the constitutive equations for the deviatoric stresses (eq. (70)) and the pressure (eq. (72)), as well with the boundary conditions (67) and (68) and the initial conditions for Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle t=0} .

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. (79). This equation is now written in an incremental form more suitable for non linear transient analysis.

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 Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle h_i} 

in Eq. (74). The stabilized FIC form of the pressure constitutive equation (Eq. (79)) is however mandatory in order to avoid volumetric locking.

The set of stabilized equations to be solved are now:

Equilibrium

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): r_i =0
(82)

Pressure constitutive equation

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\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
(83)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Delta p = p^{n+1}-p^n} 

and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Delta u_i = u_i^{n+1}-u_i^n}
are the increments of pressure and displacements, respectively. As usual Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle (\cdot )^n}
denotes values at time Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle t_n}

.

In the derivation of Eq. (83) we have accepted that Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Delta r_i =r_i^{n+1} \equiv r_i} 

as the infinitesimal equilibrium equations are assumed to be satisfied at time Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle t_n}
(and hence Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle r_i^n =0}

).

The weighted residual form of the FIC governing equations (82), (68) and (83) is

Equilibrium'

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \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
(84)

Pressure constitutive equation

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \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
(85)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \delta u_i } 

and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle q}
are arbitrary test functions representing virtual displacements and virtual pressure fields, respectively.

Integrating by parts, the terms involving Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle s_{ij}} 

and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle p}
in Eq. (84) and the term involving Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle r_i}
in Eq. (85) and neglecting the space derivatives of the intrinsic time parameters leads to

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \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
(86)

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \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
(87)

The residual Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle r_i} 

is split now as
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): r_i = \pi _i + {\partial p \over \partial x_i}
(88)

where

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \pi _i = - \frac{\partial ^2{u_i}}{\partial{t}^2} + {\partial s_{ij} \over \partial x_i}+b_i
(89)

Note that Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \pi _i} 

is the part of Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle r_i}
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 Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \pi _i}
can be considered belonging to a sub-scale space orthogonal to that of the pressure gradient terms.

In the infinitesimal limit Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle r_i=0} 

and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\partial p \over \partial x_i} + \pi _i=0}

. This limit relationship between Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\partial p \over \partial x_i}} 

and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \pi _i}
can be weakly enforced by means of a weighted residual form.

The final set of integral equations is therefore

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \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
(90)

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \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
(91)

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \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
(92)

where the Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \tau _i} 

coefficients are introduced in Eq. (92) for convenience.

We will choose Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle C^0} continuous linear interpolations of the displacements, the pressure and the pressure gradient projection Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \pi _i} 

over three-node triangles (2D) and four-node tetrahedra (3D) [3]. The linear interpolations are  written as

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {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
(93)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle n=3(4)} 

for 2D(3D) problems and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle (\bar{\cdot })}
denotes nodal variables. As usual Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle N_j}
are the linear shape functions [3]. The nodal variables are a function of the time Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle t}

. Substituting the approximations (93) into eqs. (90)–(92) gives the following system of discretized equations

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \mathbf{M} \ddot{\bar \mathbf{u}} + \mathbf{g} - \mathbf{f}=\mathbf{0}
(94)

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \mathbf{G}^T \Delta \bar \mathbf{u} -\mathbf{C} \Delta \bar \mathbf{p}- \mathbf{L} \bar \mathbf{p} - \mathbf{Q}\bar {\boldsymbol \pi }= {0}
(95)

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \mathbf{Q}^T \bar \mathbf{p}+ \bar \mathbf{C} \bar {\boldsymbol \pi }= \mathbf{0}
(96)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \ddot{\bar \mathbf{u}}} 

is the nodal acceleration vector. The expression of the different matrices in vectors appearing in Eqs. (94)–(96) is presented in the Appendix.

A four steps semi-implicit time integration algorithm can be derived from eqs. (94)–(96) as follows

step

Step 1. Compute the nodal velocities Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \dot{\bar \mathbf{u}}^{n+1/2}}


Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \dot{\bar \mathbf{u}}^{n+1/2} = \dot{\bar \mathbf{u}}^{n-1/2}+\Delta t \mathbf{M}^{-1}_d (\mathbf{f}^n - \mathbf{g}^n)
(97)

Step 2. Compute the nodal displacements Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \bar \mathbf{u}^{n+1}}


Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \bar \mathbf{u}^{n+1} = \bar\mathbf{u}^n + \Delta t \dot{\bar \mathbf{u}}^{n+1/2}
(98)

Step 3. Compute the nodal pressures Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \bar \mathbf{p}^{n+1}}


Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \bar \mathbf{p}^{n+1}=[ \mathbf{C} + \mathbf{L}]^{-1} [\Delta t \mathbf{G}^T \dot{\bar \mathbf{u}}^{n+1/2}+ \mathbf{C} \bar \mathbf{p}^n- \mathbf{Q} \bar {\boldsymbol \pi }^n]
(99)

Step 4. Compute the nodal projected pressure gradients Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \bar {\boldsymbol \pi }^{n+1}}


Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \bar {\boldsymbol \pi }^{n+1}=- \bar \mathbf{C}_d^{-1} \mathbf{Q}^T \bar \mathbf{p}^{n+1}
(100)

In above, all matrices are evaluated at Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle t^{n+1}} , Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle (\cdot )_d =\hbox{diag}\, (\cdot )} 

and
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \mathbf{g}^n =\int _{\Omega ^e} [\mathbf{B}^T {\boldsymbol \sigma }]^n\,d\Omega
(101)

where the stresses Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\boldsymbol \sigma }^n} 

are obtained by consistent integration of the adequate (non linear)  constitutive law [33].

Note that steps 1, 2 and 4 are fully explicit as a diagonal form of matrices Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \mathbf{C}} 

and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \bar\mathbf{C}}
has been chosen. The solution of step 3 with  a diagonal form for Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \mathbf{C}}
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).

A three steps approach can be obtained by evaluating the projected pressure gradient variables Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \bar {\boldsymbol \pi }^{n+1}} 

 at Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle t_{n+1}}
in a fully implicit form in Eq. (99). Eliminating then Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \bar {\boldsymbol \pi }^{n+1}}
 from the fourth step using Eq. (100) and substituting this expression into Eq. (99) leads to
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \bar\mathbf{p}^{n+1}= [\mathbf{C} + \mathbf{L}-\mathbf{S}]^{-1} [\Delta t \mathbf{G}^T \dot{\bar \mathbf{u}}^{n+1/2} + \mathbf{C} \bar\mathbf{p}^n ]
(102)

where

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \mathbf{S} = \mathbf{Q} \bar \mathbf{C}_d^{-1} \mathbf{Q}^T
(103)

Recall that for the full incompressible case Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle K=\infty } 

and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \mathbf{C}=0 }
in all above equations.

The critical time step Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \Delta t} 

is taken as that of the standard explicit dynamic scheme [28].

Explicit algorithm

A fully explicit algorithm can be obtained by computing Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \bar \mathbf{p}^{n+1}} 

from step 3 in Eq. (99)  as follows

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \bar\mathbf{p}^{n+1}=\mathbf{C}_d^{-1}[\Delta t \mathbf{G}^T \dot{\bar \mathbf{u}}^{n+1/2} + (\mathbf{C}_d+\mathbf{L}) \bar\mathbf{p}^n - \mathbf{Q} {\boldsymbol \sigma }^n]
(104)

Obviously, solution of eq. (104) breaks down for Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle K=\infty } 

as Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {C}=0 }
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.

5 DEM–FEM CONTACT ALGORITHM

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. 13. Similarly as in case of contact between two spheres (Fig. 2) the contact force between the sphere and external edge Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \mathbf{F}} 

 of a finite element is decomposed into normal and tangential components, Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \mathbf{F}_n}
and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \mathbf{ F}_T}

. In a general case the contact model between discrete elements and finite element edges can include cohesion, damping, friction, wear, heat generation and exchange.

Contact between a sphere and finite element edge
Figure 13: Contact between a sphere and finite element edge

Similarly to the relationship (10) the normal contact force Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle F_n} 

is decomposed into the elastic part Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle F_{ne}}
and the damping contact force Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle F_{nd}}


Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): F_n = F_{ne} + F_{nd}\,.
(105)

The elastic part of the normal contact force Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle F_{ne}} 

is proportional to the normal stiffness Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle k_n}
and to the penetration of the two particle surfaces Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle u_{rn}}


Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): F_{ne}=k_n u_{rn}\,.
(106)

The penetration is calculated here as

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): u_{rn} = d - r\,,
(107)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle d} 

is the distance of the particle centre Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle C}
and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle r}
is  the particle radius.

The contact damping force is assumed to be of viscous type

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): F_{nd}=c_n v_{rn}
(108)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle c_n} 

is the damping coefficient and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle v_{rn}}
is the normal component of the relative velocity at the contact point. The relative velocity at the contact point is calculated as
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \mathbf{v}_{r}= ({\mathbf{v}}_{C} +{\omega }_{C} \times \mathbf{r})- ({\mathbf{v}}_{1}H_{1} + {\mathbf{v}}_{2}H_{2}) \,,
(109)

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\mathbf{v}}_{C} +{\omega }_{C} \times \mathbf{r}} 

is the discrete element velocity at the contact point and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\mathbf{v}}_{1}H_{1} + {\mathbf{v}}_{2}H_{2}}
is the velocity of the finite element edge at the contact point expressed in terms of the nodal velocities, Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\mathbf{v}}_{1}}
and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\mathbf{v}}_{2}}

, and respective shape functions Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle H_{1}} 

and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle H_{2}}

. The relative velocity is decomposed into the normal and tangential velocity using the formulae (15) and (17).

The cohesive contact force is calculated using the elastic-perfectly brittle model described in Section 3.1. In the absence of cohesion the frictional contact reaction is calculated using the regularized Coulomb friction model described in section 2.3.3.

6 GENERAL DISCRETE ELEMENT-FINITE ELEMENT DYNAMIC FORMULATION

The general algorithm for the transient dynamic problem involving discrete elements and (stabilized) finite elements includes the following steps.

Step1. Compute the nodal velocities

Discrete elements

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \dot{\mathbf{u}}_i^{n+1/2}=\dot{\mathbf{u}}_i^{n-1/2}+\ddot{\mathbf{u}}_i^{n}{\Delta t}\,,
(110)

with

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \ddot{\mathbf{u}}_i^{n}=\frac{{\mathbf{F}}_i^n}{m_i}\,,
(111)

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {{\omega }}_i^{n+1/2}={{\omega }}_i^{n-1/2}+\dot{{\omega }}_i^{n}{\Delta t}\,.
(112)

with

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \dot{{\omega }}_i^{n}=\frac{{\mathbf{T}}_i^n}{I_i}\,,
(113)

Finite elements

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \dot{\bar \mathbf{u}}^{n+1/2} = \dot{\bar \mathbf{u}}^{n-1/2}+\Delta t \mathbf{M}^{-1}_d (\mathbf{f}^n - \mathbf{g}^n)
(114)

Step2. Compute the nodal displacements

Discrete elements

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\mathbf{u}}_i^{n+1}={\mathbf{u}}_i^{n}+\dot{\mathbf{u}}_i^{n+1/2}{\Delta t}\,.
(115)

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {{\Delta \theta }}_i={{\omega }}_i^{n+1/2}{\Delta t}\,,
(116)

Finite elements

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \bar \mathbf{u}^{n+1} = \bar\mathbf{u}^n + \Delta t \dot{\bar \mathbf{u}}^{n+1/2}
(117)

Step3. Compute the nodal pressures

Finite elements

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \bar\mathbf{p}^{n+1}= [\mathbf{C} + \mathbf{L}-\mathbf{S}]^{-1} [\Delta t \mathbf{G}^T \dot{\bar \mathbf{u}}^{n+1/2} + \mathbf{C} \bar\mathbf{p}^n ]
(118)

Step4. Update the nodal coordinates

Discrete and finite elements

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \mathbf{X}_i^{n+1} = \mathbf{X}_i^{n}+\Delta \bar \mathbf{u}_i
(119)
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \Delta \bar \mathbf{u}_i =\bar \mathbf{u}_i^{n+1} - \bar \mathbf{u}_i^n
(120)

Step 5. Check frictional contact forces

Step 6. Update residual force vector

Go to step 1

7 EXAMPLES

7.1 Uniaxial compression and tension tests of a rock material

Mechanical properties of rock materials are determined from laboratory tests such as compression, triaxial and tension tests. Figure 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.


Rock specimen after laboratory unconfined compression test
Figure 14: Rock specimen after laboratory unconfined compression test


Figure 15 presents the material sample prepared for testing. A material sample of Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle 109\times 109} 

mm is represented by an assembly of randomly compacted 2100 discs of radii 1–1.5 mm. It is shown in [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. 15 is characterized by a porosity of Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle 13}

%.


Set-up of the numerical model of uncofined compression test
Figure 15: Set-up of the numerical model of uncofined compression test


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 Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle x} 

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.

The purpose of the present study was to determine the microscopic parameters for sandstone with the following macroscopic properties: Young modulus Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle E=14} 

GPa, Poisson ratio Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \nu=0.3}
 and unconfined compressive strength Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \sigma _{c}=60}
MPa. The stress-strain relationship shown in Fig. 16 satisfying these requirements has been obtained with the following set of micromechanical parameters: contact stiffness in the normal and tangential directions Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle k_n=k_T=20}
 GPa, Coulomb friction coefficient Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \mu=0.839}
and cohesive bond strengths in the normal and tangential direction, Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle R_n=0.1}
MN/m and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle R_T=1}
MN/m, respectively. The failure mode of the specimen obtained in the simulation is shown in Fig. 17. Particles with broken cohesive bonds are marked with different colour. Comparison of Figs. 14 and 17 show that numerical analysis yields a failure mode  similar to that observed in experiments.
Stress–strain relationship for uniaxial unconfined compression test
Figure 16: Stress–strain relationship for uniaxial unconfined compression test
Failure mode obtained in the simulation of unconfined compression test
Figure 17: Failure mode obtained in the simulation of unconfined compression test

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–strain relationship obtained in this simulation is plotted in Fig. 18. The elastic modulus for tension and tensile strength obtained from the plots are the following: Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle E_t=11.5} 

GPa, Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \sigma _{t}=12.7}
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. 19.
Stress–strain relationship for uniaxial direct tension test
Figure 18: Stress–strain relationship for uniaxial direct tension test
Failure mode obtained in the simulation of direct tension test
Figure 19: Failure mode obtained in the simulation of direct tension test

7.2 Simulation of rock cutting

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. 20. The material sample Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle 109\times 109} 

mm is represented by an assembly of randomly compacted 2100 discs of radii 1–1.5 mm. Contact interface between particles representing the rock is characterized with the stiffness in the normal and tangential directions Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle k_n=k_T=20}
 GPa, the cohesive bond strengths Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle R_n=0.1}
MN/m, Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle R_T=1}
MN/m, and the friction coefficient Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \mu=0.839}

. The tool has been modelled with 6800 linear triangles. The following elasto-plastic parameters have been assumed for the tool material: Young modulus Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle E=2\cdot 10^5} 

MPa, Poisson's ratio Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \nu=0.3}

, yield stress Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \sigma _Y=600} 

MPa and hardening modulus Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle 300}
MPa. The following parameters have been assumed for the tool-rock interface: contact stiffness modulus Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle k_n=20}
GPa, Coulomb friction coefficient Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \mu = 0.839}

.

Rock cutting problem. 2100 discrete elements are used to model the rock whereas 6800 linear triangles are used to discretize the tool.
Figure 20: Rock cutting problem. 2100 discrete elements are used to model the rock whereas 6800 linear triangles are used to discretize the tool.


Cutting process has been carried out with prescribed horizontal velocity of the tool of 4 m/s. The cutting is shown in Fig. 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.

Variation of the cutting force in the initial phase of cutting is shown in Fig. 22. The Mises stress and efffective plastic strain distributions in the tool in the initial phase of cutting are shown in Figs. 23 and 24, respectively.

Draft Samper 102318540-cut-fem12-bb-0015.png Draft Samper 102318540-cut-fem12-bb-0025.png
Draft Samper 102318540-cut-fem12-bb-0065.png Rock cutting process analysed using the FEM/DEM model. The failure mode
Figure 21: Rock cutting process analysed using the FEM/DEM model. The failure mode
Variation of the cutting force in the initial phase of cutting
Figure 22: Variation of the cutting force in the initial phase of cutting
Rock cutting process analysed using the FEM/DEM model. The Mises stresses in the tool in the initial phase of cutting
Figure 23: Rock cutting process analysed using the FEM/DEM model. The Mises stresses in the tool in the initial phase of cutting
Rock cutting process analysed using the FEM/DEM model. The effective plastic strain in the tool in the initial phase of cutting
Figure 24: Rock cutting process analysed using the FEM/DEM model. The effective plastic strain in the tool in the initial phase of cutting

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

Rock cutting process analyzed with the tool discretized with 4500 discrete elements
Figure 25: Rock cutting process analyzed with the tool discretized with 4500 discrete elements

The required material constants Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle k} 

and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle H}
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 Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle k=0.0025}
has been taken in the wear simulation.

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.

Tool shapes at different stages of wear are shown in Fig. 26. As much as 40% of the tool mass has been lost due to abrasive wear. Figure 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.

Draft Samper 102318540-tshape-a.png Draft Samper 102318540-tshape-b.png
Draft Samper 102318540-tshape-c.png Evolution of the cutting tool shape due to wear during the rock cutting process
Figure 26: Evolution of the cutting tool shape due to wear during the rock cutting process
Mass loss fraction in the tool versus analysis time –- comparison of experimental and simulation wear curves
Figure 27: Mass loss fraction in the tool versus analysis time –- comparison of experimental and simulation wear curves

7.3 Strip punch test

The strip punch test investigated experimentally by Dede [41] and studied numerically using a FE model by Klerck [42] is analysed here using the DEM model. Prismatic quartzite specimens of dimensions Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle 80\times 80\times 80} 

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. 28. A confining pressure of  4.5 MPa was applied to the specimen sides.
Strip punch geometry and loading configuration
Figure 28: Strip punch geometry and loading configuration


The experimental fracture pattern of the specimen is shown in Fig. 29 and the experimental axial stress versus the specimen axial strain is plotted in Fig. 30.
Experimental fracture pattern, after [41]
Figure 29: Experimental fracture pattern, after [41]
Experimental axial stress vs. axial strain plot, after [41]
Figure 30: Experimental axial stress vs. axial strain plot, after [41]


The 2D discrete element model used in our analysis is shown in Fig. 31. The square rock specimen Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle 80\times 80} 

mm is represented by an assembly of randomly compacted 2100 discs of radii 0.8–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.
2D discrete element model
Figure 31: 2D discrete element model


Required microscopic parameters for the DEM model corresponding to the macroscopic properties of quartzite rock given in [41] are: Young's modulus Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle E=68} 

GPA, compressive strength Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle 200}
MPa and tensile strength Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle 27}
MPa. Those values have been obtained using the procedure presented in example 7.1 for the following microscopic parameters: the stiffness in the normal and tangential directions Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle k_n=k_T=100}
 GPa, the cohesive bond strengths Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle R_n=0.25}
MN/m, Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle R_T=2.5}
MN/m, and the friction coefficient Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \mu=0.839}

.

The numerical failure pattern is shown in Fig. 32. The cracking pattern obtained is similar to experimental results shown in Fig. 29. The triangular crushed zone below the punch and the axial cleavage crack are reproduced in the numerical simulation.

Quite a good agreement can also be observed between experimental (Fig. 30) and numerical (Fig. 33) critical values of axial stresses.

Draft Samper 102318540-strip-punch2-bb-035.png Draft Samper 102318540-strip-punch2-bb-040.png
Draft Samper 102318540-strip-punch2-bb-050.png Numerical failure pattern, after [41]
Figure 32: Numerical failure pattern, after [41]
Numerical axial stress vs. axial strain plot
Figure 33: Numerical axial stress vs. axial strain plot

7.4 Impact of projectile against rock plate

The impact of the projectile against a thick rock plate as defined in Fig. 34 has been studied in order to illustrate possibilities of combined FEM/DEM modelling for this kind of problems.

Impact of projectile – definition of geometry
Figure 34: Impact of projectile – definition of geometry

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 [3]. Projectile has been assumed to be made of steel with the following properties: Young's modulus Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle E=2\cdot 10^5} 

MPa,  Poisson's ratio Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \nu=0.3}

, yield stress Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \sigma _Y=400} 

MPA and isotropic linear hardening modulus Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle H=300}
MPA. The macroscopic parameters of rock has been assumed the same as in the previous example (section 7.3). Again  the same collection of discrete elements has been used.

Fracturing of the plate under the impact as well as projectile deformation is shown in Figs. 35 and 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. 35 and 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. 37.

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.

Draft Samper 102318540-rock-bara-ps3t-c2-ini.png Draft Samper 102318540-rock-bara-ps3t-c2-bb-024.png Draft Samper 102318540-rock-bara-ps3t-c2-bb-040.png
Draft Samper 102318540-rock-bara-ps3t-c2-bb-056.png Draft Samper 102318540-rock-bara-ps3t-c2-bb-080.png Impact of projectile – numerical results obtained using linear triangular elements based on FIC formulation
Figure 35: Impact of projectile – numerical results obtained using linear triangular elements based on FIC formulation
Draft Samper 102318540-rock-bara-ps3-c2-ini.png Draft Samper 102318540-rock-bara-ps3-c2-bb-024.png Draft Samper 102318540-rock-bara-ps3-c2-bb-040.png
Draft Samper 102318540-rock-bara-ps3-c2-bb-056.png Draft Samper 102318540-rock-bara-ps3-c2-bb-080.png Impact of projectile – numerical results obtained using mixed bilinear Q1/P0 quadrilateral elements
Figure 36: Impact of projectile – numerical results obtained using mixed bilinear Q1/P0 quadrilateral elements
Draft Samper 102318540-rock-bara-ps3t-c2-eps-112.png Impact of projectile – effective plastic strain distribution; comparison of numerical results obtained with FIC-based linear triangles with those obtained using mixed bilinear Q1/P0 quadrilateral elements
Figure 37: Impact of projectile – effective plastic strain distribution; comparison of numerical results obtained with FIC-based linear triangles with those obtained using mixed bilinear Q1/P0 quadrilateral elements

7.5 Study of pipe ovalization

Interaction of soil and pipe leading to pipe ovalization has been studied. This example demonstrates another possibility of combined FEM/DEM modelling in geomechanics.

Pipe ovalization – geometry definition
Figure 38: Pipe ovalization – geometry definition


The definition of the problem geometry is shown in Fig. 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 Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle Q} 

increases linearly as
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): Q =k \bar{Q}

where Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle k} 

is load factor and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \bar{Q}}
is the constant reference load defined by a set of uniformly distributed concentrated forces Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle F_i}


Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): \bar{Q}=\{ F_i\} , i=1,12

The following values of reference forces have been assumed:

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): 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
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): 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}.

The load has been transmitted to soil by a kind of diaphragm with hinges in the places of force application.

The following material properties for the pipe have been assumed: Young's modulus Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle E=2\cdot 10^5} 

MPa, Poisson's coefficient Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \nu=0.3}

, density Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \rho=7800} 

kg/mFailed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle ^3}

, yield stress Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \sigma _Y = 300} 

MPa, isotropic hardening modulus Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle H=300}
MPa. The pipe has been discretized with 42 plane strain triangular elements based on the FIC formulation.

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 Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \rho=2800} 

kg/mFailed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle ^3}

, the stiffness in the normal and tangential directions Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle k_n=k_T=200} 

 GPa, Coulomb friction coefficient Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \mu=0.8}

.

The pipe–soil interface has been defined by the stiffness in the normal and tangential directions Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle k_n=k_T=20} 

 GPa,  Coulomb friction coefficient Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \mu=0.3}

.

Numerical results for the deformed pipe shape are shown in Fig. 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. 40.

Pipe cross-sections deformation can be described using an ovalization factor parameter Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle f_{ov}} , defined in terms of the maximum and minimum pipe diameters, Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle D_{max}} 

and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle D_{min}}

, respectively as

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): f_{ov}= \frac{D_{max}-D_{min}}{D_{min}+D_{max}}\cdot 100%

Variation of the ovalization factor for different load levels is shown in Fig. 41. Numerical results agree with those given in [43].


Draft Samper 102318540-pipeovaln3s-wd-geo-0.png Draft Samper 102318540-pipeovaln3s-wd-geo-02.png
Draft Samper 102318540-pipeovaln3s-wd-geo-04.png Pipe ovalization – deformed configuration at different load levels
Figure 39: Pipe ovalization – deformed configuration at different load levels
Draft Samper 102318540-pipeovaln3s-wd-eps-0.png Draft Samper 102318540-pipeovaln3s-wd-eps-02.png
Draft Samper 102318540-pipeovaln3s-wd-eps-04.png Pipe ovalization – deformed configuration
Figure 40: Pipe ovalization – deformed configuration
Pipe ovalization factor
Figure 41: Pipe ovalization factor

Conclusions

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 [17,44].

References

[1] H. Konietzky, editor. Numerical Modeling in Micromechanics via Particle Methods. Balkema Publishers, 2002.

[2] B.K. Cook and R.P. Jensen, editors. Discrete Element Methods. Numerical Modeling of Discontinua. ASCE, 2002.

[3] O.C. Zienkiewicz and R.C. Taylor. The Finite Element Method. Butterworth-Heinemann, Oxford, V edition, 2000.

[4] E. Oñate, J. 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.

[5] E. Oñate, R.L. Taylor, O.C. Zienkiewicz, and J. Rojek. A residual correction method based on finite calculus. In FEMClass of 42 Meeting, 30–31 May, Ibiza (www.femclass42.com), 2002.

[6] P.A. Cundall and O.D.L. Strack. A discrete numerical method for granular assemblies. Geotechnique, 29:47–65, 1979.

[7] C.S. Campbell. Rapid granular flows. Annual Review of Fluid Mechanics, 2:57–92, 1990.

[8] Eng. Comput., 9(2), 1992. Special issue on Discrete Element Methods, Editor: G. Mustoe.

[9] J.R. Williams and R. O'Connor. Discrete Element Simulation and the Contact Problem. Archives Comp. Meth. Engng, 6(4):279–304, 1999.

[10] P.A. Cundall. Formulation of a Three Dimensional Distinct Element Model –- 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–116, 1988.

[11] J. Rojek, E. Oñate, F. Zarate, and J. 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.

[12] J. Argyris. An excursion into large rotations. Comput. Meth. Appl. Mech. Eng., 32:85–155, 1982.

[13] L.M. Taylor and D.S. Preece. Simulation of blasting induced rock motion. Eng. Comput., 9(2):243–252, 1992.

[14] T.J.R. Hughes. The Finite Element Method. Linear Static and Dynamic Analysis. Prentice-Hall, 1987.

[15] J.F. Archard. Contact and rubbing of flat surfaces. J. Appl. Phys., 24(8):981–988, 1953.

[16] E. Rabinowicz. Friction and wear of materials. 1995.

[17] J. Rojek, E. Oñate, F. Zarate, and J. Miquel. Thermomechanical discrete element formulation for wear analysis of rock cutting tools. (to be published).

[18] I. Babuska. The finite element method with lagrangian multipliers. Numer. Math., 20:179–192, 1973.

[19] F. 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–151, 1974.

[20] O.C. Zienkiewicz, S. Qu, R.L. Taylor, and S. Nakazawa. The patch test for mixed formulations. Int. J. Num. Meth. Eng., 23:1873–1883, 1986.

[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–544, 1997.

[22] F. Brezzi and J. Pitkäranta. On the stabilization of finite element approximations of the Stokes problem. In W. Hackbusch, editor, Efficient Solution of Elliptic Problems, Notes on Numerical Fluid Mechanics. Vieweg, Wiesbaden, 1984.

[23] O.C. Zienkiewicz, J. Rojek, R.L. Taylor, and M. Pastor. Triangles and tetrahedra in explicit dynamic codes for solids. Int. J. Num. Meth. Eng., 43:565–583, 1998.

[24] J. Rojek, O.C. Zienkiewicz, E. Oñate, and R.L. Taylor. Simulation of metal forming using new formulation of triangular and tetrahedral elements. In Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): 8^th 

Int. Conf. on Metal Forming 2000, Kraków,   Poland, 2000. Balkema.

[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–189, 1989.

[26] J. Bonet, H. Marriot, and O. Hassan. Stability and comparison of different linear tetrahedral formulations for nearly incompressible explicit dynamic applications. Int. J. Num. Meth. Eng., 50:119–133, 2001.

[27] J. Bonet, H. Marriot, and O. Hassan. An average nodal deformation tetrahedron for large strain explicity dynamic applications. Commun. Num. Meth. in Engrg., 17:551–561, 2001.

[28] R. Codina and J. Blasco. Stabilized finite element method for transient Navier-Stokes equations based on pressure gradient projection. Comput. Meth. Appl. Mech. Eng., 182:287–300, 2000.

[29] R. Codina. Stabilization of incompressibility and convection through orthogonal sub-scales in finite element methods. Comput. Meth. Appl. Mech. Eng., 190:1579–1599, 2000.

[30] M. Chiumenti, Q. Valverde, C. Agelet de Saracibar, and M. Cervera. A stabilized formulation for incompressible elasticity using linear displacement and pressure interpolations. Comput. Meth. Appl. Mech. Eng., 2002.

[31] M. Chiumenti, Q. Valverde, C. Agelet de Saracibar, and M. Cervera. A stabilized formulation for incompressible plasticity using linear triangles and tetrahedra. Int. Journal of Plasticity, 2002.

[32] E. OFailed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \tilde{\mbox{n}}} ate. Derivation of stabilized equations for advective-diffusive transport and fluid flow problems. Comput. Meth. Appl. Mech. Eng., 151:233–267, 1998.

[33] E. OFailed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \tilde{\mbox{n}}} ate. A stabilized finite element method for incompressible flows using a finite increment calculus formulation. Comput. Meth. Appl. Mech. Eng., 182:355–370, 2000.

[34] E. OFailed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \tilde{\mbox{n}}} ate, J. García, and S. Idelhson. Computation of the stabilization parameter for the finite element solution of advective-diffusive problems. Int. J. Num. Meth. Fluids., 25:1385–1407, 1997.

[35] E. OFailed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \tilde{\mbox{n}}} ate, J. García, and S. Idelhson. An alpha-adaptive approach for stabilized finite element solution of advective-diffusive problems. In P. Ladeveze and J.T. Oden, editors, New Adv. in Adaptive Comp. Meth. in Mech. Elsevier, 1998.

[36] E. OFailed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \tilde{\mbox{n}}} ate and M. Manzan. A general procedure for deriving stabilized space-time finite element methods for advective-diffusive problems. Int. J. Num. Meth. Fluids., 31:203 – 221, 1999.

[37] E. OFailed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \tilde{\mbox{n}}} ate and J. García. A finite element method for fluid–structure interaction with surface waves using a finite calculus formulation. Comput. Meth. Appl. Mech. Eng., 191:635–660, 2001.

[38] J. García and E. OFailed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \tilde{\mbox{n}}} ate. An unstructured finite element solver for ship hydrodynamic problems. J. Appl. Mech., 2002.

[39] E. OFailed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \tilde{\mbox{n}}} ate, C. Sacco, and S. Idelhson. A finite point method for incompressible flow problems. Computing and Visualization in Science, 3:67–75, 2000.

[40] H. Huang. Discrete Element Modeling of Tool-Rock Interaction. PhD thesis, University of Minnesota, 1999.

[41] T. Dede. Fracture onset and propagation in layered media. M.Sc. Dissertation, University of the Witwaterstrand, Johannesburg, South Africa, 1997.

[42] P.A. Klerck. The Finite Element Modelling of Discrete Fracture in Quasi-brittle Materials. PhD thesis, University of Wales, Swansea, 2000.

[43] http://geosim.engr.mun.ca/pipesoil.htm.

[44] E. Oñate, C. Recarey, F. Zarate, J. Miquel, and J. Rojek. Characterization of micro-macro parameters in discrete element formulation. (to be published).

APPENDIX A

A.1 Matrices and vectors in FEM equations of motion using the FIC formulation

In FEM equations of motion using the FIC formulation (94)–(96) the following vectors and matrices have been used:

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\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
(A.1)

is the mass matrix

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\bf g}= \int_\Omega {\bf B}^T {\boldsymbol\sigma} d\Omega
(A.2)

is the internal nodal force vector and the rest of matrices and vectors are

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\bf G}_{ij} = \int_{\Omega^e} ({\boldsymbol\nabla} N_i) N_j \, d\Omega a
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): 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
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\bf C} = \left[\begin{matrix} \bar {\bf C}^1 & {\bf 0}&{\bf 0}\\ {\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 (A.3)
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\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
Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\bf f}_i = \int_{\Omega^e} N_i {\bf b} \, d\Omega + \int_\Gamma N_i \bar {\bf t} \, d\Gamma \quad , \quad i,j=1, n_d

In above Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\bf b}=[b_1,b_2,b_3]^T} 

and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle \bar{\bf t} = [\bar t_1, \bar t_2,\bar t_3]^T}

,

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\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]
(A.4)

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\bf B}} 

is the standard infinitesimal strain matrix and Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\bf D}_d}
is the deviatoric constitutive matrix [3]. Note that the expression of Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://mathoid.scipedia.com/localhost/v1/":): {\textstyle {\bf g}}
of eq. (A.2) is adequate for non linear structural analysis.
Back to Top
GET EPUB

Document information

Published on 01/01/2004

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

Web of Science Core Collection® Times cited: 141
Crossref Cited-by Times cited: 146
OpenCitations.net Times cited: 67

Document Score

0

Times cited: 141
Views 214
Recommendations 0

Share this document

claim authorship

Are you one of the authors of this document?

I am an author