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 ${\textstyle i}$-th element we have

where ${\displaystyle {u}}$is the element centroid displacement in a fixed (inertial) coordinate frame ${\displaystyle \mathbf {X} }$, ${\displaystyle \omega }$– the angular velocity, ${\textstyle m}$ – the element mass, ${\textstyle I}$ – the moment of inertia, ${\displaystyle \mathbf {F} }$– the resultant force, and ${\displaystyle \mathbf {T} }$– the resultant moment about the central axes. Vectors ${\displaystyle \mathbf {F} }$ and ${\displaystyle \mathbf {T} }$ are sums of all forces and moments applied to the ${\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 ${\displaystyle \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 ${\displaystyle \mathbf {X} }$, and hence the rotational motion can be easily described in this system.

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 ${\textstyle n}$-th time step is as follows:

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

The vector of incremental rotation ${\textstyle {\Delta \theta }=\{\Delta \theta _{x}\,\Delta \theta _{y}\,\Delta \theta _{z}\}^{\mathrm {T} }}$ is calculated as

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 ${\textstyle n^{2}}$, where ${\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 ${\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 ${\textstyle {\textbf {F}}_{1}}$ and ${\textstyle {\textbf {F}}_{2}}$, which by the Newton's third law satisfy the following relation:

We take ${\textstyle {\textbf {F}}={\textbf {F}}_{1}}$ and decompose ${\displaystyle {\textbf {F}}}$into the normal and tangential components, ${\textstyle {F}_{n}}$and ${\textstyle {\textbf {F}}_{T}}$, respectively (Fig. 2)

where ${\displaystyle {\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 ${\textstyle F_{n}}$ and ${\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 ${\textstyle k_{n}}$ and ${\textstyle k_{T}}$, respectively, the Coulomb friction coefficient ${\textstyle \mu }$, and the contact damping coefficient ${\textstyle c_{n}}$.

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

Figure 3: Model of the contact interface

2.3.2 Normal contact force

The normal contact force ${\textstyle F_{n}}$ is decomposed into the elastic part ${\textstyle F_{ne}}$ and the damping contact force ${\textstyle F_{nd}}$

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 ${\textstyle F_{ne}}$ is proportional to the normal stiffness ${\textstyle k_{n}}$ and to the penetration of the two particle surfaces ${\textstyle u_{rn}}$, i.e.

The penetration is calculated as

where ${\textstyle d}$ is the distance between the particle centres, and ${\textstyle r_{1}}$, ${\textstyle r_{2}}$ their radiae. In the present study no cohesion is allowed, so no tensile normal contact forces are allowed and hence

If ${\textstyle u_{rn}<0}$, Eq. (11) holds, otherwise ${\textstyle F_{ne}=0}$.

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

where ${\textstyle v_{rn}}$ is the normal relative velocity of the centres of the two particles in contact defined by

The damping coefficient ${\textstyle c_{n}}$ can be taken as a fraction of the critical damping ${\textstyle C_{cr}}$ for the system of two rigid bodies with masses ${\textstyle m_{1}}$ and ${\textstyle m_{2}}$, connected with a spring of stiffness ${\textstyle k_{n}}$ ([13]) with

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 ${\textstyle {\textbf {F}}_{T}}$appears by friction opposing the relative motion at the contact point. The relative tangential velocity at the contact point ${\textstyle {\textbf {v}}_{rT}}$ is calculated from the following relationship

with

where ${\textstyle {\dot {\textbf {u}}}_{1}}$, ${\textstyle {\dot {\textbf {u}}}_{2}}$, and ${\textstyle {\omega }_{1}}$, ${\textstyle {\omega }_{2}}$ are the translational and rotational velocities of the particles, and ${\textstyle {\textbf {r}}_{c1}}$ and ${\textstyle {\textbf {r}}_{c2}}$ are the vectors connecting particle centres with contact points.

Figure 4: Friction force vs. relative tangential displacement a) Coulomb law, b) regularized Coulomb law

The relationship between the friction force ${\textstyle \|{\textbf {F}}_{T}\|}$and the relative tangential displacement ${\textstyle u_{rT}}$ for the classical Coulomb model (for a constant normal force ${\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, ${\textstyle {\textbf {v}}_{rT}^{\mathrm {r} }}$ and ${\textstyle {\textbf {v}}_{rT}^{\mathrm {ir} }}$, respectively as:

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

and then the slip condition is checked

If ${\textstyle \phi ^{\mathrm {trial} }\leq {0}}$, we have stick contact and the friction force is assigned the trial value

otherwise (slip contact) a return mapping is performed giving

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 ${\textstyle {F}_{i}^{\mathrm {damp} }}$ and ${\textstyle {\mathbf {T} }_{i}^{\mathrm {damp} }}$ are added to equations of motion (1) and (2)

where

• for viscous damping

• for non-viscous damping

where ${\textstyle \alpha ^{vt}}$, ${\textstyle \alpha ^{vr}}$, ${\textstyle \alpha ^{nvt}}$ and ${\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 ${\textstyle {\Delta t}}$, i.e.

where ${\displaystyle {\Delta t}_{cr}}$ is a critical time step determined by the highest natural frequency of the system ${\displaystyle \omega _{max}}$

If damping exists, the critical time increment is given by

where ${\textstyle \xi }$ is the fraction of the critical damping corresponding to the highest frequency ${\displaystyle \omega _{max}}$. Exact calculation of the highest frequency ${\displaystyle \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

where

and ${\textstyle {\mathbf {k} }_{i}}$ is the stiffness matrix accounting for the contributions from all the penalty constraints active for the ${\textstyle i}$-th particle. Equation (34) defines vectors ${\textstyle {\mathbf {m} }_{i}}$ and ${\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

Equation (33) leads to the eigenproblem

whose eigenvalues ${\textstyle \lambda _{j}}$ (${\textstyle j=1,\ldots ,6}$ in 3D case, and ${\textstyle j=1,2,3}$ for 2D case) are the squared frequencies of the free vibration frequencies:

In a 3D problem, three of six frequencies ${\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:

where it is assumed that the translational motion ${\textstyle u_{n}}$is due to the contact interaction along the normal direction (the spring stiffness ${\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 ${\textstyle k_{T}}$, the rotational stiffness ${\textstyle k_{\theta }}$ can be obtained as

where ${\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:

while the rotational frequency ${\textstyle \omega _{\theta }}$ can be obtained from

where ${\textstyle I_{i}}$ is the rotational inertia of a sphere

and ${\textstyle k_{\theta }}$ is given by Eq. (40). Substituting these expressions into Eq. (42) gives the rotational frequency as

If ${\textstyle k_{T}=k_{n}}$, the rotational frequency ${\textstyle \omega _{\theta }}$ is considerably higher than the translational frequency ${\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:

where ${\textstyle \sigma }$ and ${\textstyle \tau }$ are the normal and tangential contact force, respectively, ${\textstyle k_{n}}$ and ${\textstyle k_{t}}$ are the interface stiffness in the normal and tangential directions and ${\textstyle u_{n}}$ and ${\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:

where ${\textstyle R_{n}}$ and ${\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.

and the (positive) tangential contact force is given by

if ${\textstyle \sigma {<0}}$ or zero otherwise. The friction force is given by Eq. (50) expressing the Coulomb friction law, with ${\textstyle \mu }$ being the Coulomb friction coefficient.

Figure 5: Normal contact force in the elastic perfectly brittle model

Figure 6: Tangential contact force in the elastic perfectly brittle model (tensile normal contact force)

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:

where ${\textstyle \sigma _{i}}$ is the normal or tangential contact force, ${\textstyle u_{i}}$ the total relative displacement, ${\textstyle u_{i}^{p}}$ the plastic component of the relative displacement, ${\textstyle k_{i}}$ the elastic parameters which characterize the contact interface stiffness and ${\textstyle n,t}$ the normal and tangential directions.

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

where ${\textstyle H}$ is the plastic softening modulus (assumed to be positive) and ${\textstyle \alpha }$ is the softening parameter represented by the plastic part of the relative displacement ${\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:

where ${\textstyle \gamma }$ is the absolute value of the slip rate (${\textstyle \gamma {>0}}$), ${\textstyle {\dot {u}}_{i}^{p}}$ the derivative of the plastic relative displacement (normal or tangential) and ${\textstyle \sigma _{i}}$ the contact force in the normal or tangential direction.

Figure 8: Elasto-plastic contact law with softening for the normal direction

Figure 9: Elasto-plastic contact law with softening for the tangential direction

Figure 10: Yield surface for elasto-plastic model

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

where ${\textstyle k_{n}^{D}}$ is the elastic damaged secant modulus and ${\textstyle \omega }$ the scalar damage variable.

The scalar damage variable ${\textstyle \omega }$ is a measure of material damage. For the undamaged state ${\textstyle \omega =0}$ and for a damaged state ${\textstyle 0<\omega \leq 1}$. The scalar damage variable ${\textstyle \omega }$ can be written in the following form:

where ${\textstyle \psi (u_{n})}$ is a function of total relative displacement. For linear strain softening ${\textstyle \psi (u_{n})}$ is defined by

where ${\textstyle R_{n}}$ is the initial tensile strength and ${\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.

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 (${\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 ${\textstyle \sigma }$ is calculated from the linear relationships:

where ${\textstyle u_{n}\leq 0}$ is the normal relative displacement (penetration of one sphere against another). The frictional contact force is evaluated using the regularized Coulomb law

where ${\textstyle \mu }$ is the Coulomb friction coefficient. The frictional dissipation rate ${\textstyle {\dot {D}}}$ is calculated as

where ${\textstyle v_{t}}$ is the relative tangential velocity. The frictional dissipation is used to calculate the wear rate ${\textstyle {\dot {w}}}$ and the heat generated by friction ${\textstyle Q_{gen}}$. This is computed by

where ${\textstyle \chi }$ is the part of the friction work converted to heat (${\textstyle \chi \leq 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 ${\textstyle {\dot {w}}}$ is proportional to the contact pressure ${\textstyle \sigma }$ and to the slip velocity ${\textstyle v_{T}}$

with ${\textstyle H}$ being the hardness of worn surface and ${\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 ${\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:

in which the wear rate is related to frictional dissipation rate ${\textstyle {\dot {D}}}$ and ${\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 ${\textstyle H}$ as a function of the temperature ${\textstyle T}$

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

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]

where ${\textstyle n_{d}}$ is the number of space dimensions of the problems (i.e. ${\textstyle n_{d}=3}$ for 3D)

In (65) and (66) ${\displaystyle {\sigma }_{ij}}$and ${\displaystyle {b}_{i}}$are the stresses and the body forces, respectively and ${\displaystyle {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 ${\textstyle u_{i}}$

and the equilibrium of surface tractions

In the above ${\displaystyle {\bar {u}}_{i}}$and ${\displaystyle {\bar {t}}_{i}}$are prescribed displacements and tractions over the boundaries ${\textstyle \Gamma _{u}}$and ${\textstyle \Gamma _{t}}$, respectively, ${\textstyle n_{i}}$ are the components of the unit normal vector and ${\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

where ${\displaystyle {\delta }_{ij}}$is the Kronnecker delta function. The linear elastic constitutive equations for the deviatoric stresses ${\textstyle s_{ij}}$ are written as

where ${\textstyle G}$is the shear modulus,

The constitutive equation for the pressure ${\textstyle p}$ can be written for an arbitrary domain of finite size using the FIC formulation as [30]

where ${\textstyle K}$ is the bulk modulus of the material.

Note that for ${\textstyle h_{k}\to 0}$ the standard relationship between the pressure and the volumetric strain of the infinitesimal theory ${\textstyle (p=K\varepsilon _{v})}$ is found.

For an incompressible material ${\textstyle K\to \infty }$ and Eq. (72) yields

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

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

where ${\displaystyle {r}_{i}}$is defined by Eq. (66) and

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

where

The coefficients ${\textstyle \tau _{i}}$ in Eq. (79) are also referred to as intrinsic time parameters per unit mass (their dimensions are ${\textstyle {t^{2}m^{3} \over Kg}}$ where ${\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

where ${\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 ${\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 ${\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

Pressure constitutive equation

where ${\textstyle \Delta p=p^{n+1}-p^{n}}$ and ${\textstyle \Delta u_{i}=u_{i}^{n+1}-u_{i}^{n}}$ are the increments of pressure and displacements, respectively. As usual ${\textstyle (\cdot )^{n}}$ denotes values at time ${\textstyle t_{n}}$.

In the derivation of Eq. (83) we have accepted that ${\textstyle \Delta r_{i}=r_{i}^{n+1}\equiv r_{i}}$ as the infinitesimal equilibrium equations are assumed to be satisfied at time ${\textstyle t_{n}}$ (and hence ${\textstyle r_{i}^{n}=0}$).

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

Equilibrium'

Pressure constitutive equation

where ${\textstyle \delta u_{i}}$ and ${\textstyle q}$ are arbitrary test functions representing virtual displacements and virtual pressure fields, respectively.