This study formulates a finite element algorithm for frictional contact of solid materials, accommodating finite deformation and sliding. The algorithm uses a penalty method regularized with an augmented Lagrangian scheme to enforce contact constraints in a nonmortar surface-to-surface approach. Use of a novel kinematical approach to contact detection and enforcement of frictional constraints allows solution of complex problems previously requiring mortar methods or contact smoothing algorithms. Patch tests are satisfied to a high degree of accuracy with a single-pass penalty method, ensuring formulation errors do not affect the solution. The accuracy of the implementation is verified with Hertzian contact, and illustrations demonstrating the ability to handle large deformations and sliding are presented and validated against prior literature. A biomechanically relevant example addressing finger friction during grasping demonstrates the utility of the proposed algorithm. The algorithm is implemented in the open source software febio, and the source code is made available to the general public.

## Introduction

Frictional contact of biological tissues is fundamental in biomechanics, and consideration of friction is often necessary to fully describe relevant physical phenomena. Within the body friction plays an important role in the mechanics of tendons and ligaments, for example, in the flexor tendons of the finger [1,2], gliding resistance through the carpal tunnel [3,4], and between the cruciate ligaments and bone in the knee [5]. Externally, skin friction has been studied both to reduce workplace accidents [6] and to provide pleasing sensory feedback on commercial products [7–9]. Medical applications include friction between internal organs and surgical microdevices [10] and frictional resistance during prosthesis implantation [11]. Although some biological tissues may be treated with constitutive models involving fluid phases, this work is particularly concerned with finite element modeling of friction between solid, elastic, or inelastic materials. Other applications in biomechanics, such as friction between porous hydrated materials, e.g., articular cartilage in diarthrodial joints [12,13], will be addressed in a future study.

Due to the inherent complexity of general three-dimensional contact, finite element algorithms for frictional contact have been extensively studied over the past few decades. In contrast to early finite element efforts [14,15] which were restricted to specific discretizations, Laursen and Simo [16] provided a complete continuum formulation for contact that was independent of the numerical discretization, a framework later extended by Pietrzak and Curnier [17]. To resolve inequalities arising from contact constraints, the penalty method was widely adopted [18,19], permitting minor violations of constraints. Augmented Lagrangian schemes of first- [20] and second- [21] order were devised to avoid the ill-conditioning inherent in penalty formulations and provide more exact enforcement of contact constraints. These early efforts were generally formulated as variations of the node-to-segment or node-to-surface (NTS) approach, where contact constraints are imposed pointwise at a finite number of slave nodes [19,22], thus introducing an artificial asymmetry into the contact problem. This technique is unable to pass the patch test [23], where a constant contact pressure needs to be exactly transmitted between contacting surfaces [24]. It also suffers from geometrical difficulties related to contact detection [22]. Segment-to-segment algorithms [25,26], introduced to enforce constraints in an approximate sense in an effort to remedy errors associated with NTS algorithms, led to development of the mortar method [27–35]. An alternative approach has been the development of contact smoothing algorithms [36,37], which attempt to provide tangent plane continuity between adjacent segments, thus mitigating the errors associated with discretization at the cost of computational efficiency. More recently, the concept of isogeometric contact was proposed by Hughes et al. [38], which adopts the shape functions of 3D modeling software for contact analyses to obtain a smoother description of the contact interface. A series of frictional algorithms has been developed based on this approach [39–41].

In general, the finite element algorithms described above are implemented in custom codes and are thus unavailable to the majority of researchers. Furthermore, the simple and commonly used NTS algorithms exhibit severe convergence difficulties for complex problems, while the powerful and advanced methods (e.g., mortar methods, isogeometric analysis) are much more expensive computationally. To handle the wide variety of frictional contacts occurring in biomechanics, many of which involve finite deformation and large sliding within arbitrary geometries employing complex constitutive models, there exists a need for a simple, robust finite element algorithm for frictional contact, implemented in an open-source framework for the biomechanics community.

This study describes the formulation and implementation of a surface-to-surface finite element algorithm for frictional contact that retains the simplicity and physical interpretation of NTS methods, while robust enough to solve problems that otherwise typically require specialized contact smoothing, mortar methods, or higher order shape functions. To accomplish this objective, we propose a novel continuum approach and a new algorithmic formulation, extending our previous work on frictionless biphasic contact [42]. Here, frictional stick will be treated as a case of tied contact, with an exact return mapping to frictional slip which is independent of the single penalty parameter and based on a standard slip criterion. Unlike the traditional approach, no rate form for the frictional traction is prescribed, avoiding the need for complex numerical integration [43,44] or elastoplastic-type algorithms [18]. A main contribution of this work is the use of a contact detection algorithm that avoids many of the pitfalls of classical NTS schemes, allowing smooth results without specialized contact elements or computationally expensive smoothing algorithms [36]. Our frictional contact algorithm is formulated in the spatial frame, using a penalty method regularized with an augmented Lagrangian scheme to enforce contact constraints. To improve convergence rates, the linearization of the gap function used during contact slip is accomplished without the simplifying assumptions made in our earlier work [42], and the novel derivation of the frictional formulation utilizes a cleaner matrix notation (see Supplementary Materials, available under the “Supplemental Data” tab for this paper on the ASME Digital Collection). The implementation of this algorithm is incorporated into the free, open-source finite element software febio^{2} [45], which has been developed to address the needs of the biomechanics community. The validity of the presented formulation and algorithm is demonstrated with a series of benchmark problems that involve finite deformation and large sliding.

The paper is organized as follows: Sec. 2 describes the continuum formulation of the contact problem, for both penalty and augmented Lagrangian regularization schemes. Section 3 describes the finite element implementation and presents a series of numerical examples demonstrating the ability of the formulation to handle a wide variety of friction problems. A discussion of these results is presented in Sec. 4. The detailed linearization and discretization of the results in Sec. 2 may be found in the Supplementary Materials which are available under the “Supplemental Data” tab for this paper on the ASME Digital Collection, along with a description of the Gaussian quadrature integration scheme adopted for the present study.

## Continuum Formulation

### Notation.

The problem to be considered consists of two deformable bodies in contact. In a contact analysis, it is convenient to define one body as the contactor and the other as the target, typically referred to as the slave and master, respectively, in the literature. Each point on the slave can be assigned an intersection on the master. As this study uses a modification of the gap function and projection method described by Ateshian et al. [42], which differs from the conventional approach, the designation of surfaces as “slave” and “master” will not be employed, and surfaces will be referred to as *primary* and *secondary*. A subscripted or superscripted index (*i*), *i* = 1, 2, will be used to denote a quantity on either the primary or secondary surface, respectively. Greek indices refer to surface quantities and thus vary implicitly from 1 to 2, with repeated Greek indices implying summation. The standard derivation of the virtual work of contact forces is briefly summarized in the Appendix using the notation of this study.

### Contact Integral.

*γ*

^{(1)}only, yielding

Equation (2.1) is commonly referred to as the contact integral.

To evaluate the directional derivatives of *δG _{c}* along increments in the displacements $\Delta u(i)$ of

*γ*

^{(}

^{i}^{)}, as required for an iterative technique such as Newton’s method, it is necessary to formulate the integration over an invariant domain so that the directional derivative may be brought inside the integral without concern for variations in the domain of integration. In our approach, the integral is formulated over the invariant parametric space of

*γ*

^{(1)}[42,46].

*γ*

^{(}

^{i}^{)}is expressed in parametric form using coordinates $\eta (i)\alpha $, and material points

*X*

^{(}

^{i}^{)}are identified through their parametric coordinates. On each surface

*γ*

^{(}

^{i}^{)}, covariant basis vectors are given by

*γ*

^{(}

^{i}^{)}as it deforms over time

*t*, in terms of contravariant surface parametric coordinates $\eta (i)\alpha $. These covariant basis vectors are tangent to

*γ*

^{(}

^{i}^{)}, and it follows that the unit outward normal to

*γ*

^{(}

^{i}^{)}is

*γ*

^{(}

^{i}^{)}is evaluated as

*γ*

^{(1)}, and integration is performed over points

*X*

^{(1)}with prescribed parametric coordinates $\eta (1)\alpha $. Since $\Gamma \eta (1)$ represents a material frame, it is possible to linearize

*δG*by applying the directional derivative operator directly to the integrand

_{c}*f*

To proceed with this linearization, it is necessary to formulate the kinematics of points on *γ*^{(}^{i}^{)} and provide expressions for the contact traction **t**^{(1)} that can differentiate between frictional stick and slip.

### Slip Kinematics.

*γ*

^{(1)}and

*γ*

^{(2)}move relative to one another as their configurations evolve. Performing a contact analysis requires mapping points between these surfaces. For the spatial position $x(1)(\eta (1)\alpha ,t)$ of each material point on the primary surface, we define the intersection point $x(2)(\eta (2)\alpha ,t)$ on the secondary surface as the point intersected by a ray directed along the unit outward normal to the primary surface $n(1)(\eta (1)\alpha ,t)$

*g*is defined to be positive when the surfaces

*γ*

^{(1)}and

*γ*

^{(2)}are separated, and negative when they penetrate

The ray intersects *γ*^{(2)} at a spatial position **x**^{(2)} through which different material points *X*^{(2)} identified by parametric coordinates $\eta (2)\alpha $ may convect. Computationally, this ray intersection and contact detection is performed with an Octree method [47].

The projection approach of Eq. (2.8), described in our previous studies [42,48] and commonly termed *ray-tracing* [49], can be characterized as an inverse projection relative to the classical contact mechanics approach used for NTS contact [18,19]. Although the definition of this projection and its associated gap function is not new, it has typically been employed mostly for mortar contact (e.g., as in Tur et al. [33]). The benefits associated with a projection method such as Eq. (2.8) for nonmortar contact have been developed in detail [49]. Here, it suffices to note that avoiding a reliance on secondary surface normal vectors eliminates many contact-searching difficulties that plague NTS algorithms [22] and also serves to greatly reduce the complexity of the linearizations and the resulting stiffness matrices.

### Stick Kinematics.

*X*

^{(2)}on

*γ*

^{(2)}does not evolve during the iterative solution process, allowing the development of kinematics of sticking contact. Implicit in this condition is the assumption that the contact projection was previously resolved, and thus, contact searching is not performed again; rather, the contact point on

*γ*

^{(2)}is given by the parametric coordinates of intersection $\eta (2)\alpha $ found from the previous time point, which will be denoted as $\eta (2)p\alpha $, with the subscripted

*p*referring to the previous time. Thus, we write the spatial position of

*X*

^{(2)}in stick as

**g**

*is a vectorial gap function*

_{s}Here, **g*** _{s}* is the vectorial distance, at the current time

*t*, between material points which were in contact at the previous time-step; for perfect stick, we must have

**g**

*=*

_{s}**0**. In a finite element implementation, however, it is important to note that the intersection point $xs(2)$ is not in general the point that would be found from a ray directed from $\eta (1)\alpha $ along the unit outward normal

**n**

^{(1)}. This is because stick will not be enforced exactly, so the points

**x**

^{(1)}and $xs(2)$ will separate slightly when using a penalty method for enforcing this constraint. How to minimize this separation and enforce perfect stick behavior is the subject of Sec 2.6.

### Velocities.

The kinematics developed above can be used to formulate velocities. The formulation for stick developed in Sec. 2.6 does not rely on velocity constraints; therefore, we are only concerned with velocities of opposing contact points in slip. The development of velocities presented below anticipates the need for a relative velocity as required by Coulomb’s law of kinetic friction, which aligns the friction force with this slip direction.

*γ*

^{(1)}, the velocity

**v**

^{(1)}of these points on the primary surface is evaluated from the material time derivative in the material frame

*γ*

^{(2)}, the total velocity

**v**

^{(2)}at the intersection point on

*γ*

^{(2)}needs to be evaluated using the material time derivative in the spatial frame

Here, $\u2202x(2)/\u2202t$ represents the velocity of the intersection point on *γ*^{(2)}, whereas $\eta \u02d9(2)\alpha $ are the contravariant components of the convective velocity of material passing through this intersection point. In effect, $\eta \u02d9(2)\alpha g\alpha (2)$ represents the relative (slip) velocity between the material on *γ*^{(2)} and that on *γ*^{(1)}. Importantly, as noted below when performing time discretization and linearization, by definition $\u2202x(2)/\u2202t$ is evaluated while keeping $\eta (2)\alpha $ constant.

*γ*

^{(1)}and

*γ*

^{(2)}[49,50]

**n**

^{(1)}in the parametric material frame of

*γ*

^{(1)}, evaluated from Eq. (2.3) as

**P**

*as*

_{N}Previous authors have utilized $vr=\eta \u02d9(2)\alpha g\alpha (2)$ directly instead of the right-hand-side of Eq. (2.14); this expression requires the evaluation of time derivatives of parametric coordinates [15,16,41], which necessitates special integration algorithms to handle crossing element boundaries [44]. The relative velocity measure on the right-hand-side of Eq. (2.14) obviates the need for any such special treatment. In particular, the choice of $\u2202x(2)/\u2202t$ ensures that element boundaries will never be crossed when calculating this velocity, since it is evaluated while keeping $\eta (2)\alpha $ constant.

**s**

^{(1)}as the unit vector of the projection $PN\xb7vr$ of

**v**

*onto the tangent plane of*

^{r}*γ*

^{(1)}

These definitions of contact kinematics may now be used to formulate frictional contact.

### Coulomb Frictional Contact.

**t**

^{(}

^{i}^{)}on the opposing surfaces is determined by the sticking or slipping behavior. For Coulomb friction, the relationship between sticking and slipping is described by a slip criterion Φ, where on the primary surface

*μ*is the friction coefficient. The value of the slip criterion determines the stick–slip status

Algorithmically, stick and slip are typically based on a predictor–corrector approach derived from an analogy with elastoplasticity [43], leading to constitutive relations for the rate of the traction and thus requiring numerical integration [44]. Variations of this approach have been utilized in differing forms [15,16,41].

In contrast, this presentation proposes to treat stick and slip separately, controlled by an exact return mapping based on the slip criterion. The return mapping defines a rule for correcting a calculated traction which exceeds the slip limit and is thus not permissible. Stick will be treated as a special case of a tied interface, whereas in slip the traction will be directly prescribed. The formulation of Coulomb frictional contact is presented for both penalty and augmented Lagrangian regularization schemes.

#### Penalty Scheme.

*γ*

^{(1)}and

*γ*

^{(2)}that were in contact at the previous time must remain in contact. The stick traction may be obtained by penalizing any relative motion between such points

where *ε* is the penalty parameter, and we have employed Eq. (2.11).

Here, we have achieved an exact expression for the tangential traction in slip and remark that since *t _{n}* is strictly negative in contact, the frictional contact traction $\mu tns(1)$ acting on

*γ*

^{(1)}is in the direction opposing the motion of

*γ*

^{(1)}relative to

*γ*

^{(2)}. A trial state and return map, presented in Sec. 2.6.3, is employed to differentiate between stick and slip.

#### Augmented Lagrangian Scheme.

The augmented Lagrangian scheme employed in this study is first-order and utilizes Uzawa’s algorithm [51], where multipliers are updated outside of the Newton step, producing a double loop algorithm (see the texts by Laursen [18] and Wriggers [19] for a discussion of Uzawa’s algorithm applied to frictional contact problems). Such an approach preserves the quadratic convergence of Newton’s method near solution points. The presented scheme is a modification of the approach suggested by Simo and Laursen [20].

**g**

*,*

_{s}*g*

*λ*is the normal Lagrange multiplier. The total traction vector in slip is then defined to be

_{n}where *t _{n}* is given by Eq. (2.25). In this approach, the Lagrange multiplier $\lambda s$ augments the traction

**t**

^{(1)}in stick, but in slip only the normal component of traction

*t*is augmented by

_{n}*λ*, and the tangential traction is directly prescribed from the augmented normal component. This approach has the advantage of preserving an exact mapping to the proper tangential traction in slip, which is consistent with the augmented normal traction. As in the penalty case, a trial state and return map controlled by the slip criterion is employed to differentiate between stick and slip, presented in Sec. 2.6.3.

_{n}*λ*and $\lambda s$ are held constant during each Newton step. Outside of the Newton loop, in this study, we propose a novel update scheme where one of these multipliers is considered active and is updated from the kinematic data (

_{n}**g**

*or*

_{s}*g*), and the other is considered passive and is derived from the active multiplier. The contact status is determined via Eq. (2.20). If the current status is stick (Φ < 0), we update $\lambda s$ (active) and derive

*λ*(passive)

_{n}*λ*(active) and derive $\lambda s$ (passive)

_{n}An active-passive strategy for the multipliers ensures consistency when the contact status switches between stick and slip and is made possible due to this formulation’s use of a single penalty parameter *ε*, in conjunction with an exact return mapping for slip.

*P*

_{tol}specifies the largest allowable change. Convergence is achieved when

*L*represents the total norm of the active multipliers across the contact surface at augmentation step

^{r}*r*, calculated by summing all the individual norms

*e*th element face, and $lk,er$ is the norm of the active Lagrange multiplier at the

*k*th integration point on the

*e*th element of

*γ*

^{(1)}at augmentation step

*r*, defined by

*G*

_{tol}at every location. Convergence of the augmentations requires

where $|gs|$ is associated with points currently sticking and $|g|$ with those in slip. The tolerance *G*_{tol} has units of length, allowing enforcement of the nonpenetration and stick constraints to arbitrarily small precision.

#### Stick–Slip Algorithm.

where *t _{n}* is given by either Eq. (2.22) or (2.25). In the case of first contact, a trial stick traction cannot be calculated, as stick requires a previous intersection point. Consequently, first contact is treated as a case of slip in this framework. The alternative of treating first contact as frictionless is unsatisfying, as the lack of friction at the first iteration can lead to premature failure and thus precludes the modeling of certain problems that rely on frictional forces for stability (such as load-control analyses). After one iteration, the traction can be evaluated via the return map described above.

Computationally, care must be taken to ensure that augmentation does not unnecessarily change the stick–slip status. For normal contact, the update $\lambda n\u2190\lambda n+\epsilon g$ will augment the normal traction until the nonpenetration contact constraint is adequately satisfied, with no adverse consequence if *λ _{n}* slightly overshoots the final target value during an intermediate augmentation. For tangential contact however, augmentation of the tangential traction that overshoots the final target value may cross the boundary between stick and slip, thus changing the nature of the solution. For example, at the first augmentation step, the stick traction is $t(1)=\epsilon gs$ and the multiplier $\lambda s$ is augmented from its initial zero value to $\lambda s=\epsilon gs$ according to Eq. (2.27). Thus, at the start of the next iteration, when

**g**

*has not yet changed, the traction is calculated as $t(1)=\lambda s+\epsilon gs=2\epsilon gs$ according to Eq. (2.24), which essentially counts the gap function twice and can in some cases shift the contact from stick to slip. To circumvent potential error introduced by this step, our implementation freezes the stick–slip status until the completion of the first iteration following an augmentation step. After the first iteration, the gap function has been reduced by the augmentation, and the traction is split appropriately between the remaining gap and the multiplier. In this way, the double-counting of the gap function does not unnecessarily shift the contact status, but augmentation is able to modify the stick–slip status if accurate enforcement of the contact constraints requires it.*

_{s}## Examples and Verifications

### Constitutive Relations.

where the constants *G* and *λ* are material coefficients, $J=det\u2009F$ is the Jacobian of the deformation, and $I1=tr\u2009C$, where **C** is the right Cauchy-Green tensor. For each example, the material properties given will be Young’s modulus *E* and Poisson’s ratio *ν*, where *E* = 2*G*(1 + *ν*) and *λ* = 2*Gν*/(1 – 2*ν*). Unless otherwise noted, all benchmark problems utilize the augmented Lagrangian regularization scheme and a two-pass algorithm which enforces contact constraints twice, swapping the primary and secondary surfaces the second time and applying the same procedure [42]. All examples are quasi-static, thus references to time henceforth refer to increments in the load in quasi-static analyses.

### Finite Element Implementation.

*ε*was computed using

where the summation is taken over all *N* element faces representing the contact surfaces in the contact pair, *V _{i}* is the volume of each element,

*A*is the area of the element face on the contact surface,

_{i}*E*is a measure of the average Young’s modulus, and

_{i}*α*is a user-defined, nondimensional scale factor [42]. For all examples, the scale factor

*α*is reported, along with augmented Lagrangian tolerances

*G*

_{tol}and

*P*

_{tol}presented in Sec. 2.6.2, when applicable.

In contrast to some algorithms that automatically detect surface contact pairs that evolve in real-time as a result of penetration or destructive tests [52], febio requires users to pre-identify potential contact pairs when prescribing boundary conditions on the finite element model. In displacement-control analyses, contacting bodies may be initially separated; however, in quasi-static load control analyses, an initial overlap of contacting surfaces must be prescribed so that force balance can be satisfied at the initial time-step.

### Patch Tests.

When a contact surface between two bodies is subjected to a spatially constant stress field, this stress field should be exactly transmitted from one body to another [24], as any variations in the stress at the interface are due to errors in the formulation employed. In this study, we perform the patch test proposed by Neto et al. [37]. Two cubes of identical geometry and material properties are pressed against each other under frictionless conditions (Fig. 1(a)). The cubes are 10 mm on each side, with material properties *E* = 100 MPa and *ν* = 0.3. The bottom surface of the lower cube is fixed, and displacements of the lateral surfaces of all cubes are constrained in the normal direction. A uniform pressure is prescribed by applying a vertical displacement *u _{z}* = −1 mm to the top surface of the upper cube between

*t*= 0 and

*t*= 1 s, ramped up over ten uniform time steps.

The distribution of the vertical Cauchy stress component *T _{zz}* at the final time is shown in Fig. 1, for different choices of the primary and secondary surface. A single-pass analysis was utilized without an augmented Lagrangian regularization scheme. The patch test fails when the coarser bottom surface is taken as the primary (

*α*= 2000, Fig. 1(b)), yet the patch test is passed to five decimal places when the finer top surface is taken as the primary (

*α*= 1 × 10

^{9}, Fig. 1(c)). When keeping the top surface as the primary and utilizing an augmented Lagrangian regularization with

*P*

_{tol}= 1 × 10

^{–6}, the penalty scale factor can be reduced to

*α*= 1000 while still retaining the same accuracy.

It has been well documented that classical single-pass NTS algorithms introduce artificial asymmetry into the contact interface and cannot pass the patch test [23], although extensions of the NTS framework have been attempted. Chen and Hisada [53] developed a modified NTS approach that can pass the patch test with a single-pass algorithm, although their formulation required a Lagrange multiplier approach. More recently, Zavarise and De Lorenzis [54] were able to modify a single-pass NTS algorithm to pass the patch test with a penalty method for two-dimensional contact. In the present formulation, a single-pass algorithm with the penalty method can pass the patch test for general three-dimensional contact problems. That the patch test can be passed to such high accuracy in a single-pass analysis demonstrates the strength of the proposed algorithm and offers a powerful advantage over classical NTS formulations.

With a two-pass algorithm, similar accuracy can be obtained with *α* = 1 × 10^{8} with a penalty method and, with augmented Lagrangian regularization, *α* = 1 × 10^{4} and *P*_{tol} = 1 × 10^{−4}. Although a two-pass algorithm is not required to pass the patch test and ensure accurate transmission of contact tractions, in practice the use of a two-pass analysis often offers additional accuracy. Furthermore, augmented Lagrangian regularization allows for identical results with significantly reduced penalty factors, thus eliminating ill-conditioning associated with large penalties. However, if one contact surface is much coarser than the other, or if contact occurs between a rigid body and a deformable surface, a two-pass analysis may not be advised. Potential adverse consequences of an augmented Langrangian method will be treated in Sec 3.9. Consequently, although the following examples generally use a two-pass analysis with augmented Lagrangian regularization, certain problems may call for a different method of analysis and constraint enforcement. Fortunately, the results of the patch tests presented in this section demonstrate that the modeler has the freedom to choose between different methods without fear of compromising accuracy.

### Hertzian Contact.

Cylinder on plane Hertzian contact was used to verify the accuracy of the implementation, due to the availability of a closed form analytical solution [55] which can be compared to calculated stresses. A schematic of this problem is shown in Fig. 2, where we use *L* = 8 mm, *R* = 50 mm, and unit depth. An isotropic, linearly elastic material was used for both contacting bodies, with *E* = 115 GPa and *ν* = 0.32. The friction coefficient was specified to be *μ* = 0.7.

In this plane strain analysis, a normal load *P* = 800 N is ramped up over 1 s and then held constant in time. Subsequently, a tangential load *Q* = 504 N is applied linearly between 1 s and 2 s. The choice of *Q* = 0.9*μP* created distinct stick and slip regions, allowing evaluation of the algorithm’s ability to switch between contact cases. Since frictional problems are history-dependent, this problem was solved in two time steps in order to match the assumptions made in deriving the analytical solution. Both coarse (Fig. 2(b)) and fine (not shown) meshes were employed. The coarse mesh utilized 7680 elements, whereas the fine mesh utilized 66,240 elements. Penalty and augmented Lagrangian parameters were taken to be *α* = 0.1 and *P*_{tol} = 0.05 for the coarse mesh and *α* = 1 and *P*_{tol} = 0.05 for the fine mesh. In each case, a minimum of six augmentations was enforced, regardless of whether the tolerance condition was met at a prior augmentation.

*T _{yy}* stress contours after the final time-step are displayed on the coarse mesh (Fig. 3). Both coarse and fine meshes display strong agreement with the analytical solution (Fig. 4). The coarse mesh fails to capture the contact width exactly, although that is likely related to the absence of an element at the exact contact radius. Upon refining the mesh, the exact contact width is captured to a high degree of accuracy. Interestingly, the spatial distribution of contact traction is not exactly symmetric, in opposition to the analytical solution, although the asymmetry is very minor. This is a consequence of employing a finite deformation code to replicate a problem solved under the assumption of infinitesimal deformation.

### Stacked Blocks in Contact.

This example presents the compression of four stacked blocks with dissimilar material properties and meshes, demonstrating the ability of the algorithm to handle problems involving loss of contact. Each block is a 2 mm cube (Fig. 5(a)). The first and third blocks from the top have material properties *E* = 10 MPa and *ν* = 0.1. The second and fourth blocks have *E* = 0.3 MPa and *ν* = 0.4. The friction coefficient and penalty scale factor between the third and fourth blocks are *μ* = 1 and *α* = 10, with *μ* = 0.05 and *α* = 1 at the other two interfaces. *P*_{tol} = 0.2 was selected for convergence of the augmented Lagrangian scheme. A prescribed displacement *u _{z}* = –2 mm was applied to the top face of the top block in 10 time steps. The bottom face of the bottom block is constrained in all directions, but no other constraints are applied to the model.

As a consequence of the friction coefficients employed, the bottom two blocks stick, while sliding develops at the other interfaces (Fig. 5). This problem is adapted from an illustration provided by Laursen et al. [29], who noted that the dissimilar meshes and the transverse gap developed on the top contact surface typically leads to failure of NTS algorithms. The contact detection strategy employed in this study encounters no difficulties with the overhanging blocks.

Next, this benchmark problem is extended by introducing finite rotations. We repeat the same analysis, but now also apply a rotation *R _{z}* =

*π*/2 to the top surface of the top block simultaneously with

*u*. Due to the additional instability caused by the rotations, the friction coefficient at the top two interfaces was increased to

_{z}*μ*= 0.08, while the penalty was decreased to

*α*= 0.1. At the bottom interface, the penalty was decreased to

*α*= 5. The resulting analysis involved finite deformation edge-to-face contact, with nodes and integration points coming into and out of contact. A series of snapshots of the twisting behavior is presented in Fig. 6.

The bottom two blocks are still sticking at the end of the twisting analysis, despite the increased traction on the interface. The top three blocks are experiencing both edge-to-edge and edge-to-face contact, with large transverse gaps forming due to the twisting behavior. Despite the challenging nature of the benchmarks, the distribution of the contact traction is very clean for both examples presented here (Movies S1 and S2 and Figs. S1 and S2 are available under the “Supplemental Data” tab for this paper on the ASME Digital Collection).

### Twisting Contact Between Hemisphere and Box.

In this example described by Sauer and De Lorenzis [41], a thick-walled hollow hemisphere indents a deformable box and then rotates, developing large shear deformations and twisting contact between the hemisphere and box.

The mesh and geometry are shown in Fig. 7, where we have adopted *L* = 2 mm. The box is a cube, with material properties *E* = 10 MPa and *ν* = 0.3. The thick-walled hemisphere’s material properties are *E* = 50 MPa and *ν* = 0.3, and the friction coefficient between the hemisphere and box is taken to be *μ* = 0.5. Penalty and augmented Lagrangian parameters were taken to be *α* = 10 and *P*_{tol} = 0.05, respectively. A vertical ramp displacement *u _{z}* = −1 mm is applied to the top surface of the hemisphere between

*t*= 0 s and

*t*= 1 s, and then held constant. A rotation

*R*=

_{z}*π*is subsequently applied to the top face of the hemisphere from

*t*= 1 s to

*t*= 10 s. This problem is solved with 100 evenly spaced time steps.

During twisting, contact between the box and hemisphere transitions from predominantly sticking initially to predominantly slipping behavior by the end of the analysis (Fig. 7). Consequently, the applied torque required to rotate the hemisphere remains essentially constant beyond the transition to slip. Figure 8 displays the twisting torque as a function of prescribed rotation angle, for both this study and that of Sauer and De Lorenzis [41], demonstrating strong agreement. The slight oscillations observed in the mostly flat torque curve from 100 deg to the end of twisting are due to stick–slip instabilities in the motion. Visually, the box displays a type of slight chattering behavior during this time.

### Compression of Concentric Spheres.

The next example has been used to motivate mortar formulations for contact [31] and thus represents a challenging benchmark. We consider compression of two thick-walled concentric spheres between rigid planes (Fig. 9). Frictional sliding contact is present between the faces of the two spheres, which are modeled in octant symmetry. The use of incompatible meshes forces gap discontinuities during both contact and sliding. This problem was originally presented by Puso and Laursen [31], who noted premature failure of NTS algorithms even in frictionless contact, and was later solved by Areias et al. [56] with friction, albeit to a smaller applied displacement.

Both spheres have material properties *E* = 1 MPa and *ν* = 0.3. Sliding contact between the planes and the spheres is modeled as frictionless, with *α* = 10 and *P*_{tol} = 0.2. Since the planes are rigid, a single-pass analysis was used for these contacts, with the sphere selected as the primary surface. Frictional contact occurs between the two spheres, where we examined friction coefficients *μ* = 0, 0.5 and 2. As different friction coefficients introduced markedly different interfacial behaviors, the penalty and augmented Lagrangian parameters differed for each friction coefficient. For *μ* = {0, 0.5, 2}, we selected *α* = {1, 3, 8} and *P*_{tol} = {0.2, 0.2, 0.05}, respectively. The vertical displacement of the top plane is prescribed linearly, reaching a final displacement *u _{z}* = −10 mm in 10 s. The problem is solved in 100 evenly spaced time steps for all three values of the friction coefficient.

Snapshots of the deformed mesh are presented in Fig. 9, for the case *μ* = 0. For *μ* = 0.5, the buckling of the inner sphere occurs slightly later, and in the *μ* = 2 case, not at all. This can be observed in plots of the vertical reaction force during compression (Fig. 10(a)), where the force is much larger in the *μ* = 2 case due to lack of this buckling behavior. As Areias et al. [56] only performed their analysis to *u _{z}* = −9 mm, these results were also compared to those obtained with the commercial finite element software abaqus

^{3}for all three values of friction coefficient (Figs. 10(b)–10(d)); thus, this benchmark also served to evaluate the agreement between the present algorithm implemented in febio and proprietary commercial software. Notably, the data from Areias et al. display a much smaller reaction force for

*μ*= 2, although agreement between febio and abaqus is quite good. Areias et al. do not present details on buckling of the inner sphere, but it is possible that the discrepancy may have emerged if their model showed buckling at

*μ*= 2, which would have reduced the magnitude of the normal compression force.

### Deep Indentation of Square Blocks.

This benchmark problem, originally proposed by Temizer [34] in a mortar framework, introduces sharp corners and edges into a three-dimensional finite deformation problem. In this example, a small stiff block is fully compressed into a larger, softer body (Fig. 11), developing contact between the sides of the small block and the surface of the larger block.

The top block has material properties *E* = 100 MPa and *ν* = 0.3, whereas those of the bottom block are *E* = 1 MPa and *ν* = 0.3. A linearly increasing vertical displacement *u _{z}* = −0.6 mm was applied to the top face of the top block over 1 s. This example was solved to a slightly larger final displacement, as Temizer begins with the blocks separated by 0.0375 mm [34]. Although the problem was initially presented as frictionless, we consider

*μ*= 0.2 here to evaluate the performance of the proposed friction algorithm. The penalty scale factor was

*α*= 150 and augmentation was performed on the gap, with

*G*

_{tol}= 0.0009 mm. Since the large indentation can cause the softer material to rise up against the indenting small block, and since the small block has a coarse mesh in the vertical direction, a single-pass analysis is employed for this benchmark, using the surface of the large block as the primary surface. 100 evenly spaced time steps were utilized in this analysis.

This example is challenging due to the sharp corners in the indenting block, which tend to introduce singularities in the contact pressure. Indeed, examples involving indentation (and sliding) with sharp corners have long been used to motivate mortar formulations for contact [32,35]. Our implementation in febio is able to complete this analysis without difficulty. Examination of the distribution of contact pressure (Fig. 12) displays a large peak in the corners and elevated values along the edges, as may be expected. Comparing the vertical reaction force to results obtained by Temizer [34] (scaled by the correct factor of 10), a reasonable agreement is seen (Fig. 13), although the slopes differ slightly. This may reflect differences between the mortar strategy employed by Temizer and the formulation of this study, differences in the constitutive model adopted for neo-Hookean behavior, errors in the febio solution introduced by the sharp corners in this example, or the effect of friction on the reaction force. To address the latter concern, an identical analysis with *μ* = 0 was performed (results not shown), and the vertical reaction force from the frictionless case was virtually indistinguishable from that of *μ* = 0.2. In an effort to understand the discrepancies between the results obtained Temizer et al. and those from febio, this problem was solved in abaqus as well, with *μ* = 0.2. The vertical reaction force obtained from abaqus (Fig. 13), although displaying the same trend as the other models, lies in between the two curves, suggesting that results from this problem may be influenced by differences in contact formulations or constitutive models.

### Shallow Ironing.

This example demonstrates shallow ironing of a rectangular block by a smaller indenter and represents a classic benchmark that can be used to evaluate contact smoothness. Various forms of this example have been presented by many authors, e.g., see Refs. [28], [33], [41], and [49]. For this study, we adopt the geometry of Neto et al. [37]. In this plane strain example, a stiff indenter with rounded corners is pressed vertically into a larger, soft block, and then displaced horizontally. The discretization employed for both surfaces has been noted to lead to chatter [28], or oscillations in the response, as contact shifts from one discretized facet to another with no tangent plane continuity. This example demonstrates the ability of the presented formulation to ameliorate discretization-induced noise in a framework that does not explicitly adopt smoothing techniques.

The geometry and mesh are given in Fig. 14. The material properties for the indenter are *E* = 1000 MPa and *ν* = 0.3, while those of the block are *E* = 100 MPa and *ν* = 0.3. The friction coefficient was taken to be *μ* = 0.3. A downward displacement *u _{y}* = −10 mm is applied to the top surface of the indenter over 1 s and then held while a horizontal displacement

*u*= 250 mm is prescribed between

_{x}*t*= 1 s and

*t*= 2 s. As Neto et al. [37] do not provide information on the number of time steps utilized in their study, we performed simulations using 200, 400, 800, and 1600 uniformly spaced time steps. 200 time steps appeared to introduce artificial smoothing, but results obtained with 800 and 1600 time steps did not differ appreciably from results with 400 time steps. Therefore, for what follows, 400 uniformly spaced time steps were used. In keeping with the presentation of Neto et al., the surface of the block was selected as the primary surface and a single-pass analysis was employed. Contact constraints were enforced with both penalty (

*α*= 5) and augmented Lagrangian (

*α*= 2,

*P*

_{tol}= 0.2) methods. Snapshots of the deformed mesh at the end of indentation, halfway through sliding, and upon completion of sliding are presented in Fig. 15 for the augmented Lagrangian analysis.

The study of Neto et al. [37] examined the effect of contact smoothing with Nagata patches on the behavior of an NTS algorithm and compared results for faceted and smoothed surfaces. The evolution of the normal (*F _{N}*) and tangential (

*F*) components of the contact force, for this study and that of Neto et al., is presented in Fig. 16 for both penalty and augmented Lagrangian regularizations. The results obtained by the penalty method, although generally in agreement with those of Neto et al., display weaker enforcement of sticking behavior during the indentation phase (indicated by the arrow in Fig. 16), necessitating the use of augmentations. Both normal and tangential contact forces predicted by febio with an augmented Lagrangian method are in excellent agreement with the smoothed results of Neto et al. (in that study, the faceted surface of the indenter caused the iterative procedure to diverge around

_{T}*u*= 85 mm if smoothing was not performed). If we calculate the friction coefficient from these data as $\mu =FT/FN$, we find

_{x}*μ*= 0.298 ± 0.002 for the penalty method and

*μ*= 0.296 ± 0.001 for the augmented Lagrangian method between

*t*= 1.2 s and

*t*= 1.7 s, corresponding to the pure sliding regime. This result is in excellent agreement with the prescribed value

*μ*= 0.3.

A close-up of the tangential forces in Fig. 16 reveals substantial differences between schemes during sliding (Fig. 17). Oscillations in the tangential force during sliding are numerical and caused by the discretization. The faceted data from Neto et al. show large chatter, as elements on their master surface come into and out of contact. Their smoothing technique, employing Nagata patches to reduce discretization error on the master surface, is successful and greatly reduces the magnitude of the oscillations. In contrast, the present formulation, which employs no smoothing algorithms and relies on a wholly different contact detection approach, exhibits oscillations of much smaller magnitude. The present formulation is independent of the secondary surface normal, but the projection is still between two bilinearly interpolated surface facets and involves a degree of discretization error. The use of very small time steps amplifies the nonsmoothness of the gap function, which is due to the contact associated with a primary surface integration point intersecting a linearly discretized surface. Although the augmented Lagrangian scheme reduced the magnitude of the oscillations during the steady sliding regime, it created larger spikes in the response as sliding was initiated; these spikes were not present in the penalty method. Nearly identical spikes were observed regardless of the number of time steps. Although these larger spikes are still smaller than the magnitude of the oscillations in the smoothed data of Neto et al., and are quickly damped, it is instructive to consider their origin and whether they reveal flaws in the augmented Lagrangian scheme presented in Sec. 2.6.2.

As we shall see, the origin of these spikes reveals potential pitfalls associated with employing augmented Lagrangian methods for contact of curved surfaces modeled with linear elements. Penalty and augmented Lagrangian methods are used to enforce contact constraints, chiefly the nonpenetration constraint for normal contact and the constraint of no relative motion during stick. For this benchmark, a theoretically perfect enforcement of these contact constraints is actually not desirable. Since we are approximating a smoothly curved surface with linear elements, a degree of discretization error is introduced in the secondary surface; additionally, although the rectangular primary surface suffers no discretization error relating to its geometry, deformations of this surface are limited by the linear elements employed. Thus, for this benchmark, a perfect enforcement of contact constraints would in fact correspond physically to sliding contact between a polygonal indenter and a surface comprised of straight links connected by flexible joints. The penalty method, as it enforces the constraints imperfectly, actually implicitly provides for a degree of smoothing and comes closer to replicating the physical scenario of interest. In contrast, augmented Lagrangian methods, though a powerful tool to reduce unphysical constraint violations, can dramatically increase discretization-induced errors under special circumstances. This mechanism is not purely a frictional phenomenon. Indeed, setting *μ* = 0 in the penalty and augmented Lagrangian models presented above shows a similar increase in oscillatory noise as augmentation tolerances are tightened (results not shown).

To evaluate this hypothesis more fully, the geometry of Fig. 14 was reproduced with 20-node quadratic hexahedral elements with eight integration points per surface (hex20), although the number of elements was halved along each side. Consequently, 200 evenly spaced time steps were utilized in order to obtain the same rate of crossing element boundaries seen in the linear case. For this analysis, *α* = 2 and *P*_{tol} = 0.2. As shown in Fig. 18, utilizing hex20 elements nearly entirely eliminates these initial spikes and ameliorates the oscillations during sliding, although they still remain present. The remaining oscillations occur because although quadratic elements allow for individual element faces to conform to one another smoothly, they still do not enforce tangent or normal continuity between elements. For the hex20 simulation, we find *μ* = 0.299±0.0005 between *t* = 1.2 s and *t* = 1.7 s, virtually identical to the prescribed *μ* = 0.3. As the number of elements was halved along each direction, the number of time steps was also halved, and a 20-node hexahedral variant with only eight integration points per surface was used, there was no extra computational cost associated with this analysis.

### Pinching of a Hollow Ball.

This last example simulates compression of a hollow ball between two fingers, demonstrating the biomechanical relevance of incorporating frictional interactions into computational models (Fig. 19(a)). Two hollow tubes represent deformable fingers. Each tube has an outer radius of 1 cm, an inner radius of 0.5 cm, a length of 10 cm, and is capped by a hemisphere with outer radius 1 cm and inner radius 0.5 cm. To obtain a reasonable response for the soft pads of the fingers, material properties for each tube were taken to be *E* = 1 MPa and *ν* = 0.3, which may be considered average values for biological tissues [57]. To take into account the rigid bone underlying the soft tissue of the finger, a rigid body was attached to the inner surfaces of each tube, preventing deformation. A revolute joint connects the tubes at the intersection of their centerlines, enabling them to come together in an approximation of a pinching motion. A hollow sphere is placed initially between the two fingers, in grazing contact with the finger surfaces. The sphere has an outer radius of 2.5 cm and an inner radius of 2.25 cm, and is modeled as a nearly incompressible Mooney–Rivlin material where *c*_{1} = 1.25 MPa, *c*_{2} = 0, and the bulk modulus is *k* = 1250 MPa. Contact was defined separately between each finger and the sphere and a single-pass analysis was used, where the finger was the primary surface for each contact interface. Augmented Lagrangian parameters were *α* = 4 and *G*_{tol} = 0.01 cm, and the friction coefficient was taken to be *μ* = 0.9. A 15 degree rotation of the revolute joint comprised the pinching behavior and was applied in 20 uniformly spaced time steps over 1 s.

This model displays several interesting behaviors. During compression, the ball sticks to the finger surfaces and deforms (Figs. 19(b) and 19(c)). Around *t* = 0.7 s, the sphere buckles inward and large sections lose contact with each finger. However, the finger has deformed, as it is sandwiched between the rigid material simulating bone on the inside and the ball surface on the outside, and this deformation increases the contact surface area and thus the friction force holding the ball in place. Such behavior cannot be observed in the absence of friction. In a frictionless simulation (results not shown), the model managed two time steps and failed at *t* = 0.06 s, when the ball lost contact with the fingers and jumped away. The inclusion of friction allows this model to develop deformation of the finger (Fig. 19(d)), which further traps the ball by adding a plowing component to the friction force, and demonstrates the necessity of including frictional interactions if one wants to model certain physically relevant behavior.

## Discussion

This paper has presented the formulation, implementation, and verification of a surface-to-surface finite element algorithm for frictional contact of solid materials that can accommodate large deformation and sliding within a conceptually simple framework. A distinguishing characteristic of the proposed formulation is the ability to solve challenging problems without recourse to mortar methods or contact smoothing, thus retaining the simple physical interpretation of classical NTS formulations and avoiding computational overhead associated with more complex algorithms. Consequently, special contact elements and large numbers of integration points are not required, and the contact algorithm functions efficiently on faces of standard linear elements, such as the eight-node hexahedral bricks used in the examples presented in Sec. 3. This formulation endows the presented algorithm with a robustness that allows solution of complex three-dimensional finite deformation problems without significantly increased computational demands, facilitating its incorporation into challenging contact problems typically encountered in biomechanics.

The numerical examples presented in Sec. 3 were selected to verify the accuracy of the formulation and validate it against common test problems in the literature. The results demonstrate the ability of the proposed algorithm to handle a wide variety of challenging frictional contact problems, and a robust mechanism to differentiate between areas of stick and slip. Patch tests were passed to five decimal places with a single-pass algorithm using only the penalty method (Fig. 1), and the use of two-pass algorithms and augmented Lagrangian regularization achieved the same accuracy with significantly reduced penalty values. The algorithm was verified with Hertzian contact, where analytical solutions were in excellent agreement with finite element results (Fig. 4). Additional examples demonstrated the ability to handle three-dimensional large-deformation contact problems (Figs. 7 and 9), including the biomechanically relevant problem of object pinching (Fig. 19). Importantly, the present approach was able to solve problems typically associated with failure of NTS algorithms (Figs. 5 and 11), and also provided results which could previously only be obtained with contact smoothing techniques (Fig. 17). The proposed augmented Lagrangian scheme of Sec. 2.6.2 was shown to be robust, generally needing two or fewer augmentations on average to meet the required tolerances. Exceptions include the patch tests where an extremely demanding tolerance was required and the Hertzian contact benchmark where we enforced a minimum of six augmentations. A full summary of the average number of augmentations required for convergence is available in Table S1 available under the “Supplemental Data” tab for this paper on the ASME Digital Collection, for all benchmarks presented in Sec. 3 utilizing an augmented Lagrangian scheme.

One weakness of the proposed algorithm derives from the formulation, which prioritized a robust, fast algorithm over the ability to solve the most demanding problems. To maintain speed and robustness, the contact algorithm was developed to provide excellent results with standard linear finite elements of the type most commonly adopted for complex problems. Consequently, although able to solve numerous problems typically used in the literature to motivate more complex algorithms, the proposed formulation did encounter difficulties with the problem of deep ironing with a square indenter (see, e.g., Ref. [35] for the problem description). The use of bilinear interpolation and linear finite elements led to difficulties in wrapping individual surface facets around the sharp corner of the indenter and produced convergence issues leading to failure at *t* = 1.6 s, although the problem could be solved if a slight radius was applied to the sharp corners of the indenting block (results not shown). Although the presented algorithm develops a powerful, robust method to solve a broad variety of challenging contact problems in a computationally efficient manner, it may not represent the most technically advanced or complex framework. For the purpose of advancing the capabilities of an open-source software tool, this type of general purpose algorithm is desirable, although the tradeoff is a lack of ability to solve certain highly specialized contact problems and may be considered a limitation.

This study addresses an important computational need in the biomechanics community. Importantly, the algorithm is implemented in the open source finite element code febio, which is available to the general public and contains specialized constitutive models and boundary conditions specific to biomechanics problems. Incorporation of a robust frictional contact algorithm within this framework provides a valuable computational tool to the biomechanics community and represents a necessary first step toward the development of a finite element algorithm for frictional contact of biphasic materials.

## Disclaimer

The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health or the National Science Foundation.

## Funding Data

National Institute of General Medical Sciences (Grant No. GM083925).

National Science Foundation Graduate Research Fellowship Program (DGE-11-44155).

### Appendix: Principle of Virtual Work

*b*is the domain of interest in the current configuration,

**T**is the Cauchy stress,

*δ*

**v**is a virtual velocity, and

*dv*is an elemental volume [46]. By applying the divergence theorem, this expression may be rearranged as

*∂b*,

**n**represents the unit outward normal to

*∂b*, and

*da*represents an elemental area of

*∂b*. For convenience, we split the virtual work into its internal and external parts, $\delta W=\delta Wint\u2212\delta Wext$, where

*δW*at some trial solution for the motion $\chi $, along an increment in displacement $\Delta u$ in $\chi $

where $Df[\Delta u]$ represents the directional derivative of *f* along $\Delta u$. The directional derivatives of the external virtual work depend on the type of boundary conditions present. Conditions relevant to contact analyses are discussed next.

*b*consisting of two bodies

*b*

^{(1)}and

*b*

^{(2)}with respective boundaries

*∂b*

^{(1)}and

*∂b*

^{(2)}. The two bodies are in contact over portions of

*∂b*

^{(1)}and

*∂b*

^{(2)}, respectively denoted

*γ*

^{(1)}and

*γ*

^{(2)}. According to Eq. (A4), the contribution of contact to the external virtual work may be written as

where $\delta v(i)$ is a virtual velocity, **t**^{(}^{i}^{)} represents the traction on *γ*^{(}^{i}^{)}, and *da*^{(}^{i}^{)} is an elemental area of *γ*^{(}^{i}^{)}. In contact analyses, the tractions on *γ*^{(}^{i}^{)} are equal and opposite, $t(1)=\u2212t(2)$, and the contact surfaces are shared, hence we may select one surface to perform the integration over.