Coulomb friction

Abaqus provides an extended version of the classical isotropic Coulomb friction model for use with all contact analysis cababilities. The extensions include an additional limit on the allowable shear stress, anisotropy, and the definition of a “secant” friction coefficient.

The following topics are discussed:

Elastic stick formulation
Frictional stress due to forced sliding in static analysis
Friction contributions for complex eigenvalue extraction analysis
Viscous stick formulation for steady-state transport
Exact stick formulation
Anisotropic friction
Viscous damping

Elastic stick formulation

In the elastic stick formulation in Abaqus/Standard, the “elastic” tangential slip $γ_{i}^{e l}$ is defined as the reversible relative tangential motion from the point of zero shear stress. The elastic slip is related to the interface shear stress with the relation

τ_{i} = k_{s} γ_{i}^{e l},

where $k_{s}$ is the (current) “stiffness in stick,” which follows from the relation

k_{s} = \frac{τ_{c r i t}}{γ_{c r i t}} .

Since $τ_{c r i t}$ may be dependent on contact pressure, slip rate, average surface temperature at the contact point, and field variables, $k_{s}$ may change during the analysis. The behavior remains elastic as long as the equivalent stress does not exceed the critical stress; hence,

γ_{i}^{e l} (t + Δ t) = γ_{i}^{e l} (t) + Δ γ_{i} .

Consistent linearization of this expression yields

\begin{matrix} (1) & \begin{matrix} d τ_{i} & = k_{s} d γ_{i} if τ_{c r i t} = τ_{m a x}, \\ d τ_{i} & = k_{s} d γ_{i} + \frac{τ_{i}}{τ_{c r i t}} (μ + p \frac{\partial μ}{\partial p}) d p if τ_{c r i t} = μ p . \end{matrix} \end{matrix}

The contributions from the contact pressure are nonsymmetric for the second case. Since the slip rate is zero in the elastic stick formulation, derivatives with respect to the slip velocity are not needed.

The above expressions hold if the equivalent shear stress remains less than the critical stress. If the equivalent stress exceeds the critical stress, slip must be taken into consideration so that the condition $τ_{e q} = τ_{c r i t}$ is satisfied. Let the starting situation be characterized by the elastic slip ${\bar{γ}}_{i}^{e l}$ . The critical stress at the end of the increment follows from the contact pressure, $p$ , and the slip rate, ${\dot{γ}}_{e q}^{s l}$ .

Let the (as yet unknown) elastic slip at the end of the increment be $γ_{i}^{e l}$ and the slip increment be $Δ γ_{i}^{s l}$ . Consistency requires that

Δ γ_{i} = γ_{i}^{e l} - {\bar{γ}}_{i}^{e l} + Δ γ_{i}^{s l},

and the shear stress at the end of the increment follows from the elasticity relation,

τ_{i} = k_{s} γ_{i}^{e l} = \frac{τ_{c r i t}}{γ_{c r i t}} γ_{i}^{e l} .

The slip increment is related to the stress at the end of the increment with the backward difference approach:

Δ γ_{i}^{s l} = \frac{τ_{i}}{τ_{c r i t}} Δ γ_{e q}^{s l} .

With these equations and the critical stress equality $τ_{e q} = τ_{c r i t},$ it is possible to solve for $γ_{i}^{e l}$ , $Δ γ_{i}^{s l}$ , and $τ_{i}$ . Elimination of $γ_{i}^{e l}$ and $Δ γ_{i}^{s l}$ from the above equations yields,

Δ γ_{i} = \frac{τ_{i}}{τ_{c r i t}} γ_{c r i t} - {\bar{γ}}_{i}^{e l} + \frac{τ_{i}}{τ_{c r i t}} Δ γ_{e q}^{s l}

τ_{i} = \frac{{\bar{γ}}_{i}^{e l} + Δ γ_{i}}{γ_{c r i t} + Δ γ_{e q}^{s l}} τ_{c r i t} .

It is convenient to define the “elastic predictor strain”

γ_{i}^{p r} \overset{def}{=} {\bar{γ}}_{i}^{e l} + Δ γ_{i},

which simplifies the expression for the stress to

τ_{i} = \frac{γ_{i}^{p r}}{γ_{c r i t} + Δ γ_{e q}^{s l}} τ_{c r i t} .

Substitution in the critical stress equality yields

Δ γ_{e q}^{s l} = γ_{e q}^{p r} - γ_{c r i t},

where

γ_{e q}^{p r} \overset{def}{=} \sqrt{{(γ_{1}^{p r})}^{2} + {(γ_{2}^{p r})}^{2}} .

Substitution in the expression for $τ_{i}$ and introduction of the normalized slip direction $n_{i} = γ_{i}^{p r} / γ_{e q}^{p r}$ furnishes the final result

\begin{matrix} (2) & τ_{i} = n_{i} τ_{c r i t} . \end{matrix}

Here $τ_{c r i t}$ is a function of the slip rate, which is obtained with

{\dot{γ}}_{e q}^{s l} = Δ γ_{e q}^{s l} / Δ t,

where $Δ t$ is the time increment in a static analysis. In the case of dynamics with the Hilber-Hughes-Taylor time integration operator, $Δ t$ is scaled by the Hilber-Hughes-Taylor time integration operator parameters, $γ$ and $β$ .

For the iterative solution scheme this equation must be linearized. Some straightforward algebra yields

\begin{matrix} d τ_{i} = & \frac{d γ_{i}^{p r}}{γ_{e q}^{p r}} τ_{c r i t} - n_{i} n_{j} \frac{τ_{c r i t}}{γ_{e q}^{p r}} d γ_{j}^{p r} + n_{i} d τ_{c r i t} \\ = & (δ_{i j} - n_{i} n_{j}) \frac{τ_{c r i t}}{γ_{e q}^{p r}} d γ_{j}^{p r} + n_{i} (\frac{\partial τ_{c r i t}}{\partial p} d p + \frac{\partial τ_{c r i t}}{\partial {\dot{γ}}_{e q}^{s l}} d {\dot{γ}}_{e q}^{s l}) . \end{matrix}

With the expression for the equivalent slip the final result is

\begin{matrix} (3) & \begin{matrix} d τ_{i} = & (δ_{i j} - n_{i} n_{j}) \frac{τ_{c r i t}}{γ_{e q}^{p r}} d γ_{j} if τ_{crit} = τ_{\max}; \\ d τ_{i} = & (δ_{i j} - n_{i} n_{j}) \frac{τ_{c r i t}}{γ_{e q}^{p r}} d γ_{j} + n_{i} (μ + p \frac{\partial μ}{\partial p}) d p + n_{i} n_{j} \frac{p}{Δ t} \frac{\partial μ}{\partial {\dot{γ}}_{e q}} d γ_{j} if τ_{crit} = μ p . \end{matrix} \end{matrix}

In this case the unsymmetric terms may have a strong effect on the speed of convergence of the Newton scheme. Hence, use of the unsymmetric equations solver is strongly recommended for the analysis of problems in which sliding friction occurs. In the case of dynamics with the Hilber-Hughes-Taylor time integration operator, $Δ t$ is scaled by the Hilber-Hughes-Taylor time integration operator parameters.

Frictional stress due to forced sliding in static analysis

If the user specifies different velocities for two bodies in contact as a predefined field, it is assumed that the slip velocity follows from the prescribed motion and is independent of the displacement values. The frictional stresses then follow from

τ_{i} = τ_{e q} \frac{{\dot{γ}}_{i}}{{\dot{γ}}_{e q}},

where

τ_{e q} = μ p,

and

{\dot{γ}}_{e q} = \sqrt{{({\dot{γ}}_{1})}^{2} + {({\dot{γ}}_{2})}^{2}} .

${\dot{γ}}_{i}$ are the tangential slip velocities derived from the user-defined fields ( $i = 1, 2$ ).

Friction contributions for complex eigenvalue extraction analysis

At the contact nodes at which the velocity differential is imposed, the linearized shear stress can be expressed in the following form:

\begin{matrix} d τ_{i} = (δ_{i j} - n_{i} n_{j}) \frac{τ_{c r i t}}{γ_{e q}^{p r}} d {\dot{γ}}_{j} + n_{i} (μ + p \frac{\partial μ}{\partial p}) d p + n_{i} n_{j} p \frac{\partial μ}{\partial {\dot{γ}}_{e q}} d {\dot{γ}}_{j} . \end{matrix}

In the equation above the first term is generated by the friction forces acting in the direction perpendicular to the direction of slip. The third term exists if the friction coefficient depends on velocity. In the complex eigenvalue extraction analysis both terms contribute to the damping matrix. The second term contributes to the stiffness matrix.

Viscous stick formulation for steady-state transport

For steady-state transport a viscous stick formulation is used. In this case the “viscous” tangential slip rate, ${\dot{γ}}_{i}^{v}$ , is related to the interface shear stress by

τ_{i} = κ_{s} {\dot{γ}}_{i}^{v},

where $κ_{s}$ is the “stick viscosity,” which follows from the relation

κ_{s} = \frac{τ_{c r i t}}{{\dot{γ}}_{c r i t}} .

The allowable viscous slip is defined as a fraction of the relative tangential velocity

{\dot{γ}}_{c r i t} = 2 F_{f} ω R,

where $F_{f}$ is a user-defined slip tolerance, $ω$ is the angular spinning velocity, and $R$ is the radius of a point on the contact surface in the undeformed configuration.

Exact stick formulation

In Abaqus/Standard Lagrange multipliers are used to enforce exact sticking conditions. A constraint term enforced with Lagrange multipliers $δ q_{i}$ is added to the virtual work statement:

\begin{matrix} (4) & δ Π^{*} = \int_{S} (τ_{i} δ γ_{i} + Δ γ_{i} δ q_{i}) d S, \end{matrix}

where $δ q_{i}$ is a Lagrange multiplier. By taking the rate of change of Equation 4, the rate of virtual work is obtained:

\begin{matrix} (5) & d δ Π^{*} = \int_{S} (d τ_{i} δ γ_{i} + d γ_{i} δ q_{i} + τ_{i} d δ γ_{i}) d S, \end{matrix}

where the last term results from a nonlinear relation between nodal displacements and relative motion $γ_{i}$ . This term is nonzero only for contact with finite sliding (see Finite-sliding interaction between deformable bodies and Finite-sliding interaction between a deformable and a rigid body).

To obtain a complete formulation, a relationship between $Δ γ_{i}$ , $q_{i}$ , and $τ_{i}$ must be defined. For sticking conditions a suitable relation is

τ_{i} = q_{i} + k_{0} Δ γ_{i},

with $k_{0}$ a reference stiffness selected internally in Abaqus. The reference stiffness can be considered to model a spring that is in parallel to the sticking constraint. Since the constraint prevents relative motion, the reference stiffness has no physical significance and is added only to eliminate zero terms on the diagonal of the stiffness matrix that could cause equation solver problems. The variational and rate forms of the equation for $τ_{i}$ are readily obtained as

δ τ_{i} = δ q_{i} + k_{0} δ γ_{i} and d τ_{i} = d q_{i} + k_{0} d γ_{i} .

Substitution in Equation 5 yields

\begin{matrix} (6) & d δ Π^{*} = \int_{S} (k_{0} δ γ_{i} d γ_{i} + δ γ_{i} d q_{i} + δ q_{i} d γ_{i} + τ_{i} d δ γ_{i}) d S . \end{matrix}

If the element is slipping, $Δ γ_{i}$ can take an arbitrary value. Consequently the terms involving $δ q_{i}$ in Equation 4 and Equation 5 vanish. With the backward difference method, the frictional stress is obtained as

\begin{matrix} (7) & τ_{i} = τ_{c r i t} \frac{Δ γ_{i}}{Δ γ_{e q}}, \end{matrix}

and the linearized form is

\begin{matrix} (8) & \begin{matrix} d τ_{i} = & \frac{τ_{c r i t}}{Δ γ_{e q}} (δ_{i j} - n_{i} n_{j}) d γ_{j} if τ_{crit} = τ_{\max}, \\ d τ_{i} = & \frac{τ_{c r i t}}{Δ γ_{e q}} (δ_{i j} - n_{i} n_{j}) d γ_{j} + n_{i} (μ + p \frac{\partial μ}{\partial p}) d p + n_{i} n_{j} \frac{p}{Δ t} \frac{\partial μ}{\partial {\dot{γ}}_{e q}} d γ_{j} if τ_{crit} = μ p . \end{matrix} \end{matrix}

The rate of virtual work can, hence, be written in the form

\begin{matrix} (9) & \begin{matrix} d δ Π^{*} = & \int_{S} [\frac{τ_{c r i t}}{Δ γ_{e q}} (δ_{i j} - n_{i} n_{j}) δ γ_{i} d γ_{j} + τ_{i} d δ γ_{i}] d S if τ_{crit} = τ_{\max}, \\ d δ Π^{*} = & \int_{S} [\frac{τ_{c r i t}}{Δ γ_{e q}} (δ_{i j} - n_{i} n_{j}) δ γ_{i} d γ_{j} + (μ + p \frac{\partial μ}{\partial p}) n_{i} δ γ_{i} d p \\ + \frac{p}{Δ t} \frac{\partial μ}{\partial {\dot{γ}}_{e q}} n_{i} n_{j} δ γ_{i} d γ_{j} + τ_{i} d δ γ_{i}] d S if τ_{crit} = μ p . \end{matrix} \end{matrix}

Observe that for two-dimensional problems the term involving $(δ_{i j} - n_{i} n_{j})$ vanishes. In the case of dynamics with the Hilber-Hughes-Taylor time integration operator, $Δ t$ is scaled by the Hilber-Hughes-Taylor time integration operator parameters.

To determine whether an element sticks or slips, the following tests are used. If the element is currently sticking, it is checked whether the magnitude of the frictional stress exceeds the critical value. Hence, if

τ_{e q} > τ_{c r i t},

the state is changed to slipping.

If the element is slipping, the direction of slip is compared with the frictional stress $τ_{i} (t)$ applied at the start of the iteration. If at the end of the iteration

Δ γ_{i} τ_{i} (t) < 0,

the state is changed to sticking.

Anisotropic friction

In anisotropic friction the friction coefficients in the two (principal) directions are different; hence, the critical shear stresses are different:

τ_{1}^{c r i t} = μ_{1} p and τ_{2}^{c r i t} = μ_{2} p .

The two critical shear stresses are at the extreme points of the friction ellipse, given by the equation

{(\frac{τ_{1}}{μ_{1}})}^{2} + {(\frac{τ_{2}}{μ_{2}})}^{2} = p^{2} .

This suggests definition of the scaled shear stresses,

\begin{matrix} (10) & \bar{τ_{1}} \overset{def}{=} \frac{\bar{μ}}{μ_{1}} τ_{1} and \bar{τ_{2}} \overset{def}{=} \frac{\bar{μ}}{μ_{2}} τ_{2}, \end{matrix}

where $\bar{μ}$ is the average friction coefficient defined by

\bar{μ} = \sqrt{\frac{1}{2} (μ_{1}^{2} + μ_{2}^{2})} .

Hence, the (scaled) equivalent frictional stress is defined by

{\bar{τ}}_{e q} = \sqrt{{\bar{τ_{1}}}^{2} + {\bar{τ}}_{2}^{2}},

and the average critical stress

{\bar{τ}}_{c r i t} = m i n (\bar{μ} p, τ_{m a x}),

where $τ_{m a x}$ is the user-specified limiting shear stress, which is applied to the scaled shear stresses. No slip will occur if ${\bar{τ}}_{e q} < {\bar{τ}}_{c r i t}$ , and slip can occur if ${\bar{τ}}_{e q} = {\bar{τ}}_{c r i t}$ .

If slip occurs, it is assumed that the direction of slip is governed by an associated slip law:

{\dot{γ}}_{1} = λ \frac{\partial {\bar{τ}}_{e q}}{\partial τ_{1}} = λ \frac{\partial {\bar{τ}}_{e q}}{\partial {\bar{τ}}_{1}} \frac{d {\bar{τ}}_{1}}{d τ_{1}} = λ \frac{{\bar{τ}}_{1}}{{\bar{τ}}_{e q}} \frac{\bar{μ}}{μ_{1}} and {\dot{γ}}_{2} = λ \frac{\partial {\bar{τ}}_{e q}}{\partial τ_{2}} = λ \frac{\partial {\bar{τ}}_{e q}}{\partial {\bar{τ}}_{2}} \frac{d {\bar{τ}}_{2}}{d τ_{2}} = λ \frac{{\bar{τ}}_{2}}{{\bar{τ}}_{e q}} \frac{\bar{μ}}{μ_{2}} .

Hence, scaled rates of slip are defined by

\begin{matrix} (11) & {\dot{\bar{γ}}}_{1} \overset{def}{=} \frac{μ_{1}}{\bar{μ}} {\dot{γ}}_{1} and {\dot{\bar{γ}}}_{2} \overset{def}{=} \frac{μ_{2}}{\bar{μ}} {\dot{γ}}_{2}, \end{matrix}

and the (scaled) equivalent slip rate

{\dot{\bar{γ}}}_{e q} = \sqrt{{\dot{\bar{γ}}}_{1}^{2} + {\dot{\bar{γ}}}_{2}^{2}} .

As a result, the direction of the scaled slip and the scaled shear stress coincide:

\frac{{\bar{τ}}_{i}}{{\bar{τ}}_{e q}} = \frac{{\dot{\bar{γ}}}_{i}}{{\dot{\bar{γ}}}_{e q}} .

Since these equations have exactly the same structure as the equations for isotropic friction, the same solution algorithms can be used to get the incremental and linearized equations in terms of the scaled stress and slip. The actual stresses are readily obtained from the scaled stresses with Equation 10, and the linearized equations are scaled with the same factors.

In the anisotropic elastic stick formulation, it is assumed that the stiffness in stick in terms of the scaled stress and slip is constant:

{\bar{k}}_{s} = \frac{{\bar{τ}}_{c r i t}}{γ_{c r i t}};

hence, the scaled stress is related to the scaled elastic slip by the relation

{\bar{τ}}_{i} = {\bar{k}}_{s} {\bar{γ}}_{i}^{e l} .

In terms of the actual stress and strain in the $1$ - and $2$ -directions, this implies

τ_{1} = k_{1} γ_{1}^{e l} and τ_{2} = k_{2} γ_{2}^{e l},

with

k_{1} = {(\frac{μ_{1}}{\bar{μ}})}^{2} {\bar{k}}_{s}, k_{2} = {(\frac{μ_{2}}{\bar{μ}})}^{2} {\bar{k}}_{s} .

Hence, the anisotropy in terms of the stiffness in stick is more pronounced than the anisotropy in terms of critical stress.

Viscous damping

In addition to the friction models described above, Abaqus allows for the definition of a “viscous” shear stress $τ_{i}^{v}$ that is proportional to the relative tangential velocity ${\dot{γ}}_{i}$ . This viscous damping in the tangential direction occurs if

damping in the contact direction is included and the tangent fraction is nonzero; or
contact controls are used to stabilize rigid body modes automatically in multi-body contact analysis.

The viscous damping stress is proportional to the tangential damping coefficient $μ^{v}$ , which is a function of the overclosure as follows:

μ^{v} (h) = {\begin{matrix} 0 & for h < - c_{0} \\ μ_{0}^{v} (h + c_{0}) / (1 - η) c_{0} & for - c_{0} \leq h < - η c_{0} \\ μ_{0}^{v} & for h \geq - η c_{0}, \end{matrix}

where $μ_{0}^{v}$ is the value of the tangential damping coefficient at zero overclosure and $η$ is the fraction of the overclosure interval $[- c_{0}, 0]$ over which the damping coefficient is equal to $μ_{0}^{v}$ .

The virtual work contribution associated with the viscous shear stress is

δ Π = τ_{i}^{v} δ γ_{i} .

The contribution to the stiffness matrix for the Newton solution is given by the linearized form of the virtual work contribution:

\begin{matrix} d δ Π & = δ γ_{i} d τ_{i}^{v} = δ γ_{i} (\frac{\partial τ_{i}^{v}}{\partial h} d h + \frac{\partial τ_{i}^{v}}{\partial {\dot{γ}}_{j}} \frac{\partial {\dot{γ}}_{j}}{\partial γ_{k}} d γ_{k}) \\ = δ γ_{i} \frac{d μ^{v}}{d h} {\dot{γ}}_{i} d h + δ γ_{i} μ^{v} \frac{\partial {\dot{γ}}_{i}}{\partial γ_{k}} d γ_{k}, \end{matrix}

where

\frac{d μ^{v}}{d h} = {\begin{matrix} 0 & for h < - c_{0} \\ μ_{0}^{v} / (1 - η) c_{0} & for - c_{0} \leq h < - η c_{0} \\ 0 & for h \geq - η c_{0} . \end{matrix}

In static analysis the velocity is defined as the displacement increment divided by the time increment. Therefore, ${\dot{γ}}_{j} = Δ γ_{j} / Δ t$ , and the stiffness contribution reduces to

d δ Π = δ γ_{i} \frac{d μ^{v}}{d h} \frac{Δ γ_{i}}{Δ t} d h + δ γ_{i} \frac{μ^{v}}{Δ t} d γ_{i} .

The previous expression also applies to dynamics with the backward Euler time integration operator. In the case of dynamics with the Hilber-Hughes-Taylor time integration operator, $\partial {\dot{γ}}_{j} / \partial γ_{k}$ is defined by the dynamic time integration operator and the stiffness contribution can be written as

d δ Π = δ γ_{i} \frac{d μ}{d h} {\dot{γ}}_{i} d h + δ γ_{i} \frac{γ μ}{β Δ t} d γ_{i},

where $γ$ and $β$ are the Hilber-Hughes-Taylor time integration operator parameters. Viscous damping cannot be used in a Riks analysis since velocity is not defined.