## Abstract

The frictional response of porous and permeable hydrated biological tissues such as articular cartilage is significantly dependent on interstitial fluid pressurization. To model this response, it is common to represent such tissues as biphasic materials, consisting of a binary mixture of a porous solid matrix and an interstitial fluid. However, no computational algorithms currently exist in either commercial or open-source software that can model frictional contact between such materials. Therefore, this study formulates and implements a finite element algorithm for large deformation biphasic frictional contact in the open-source finite element software FEBio. This algorithm relies on a local form of a biphasic friction model that has been previously validated against experiments, and implements the model into our recently-developed surface-to-surface (STS) contact algorithm. Contact constraints, including those specific to pressurized porous media, are enforced with the penalty method regularized with an active–passive augmented Lagrangian scheme. Numerical difficulties specific to challenging finite deformation biphasic contact problems are overcome with novel smoothing schemes for fluid pressures and Lagrange multipliers. Implementation accuracy is verified against semi-analytical solutions for biphasic frictional contact, with extensive validation performed using canonical cartilage friction experiments from prior literature. Essential details of the formulation are provided in this paper, and the source code of this biphasic frictional contact algorithm is made available to the general public.

## 1 Introduction

Frictional contact between soft hydrated biological tissues is of critical importance in biomechanics, and has particular relevance in understanding mechanics of diarthrodial joints, where opposing articular cartilage surfaces undergo reciprocal contact loading for thousands of cycles per day. Articular cartilage forms an extremely low friction bearing surface covering the ends of bones in articulating joints, which reduces wear and minimizes the energy cost of motion. As reviewed previously [1], cartilage friction is known to be dependent on the load sharing between fluid and solid constituents, leading to a distinction between the drained friction coefficient $μeq$, which is a boundary friction property of the load-bearing surfaces of the solid matrix alone, and the experimentally recorded effective friction coefficient $μeff$, which is calculated as the ratio of tangential to normal forces and is inclusive of load-sharing mechanisms. In the cartilage literature as well as our treatment here, the term “friction coefficient” refers to $μeff$ unless otherwise specified. Due to interactions between fluid and solid phases, time- and load-dependent effects produce values of $μeff$ that can vary by two orders of magnitude over relatively short periods of time. The subject of intense study over several decades, it is now clearly understood that the mechanisms by which cartilage maintains low friction derive from its fluid-filled porous fibrillar structure.

Based on the load-sharing concept initially proposed by McCutchen [2,3] and supported by the experiments of Forster and Fisher [4], Ateshian et al. [5,6] developed a nonlocal boundary friction model which described the dependence of the measured friction coefficient $μeff$ on the interstitial fluid pressurization (IFP); subsequent experimental studies confirmed the validity of this model, which was able to predict observed frictional behaviors not seen in traditional engineering materials [712]. Due to the low permeability of the porous matrix of cartilage, upon loading the interstitial fluid pressurizes for a substantial duration. This pressurized fluid bears a significant fraction of the applied load, thereby greatly reducing the contribution of contact forces transmitted by the porous solid matrix of the opposing cartilage surfaces. As viscous shear stresses are comparatively negligible in thin films of sheared synovial fluid or in the interstitial fluid of cartilage, frictional forces arise primarily from those solid-on-solid contact interactions. As a result, friction remains low when IFP is high. Over time, fluid leaves the loaded regions of the tissue as they consolidate, thus increasing the load supported by the solid matrix and consequently the friction due to solid-on-solid contacts. This load transfer explains the large temporal increase in cartilage friction seen under certain testing configurations [24]. It was subsequently shown that migrating contact areas (MCA; e.g., glass lens on cartilage strip) maintained low friction over lengthy tests [10,13], consistent with theoretical predictions of sustained IFP [14]. In contrast, stationary contact areas (SCA; e.g., glass slide on cartilage plug) rapidly depleted fluid pressurization and led to 50-fold increases in friction within an hour.

In musculoskeletal biomechanics, it is common to examine the interactions between tissues and other structures using the finite element method. However, to the best of our knowledge, there have been no efforts at creating finite element algorithms for modeling frictional contact between biphasic materials, despite the importance of frictional effects in joint mechanics. Algorithms for frictionless biphasic contact have been developed [1518], including our own prior work [19]. Commercially, ABAQUS (Dassault Systemes) can handle frictionless poroelastic (biphasic) contact, although a user-written subroutine is needed to correctly apply boundary conditions that enforce continuity of fluid pressure and normal fluid flux across the contact interface [20,21]. Finite element algorithms for frictional contact of single-phase (elastic or inelastic solid) materials have been well documented [2224]. However, the coupling between fluid pressure and frictional tractions is not trivial; a well-posed set of interfacial boundary conditions must be formulated to satisfy mass, momentum, and energy balance across the interface. Existing algorithms for single-phase friction do not provide the complete set needed for modeling friction between porous-permeable hydrated tissues. The inability to model frictional contact between such materials represents a major obstacle in computational biomechanics [25].

This study describes the formulation of a frictional contact algorithm for biphasic materials, and its implementation into the free, open-source finite element code FEBio [26]; this software has been developed specifically for the biomechanics and biophysics communities, and the inclusion of biphasic frictional algorithms in FEBio represents a significant advancement in computational modeling of biological tissues. To handle the complexities of biological contact, we combine our recently developed nonsymmetric STS algorithm for frictional contact of single-phase materials [27] with the previously validated boundary friction theory for biphasic friction [6]. This STS algorithm retains the mathematical simplicity of classical node-to-segment (NTS) formulations yet is able to solve challenging problems typically requiring the use of mortar or contact smoothing methods, allowing us to accurately simulate complex contact conditions without associated prohibitive computational expenses.

The structure of this paper is as follows: Section 2 develops a local form of the boundary friction model of Ateshian et al. [6], anticipating the need for pointwise equations in a finite element framework. Section 3 derives the continuum formulation of the finite element algorithm for frictional biphasic contact. The full discretization and linearization of the friction model can be found in the supplementary materials on the ASME Digital Collection, due to its length. Section 4 presents a validation of the biphasic friction model against classical experimental results from the cartilage friction literature, and Sec. 5 closes the paper with a discussion. A brief overview of contact kinematics from our previous work [27] is provided in Appendix  A and a summary of the constitutive models employed in this study is provided in Appendix  B.

## 2 Boundary Friction Formulation

The theoretical formulation for boundary friction in biphasic materials summarized in this section is drawn from Refs. [1], [5], and [6]; a similar concept was proposed earlier by Forster and Fisher [4]. Here, we summarize the salient equations and details of this model. The biphasic boundary friction model quantifies the load sharing between the fluid and solid constituents at the contact interface. The model considers the IFP and the fraction of the area over which the opposing porous solid matrices are in contact, denoted by $φ$. The total stress T in a biphasic material is given by
$T=−pI+Te$
(2.1)

where $Te$ is the effective (usually elastic) stress resulting from deformation of the solid matrix, p is the interstitial fluid pressure, and I is the identity tensor [28]. The normal load transmitted across the interface is $W=∫An·T·n dA$, where n is the unit outward normal at the interface and A is the area over which the surfaces are in contact. Integrating the fluid pressure p over the contact interface produces a total fluid load $Wp=−∫Ap dA$, and the ratio $Wp/W$ is defined as the fraction of interstitial fluid load support [1]. However, the actual load supported by the fluid varies due to the porous nature of the contacting surfaces, and can be written as $(1−φ)Wp$, where $1−φ$ is the fraction of the apparent contact area where pressurized fluid is present and able to support load [29]. It then follows that the total component of the load supported by frictional solid-on-solid contact is $Wss=W−(1−φ)Wp$. As the fluid pressure changes over time, all load sharing components may experience significant temporal changes. When the pressure subsides, $Wp=0$, we find that Wss becomes equal to the total applied load W.

The friction model assumes the frictional force F at the interface results primarily from solid-on-solid contact, and thus is directly proportional to Wss,
$F=μeqWss=μeq[W−(1−φ)Wp]$
(2.2)

where $μeq$ is the equilibrium friction coefficient achieved when $Wp=0$. For two idealized perfectly smooth surfaces in contact, $φ$ is the product of the solid area fractions (which are equal to the solid volume fractions by Delesse's law) of the two surfaces, e.g., $φ=1$ in the limiting case of two perfectly smooth nonporous surfaces. For two nonsmooth porous surfaces being slid against one another, as in contact between two cartilage layers, the product of the solid volume fractions only represents an idealized upper bound for $φ$. Thus, for nonsmooth surfaces contacting at asperities, or in the case of fluid film lubrication or trapped lubricant pools, the contact area fraction will be less than the product of the two solid area fractions [29]. Consequently, $φ$ must be obtained experimentally, along with $μeq$ [6]; both parameters will therefore need to be specified by the user in computational models (Sec. 4.1).

Given the relation of Eq. (2.2), the model for the experimentally measured effective friction coefficient F/W becomes
$μeff=μeq[1−(1−φ)WpW]$
(2.3)

and it is apparent that in this formulation, $μeff$ is a linear function of the interstitial fluid load support, reaching a minimum value of $μeff=φμeq$ when $Wp/W=1$. As this model neglects any contribution from viscous fluid shear, when the surfaces are completely separated by a fluid film ($Wp/W=1, φ=0)$ the friction coefficient reduces to $μeff=0$. Furthermore, in the idealized case of two perfectly smooth surfaces in contact ($φ=1$), the model of Eq. (2.3) correctly reduces to Coulomb's law and produces $μeff=μeq$.

The original development of this friction model was accomplished by integrating tractions over the entire contact surface, producing a nonlocal theory. Anticipating the implementation of this friction theory in a finite element code necessitates localizing this model to produce the pointwise counterpart. Locally, the fluid load support becomes
$WpW→−ptn$
(2.4)
where $tn=n·T·n$ is the normal component of the contact traction, and the negative sign arises from the convention in mechanics that a positive fluid pressure represents a compressive stress. A local form of the model in Eq. (2.3) is then written as
$μloc=μeq[1+(1−φ)ptn]$
(2.5)

where the negative sign associated with the pressure has been distributed. It must be noted here that there is no theoretical principle that limits the local fluid load support to the range $−p/tn≤1$. However, on thermodynamic grounds, it is impossible to have a negative friction coefficient $μloc$ in a continuum framework. Therefore, the theoretical upper bound on $−p/tn$ is $(1−φ)−1$. It must be emphasized that $μloc$ is generally not available experimentally.

The theory derived above represents a boundary contact model as reviewed in Ref. [1], as it precludes the existence of a full-thickness lubricant film that completely separates the load-bearing surfaces. This framework also accommodates mixed lubrication (a mix of boundary contact of surface asperities and pressurized fluid film pools) via the parameter $φ$ as explained in greater detail in Ref. [29]. Mixed lubrication occurs when the film thickness cannot sufficiently exceed surface roughness. During mixed lubrication, the pressure in lubricant pools becomes essentially equal to the biphasic contact pressure achieved under boundary contact, as demonstrated using our recent formulation of biphasic-fluid-structure interactions [30] in a case study posted online.1 Further discussion along these lines is deferred to Sec. 5.

## 3 Continuum Formulation

The purpose of this treatment is to consider contact of hydrated materials. A frictional algorithm for biphasic contact must address load sharing between the fluid and solid constituents and the effect of IFP on the contact tractions, while enforcing continuity of the normal component of the fluid flux.

In a series of previous papers, we have developed finite element algorithms for frictionless biphasic [19], biphasic-solute [31], and multiphasic [32] contact, although these contact algorithms were derived using certain simplifying assumptions which we choose to no longer adopt here, to produce better numerical convergence. Consequently, parts of the present formulation use notation and terminology similar to that in these referenced works, but a direct comparison should reveal certain differences. The formulation of the frictional contact problem comes directly from our previous work on frictional contact of elastic materials [27], and the present work directly adopts the notation and definitions therein. We refer the reader to that work for a comprehensive discussion of frictional contact and the finite element method. We briefly review salient contact kinematics and notation in Appendix  A.

### 3.1 Biphasic Frictional Contact.

The problem to be considered in this section consists of frictional contact between the surfaces of biphasic bodies, denoted as primary and secondary surfaces. Quantities on the primary and secondary surfaces are denoted with indices (1) and (2), respectively. Greek indices denote parametric surface coordinates and thus vary implicitly from 1 to 2, with repeated Greek indices implying summation. Following standard finite element contact procedures [27,33], the external virtual work $δGc$ arising from biphasic contact can be written as an integral over the primary surface $γ(1)$ only, yielding [19]
$δGc=∫γ(1)((δv(1)−δv(2))·t(1)+(δp(1)−δp(2))wn)da(1)$
(3.1)
where $t(1)$ is the contact traction on the primary surface, which is equal and opposite to that on the secondary surface, $t(1)=−t(2)$; $wn≡wn(1)$ is the outward normal component of the fluid flux $w(1)$ from the primary surface, $δv(i)$ are virtual velocities, $δp(i)$ are virtual fluid pressures, and $da(1)$ is an elemental area on the primary surface $γ(1)$. The contact integral is then written over the invariant parametric space of $γ(1)$, which is denoted by $Γη(1)$ [33], facilitating its linearization as required for use with an iterative solution method such as Newton's method. The elemental area is given by
$da(1)=Jη(1)dη(1)1dη(1)2, Jη(1)=|g1(1)×g2(1)|$
(3.2)
using the covariant basis vectors
$gα(i)=∂x(i)∂η(i)α$
(3.3)
where $x(i)(η(i)α,t)$ is the spatial representation of surface $γ(i)$ as it deforms in time t, in terms of contravariant surface parametric coordinates $η(i)α$. Casting Eq. (3.1) into convenient matrix notation and switching the domain of integration to the parametric frame yields the invariant biphasic contact integral
$δGc=∫Γη(1)([δv(1)δp(1)]·[t(1)wn]+[δv(2)δp(2)]·[−t(1)−wn])Jη(1)dη(1)1dη(1)2$
(3.4)
Since the parametric frame $Γη(1)$ is also a valid invariant material frame, the linearization of $δGc$ can be accomplished by applying the directional derivative operator directly to the integrand,
$DδGc=∫Γη(1)D([δv(1)δp(1)]·[t(1)wn]Jη(1)+[δv(2)δp(2)]·[−t(1)−wn]Jη(1))dη(1)1dη(1)2$
(3.5)
where it is understood that for any function f in this biphasic analysis,
$Df≡∑i=12(Df[Δu(i)]+Df[Δp(i)])$
(3.6)

To proceed with this linearization, it is necessary to provide expressions for the contact traction $t(1)$ that can differentiate between frictional stick and slip, using the kinematic framework laid out in our previous work [27]. Following the development of biphasic frictional contact boundary conditions and algorithmic details in Sec. 3.2, the full discretization and linearization of the governing equations, leading to expressions for the residual vectors and stiffness matrices, may be found in the supplementary materials on the ASME Digital Collection.

### 3.2 Biphasic Friction Formulation.

An examination of Eq. (3.4) reveals that a contact framework must provide expressions for both the normal fluid flux wn and the contact traction $t(1)$. Expressions for the normal fluid flux were presented previously [19], but will be repeated below briefly for completeness. As for the contact traction, this work implements our previously validated model of biphasic boundary friction, summarized above in Sec. 2, into a Coulomb-like framework for frictional contact [27]. In this framework, the relationship between sticking and slipping, and thus the contact traction $t(i)$ on the opposing surfaces, is described by a slip criterion $Ψ$, where on the primary surface
$Ψ=‖tT(1)‖−μloc|tn|$
(3.7)
Here, $tn=t(1)·n(1)$ is the normal component of the contact traction (negative in compression), $tT(1)=t(1)−tnn(1)$ is the tangential component of the contact traction whose scalar magnitude is $‖tT(1)‖=(tT(1)·tT(1))1/2$, and the effective friction coefficient $μloc$ is given by Eq. (2.5) with no distinction made between static and kinetic coefficients of friction. We are able to adopt the absolute value $|tn|$ as the friction multiplier in the slip criterion since $μloc$ in Eq. (2.5) has already incorporated the sign of both tn and p. The value of the slip criterion determines the stick-slip status via
$Ψ{<0sticking=0slipping$
(3.8)

Following our prior study [27], this work treats stick and slip separately, controlled by an exact return mapping predicated on the value of the slip criterion. The return mapping defines a rule for correcting a calculated stick traction which exceeds the slip limit and is thus not permissible. Stick is treated as a special case of a tied biphasic interface [32], whereas in slip, the traction is directly prescribed as a natural boundary condition. The formulation of biphasic frictional contact is presented for both penalty and augmented Lagrangian regularization schemes. As the basic framework adhered to in this paper was developed previously [27], only the necessary equations and descriptions are given below. We refer the reader to that work for a complete development of the frictional contact framework.

#### 3.2.1 Penalty Scheme.

The defining characteristic of frictional stick is a lack of relative motion between points which were previously in contact. Consequently, the stick traction is obtained by penalizing relative motion between such points,
$t(1)=εgs$
(3.9)

where ε is the contact penalty parameter and we have utilized Eq. (A5) to define the gap vector $gs$ during stick.

During slip, the normal component of the contact traction is first calculated by penalizing the normal component g of the gap, given by Eq. (A3),
$tn=εg$
(3.10)
The total traction vector in slip is then directly prescribed as
$t(1)=tn(n(1)+μlocs(1))$
(3.11)

where $s(1)$ is the unit vector in the slip direction, given by Eq. (A12). We note that Eq. (3.11) has the same form as the classical Coulomb friction considered in our previous work [27], with the standard Coulomb friction coefficient replaced by $μloc$, as evaluated from Eq. (2.5). A trial state and return map, adapted from our previous work [27] and presented in Sec. 3.2.3, is employed to differentiate between stick and slip.

As discussed in Sec. 2, the jump conditions for a biphasic mixture require continuity of the fluid pressure across a contacting interface and therefore an expression for the normal fluid flux wn can be obtained by penalizing the fluid pressure gap π between contacting points [19] when contact occurs $(tn<0)$. When there is no contact $(tn=0)$ the surfaces have pulled apart and the essential boundary conditions on fluid pressures $p(i)$ must be enforced. The corresponding boundary equations are then
$wn=εpπ, tn<0p(i)=0, tn=0$
(3.12)
where $p(i)=0$ prescribes zero fluid pressure (free-draining conditions) on the portions of the boundary where no contact takes place. We define the fluid pressure gap as
$π=p(1)−p(2)$
(3.13)

and $εp$ is the pressure penalty parameter, which has units of hydraulic permeability per unit length (e.g., m3/N$·$s, similar to the hydraulic permeability of a membrane). Equation (3.12) is valid for both stick and slip. Note that $p(1)=p(1)(η(1)α)$ is evaluated at a point with parametric coordinates $η(1)α$ on the primary surface, whereas the pressure on the secondary surface, $p(2)=p(2)(η(2)α,t)$, is evaluated at the parametric coordinates of the intersection of a ray issued from that primary surface point and normal to $γ(1)$. As detailed previously [27] and summarized in Appendix  A, the parametric coordinates of intersection are dependent on the stick-slip status, so care must be taken to evaluate wn from Eq. (3.12) using the contact kinematics determined by Eq. (3.8). In practice, this is accomplished by calculating wn once the stick-slip status has been resolved for each iteration.

#### 3.2.2 Augmented Lagrangian Scheme.

The augmented Lagrangian scheme presented herein was developed in our previous work [27] as a modification of the approach proposed by Simo and Laursen [34]. Briefly, this is a first-order augmentation scheme that utilizes Uzawa's algorithm [35], where the multipliers are updated outside of the Newton step, producing a double loop algorithm [23,24] and preserving quadratic convergence of Newton's method near solution points.

During stick, the traction is calculated by augmenting the vector gap $gs$,
$t(1)=λs+εgs$
(3.14)
where $λs$ is the vectorial Lagrange multiplier in stick. In slip, the normal component of the contact traction is first calculated by augmenting the normal gap g,
$tn=λn+εg$
(3.15)
where λn is the normal Lagrange multiplier. The total traction vector in slip is then directly prescribed as
$t(1)=tn(n(1)+μlocs(1))$
(3.16)

where Eq. (3.15) has been used. The update formulas for the Lagrange multipliers, as well as a detailed discussion on the benefits of the active-passive multiplier strategy employed here, can be found in our previous work [27]. Here, it suffices to note that by augmenting only the normal gap and employing Eq. (3.16), we have ensured 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, presented in Sec. 3.2.3 and controlled by the slip criterion, is used to differentiate between stick and slip.

In this augmentation scheme, the normal fluid flux is given by
$wn=λp+εpπ ,tn<0p(i)=0 ,tn=0$
(3.17)
where λp is the fluid pressure Lagrange multiplier, and the update formula for λp has been detailed in our prior work on frictionless biphasic contact [19] and is given by
$λp←λp+εpπ$
(3.18)

Unlike $λs$ and λn [27], λp is not dependent on the stick-slip status. However, as noted before, Eq. (3.17) must be evaluated after the stick-slip status has been resolved.

As detailed above, augmented Lagrange iterations calculate a Lagrange multiplier at every integration point, which is used to enforce certain contact criteria to a user-specified tolerance as the augmented Lagrangian contact model is solved pointwise at integration points. This calculation occurs independently at each integration point on the primary surface, or on both surfaces if a two-pass algorithm is used, and thus enforces no sense of smoothness or continuity between adjacent surface elements. The penalty method does allow for a certain smoothing which is often desirable, as we detailed in our prior work [27]. To avoid this limitation of the augmented Lagrange method, we developed a smooth augmentation method, which has been implemented as the “smooth_aug” contact flag in FEBio. This method takes advantage of the tendency of solid materials to act as a low-pass filter, which manifests in a finite element framework as the effect where element-level stresses, evaluated at integration points, filter out the spatial noise of surface-level contact stresses. This becomes most evident in a scenario where a coarsely discretized surface with poor contact enforcement may display noisy contact tractions and forces, yet produce far smoother element-level stresses. This one-way, low-pass filtering effect implies that, in a numerical contact algorithm, which relies on augmentation, there is no unique solution for the augmented contact traction as long as it satisfies the force equilibrium embodied in the contact integral, leading to potentially noisy contact traction solutions. Instead, the smooth augmentation method sets the Lagrange multipliers in Eqs. (3.14) and (3.15) equal to the surface traction and surface contact pressure evaluated by projecting the stress tensor from integration points in the solid element attached to the surface, thereby replacing the standard update formulas presented in Ref. [27]. When adopting smooth augmentation, the tractions are calculated according to Eqs. (3.14)(3.16) as appropriate; the only difference is that the update formulas in Ref. [27] are replaced by the following. If the slip criterion $Ψ<0$ (stick), we update the active multiplier $λs$ and derive the passive multiplier λn according to
$λs←Tel·n(1)λn=λs·n(1)$
(3.19)
Alternatively, if $Ψ=0$ (slip), we update the active multiplier λn and derive the passive multiplier $λs$
$λn←n(1)·Tel·n(1)λs=λn(n(1)+μlocs(1))$
(3.20)

In the above formulas, $Tel$ is the total Cauchy stress in the solid element attached to the contact face2 and $n(1)$ has the usual meaning of the outward unit normal at the current integration point.

Assigning these smoothed values to the Lagrange multipliers prevents solutions at integration points from evolving independently of neighboring values and restores a smoothness typically absent from augmented Lagrange methods. More concisely, the faceted nature of finite element surfaces manifests as noise in the contact traction when an augmentation scheme enforces contact constraints to a high degree of accuracy; this enforcement emphasizes the nonphysical discretization of a physically smooth and continuous surface. Rather than apply more expensive smoothing algorithms to the faceted surface, we developed this method to smooth the Lagrange multipliers acting on the surface instead.

#### 3.2.3 Stick-Slip Algorithm.

Determining whether stick or slip is occurring at a given instant is accomplished by a trial state and return map, and has the same form for both penalty and augmented Lagrangian regularization schemes [27]. We initially assume stick and calculate a trial traction $t̃(1)$, utilizing either Eq. (3.9) or Eq. (3.14). The trial normal and tangential components $t̃n$ and $t̃T(1)$ are evaluated from $t̃(1)$ and inserted into the slip criterion $Ψ$
$Ψ=t~T(1)−μloc|t~n|$
(3.21)
Based on the slip criterion and trial traction vector, the traction vector is calculated from the return mapping as
$t(1)={t̃(1)Ψ<0, stickingtn(n(1)+μlocs(1))Ψ=0, slipping$
(3.22)

where tn is given by either Eq. (3.10) or Eq. (3.15).

## 4 Verifications and Validations

### 4.1 Finite Element Implementation.

The presented algorithm was implemented in the custom, open-source finite element code FEBio as “sliding-biphasic” contact. All finite element analyses were performed with three-dimensional models and eight-node hexahedral isoparametric brick elements. As a consequence, contact surfaces consisted of four-node quadrilateral isoparametric element faces. Per the discussion in our previous works [19,27], Gaussian quadrature was used for integration, with four integration points per element face. To assist in calculating the penalty parameters, an autopenalty feature was added, which uses geometrical and material properties to calculate a suitable penalty parameter. In this feature, the user's input parameter acts as a scale factor. The autopenalties are given by
$ε=α1N∑i=1NEiAiViεp=αp1N∑i=1NkiAiVi$
(4.1)

where the summation is taken over all N element faces representing the surfaces in the contact pair, Vi is the volume of each element, Ai is the area of the element face on the contact surface, Ei is a measure of the average Young's modulus, and ki is a measure of the average hydraulic permeability, in the direction normal to the surface. The parameters α and αp are user-defined scale factors. When the autopenalty feature is not used, the penalties reduce to these user-defined parameters, $ε=α$ and $εp=αp$. For elastic (single-phase) friction, the autopenalties are calculated once during model initiation. For hydrated media, however, much of the effective stiffness arises from the pressurized interstitial fluid, and as the tissue depressurizes the material becomes softer. As this may result in the penalty becoming too severe with fluid pressure loss and concomitant compaction, an “update_penalty” contact flag was added, which updates the autopenalties at the beginning of each time-step. In all examples in the following sections, the autopenalty feature is used. Since a rigid body cannot be the primary surface, only the MCA example of Sec. 4.6 uses a two-pass contact algorithm [19,27]. To avoid numerical error at the first time-step, contacting surfaces were positioned to have a minor overlap at the beginning of the simulation in all force-driven analyses.

The primary user-supplied parameters for biphasic friction analyses are the equilibrium friction coefficient $μeq$ and the solid–solid contact area fraction $φ$. For users unfamiliar with these concepts, we briefly remark upon best practices for determining the correct values to apply. Experimentally, $μeq$ may only be measured in the context of SCA friction, where typically one counterface is impermeable; once measured, however, values of $μeq$ may be applied to any contact modality, as shown by experimental evidence and the validations below. The bulk of experimental evidence supports approximately $0.15≤μeq≤0.30$ for articular cartilage, though we emphasize that there are no restrictions upon this parameter beyond $μeq>0$. Similarly, $φ$ may be approximated as the product of the solidities of the contacting surfaces, but will generally be lower due to asperity contact and the potential presence of fluid film lubrication or trapped lubricant pools. For typical cartilage samples, where the solid volume fraction is often no more than 0.25, $φ∽0.1$ for contact between cartilage and an impermeable surface, and $φ∽0.01$ for cartilage-on-cartilage contact. As it represents an area fraction, the only restriction is $0≤φ≤1$.

Modeling challenging friction experiments computationally necessitates a discussion of the best ways to extract physically relevant parameters. In FEBio, contact forces are available at a contact surface, calculated from integrating contact tractions. In much the same manner, the fluid force Wp is available and obtained by integrating fluid pressures. These contact forces are identical to reaction forces, such as those on rigid bodies interfacing with the cartilage. Experimental friction measurements are generally obtained by taking the ratio of tangential to normal forces, using values directly supplied by or calculated from load cell data. In a similar manner, all computational friction values reported in this study are obtained from the ratio of the surface tangential force to normal force,
$μ=F/W$
(4.2)

Location-dependent friction coefficients $μloc$ (and fluid load support $−p/tn$) may of course be found by using contact tractions and fluid pressure, and the contact variables “effective friction coefficient” and “local fluid load support” in FEBio provide those measures. However, comparison to experimental data is best achieved by applying experimental data collection practices to the computational model. Following this line of thinking, finite element results for the calculated friction coefficient employ Eq. (4.2) for the purpose of comparison with experimental data, whereas predictions of the friction coefficient using our friction model employ $μeff$ in Eq. (2.3), using the interstitial fluid load support fraction ($Wp/W$) calculated from the finite element analysis. It should be understood that the finite element results internally rely on using the local value $μloc$ in Eq. (2.5) in the frictional biphasic contact algorithm. Using the model for $μeff$ to compare against experimental μ is also in agreement with experimental protocols that measure fluid load support through pressure transducers that report surface-averaged results [7,8,36,37]. To ensure the most accurate fluid pressure solution, mesh biasing is employed for all models near boundaries where fluid pressure may change rapidly.

Reporting the results of a study also requires additional care. In reciprocal sliding experiments, a transition occurs between slip and stick in the vicinity of the reciprocating motion; during this stick period, the tangential load response represents the reciprocal shearing of the loaded bearing materials, whose temporal profile is not relevant to the determination of the sliding friction response. In practice, experimentalists either report trimmed data, with a region of the wear track removed, or cycle-averaged data, which has the additional benefit of correcting for curvature errors in sample alignment [38]. Similarly, in all finite element analyses reported below, data from 10% of the travel range from each side of a reciprocal sliding cycle have been removed in accordance with experimental practice. For both SCA under a constant normal load and MCA configurations below, friction coefficients, displacements, and fluid load support fractions are reported on a cycle-averaged basis. For SCA under a cyclical normal load, trimmed values are reported. When direct comparisons are made to experimental data in Sec. 4.5, the raw experimental data have been reprocessed in the same manner.

### 4.2 Smoothing Fluid Load Support.

The presented frictional algorithm critically depends on the local fluid load support, originally evaluated as $−p/tn$ at the contact surface integration points. However, this evaluation at the surface is susceptible to numerical noise for the reasons given above for tn, especially when using the augmented Lagrange method. This spatial noise becomes problematic when evaluating $μloc$, especially when a noisy $−p/tn$ exceeds the theoretical upper bound $(1−φ)−1$; in those cases, $μloc$ must be set to zero at those contact surface integration points. To help resolve this issue in practice, we implemented a procedure which calculates the fluid load support used to find $μloc$ in a smoothed manner by projecting element-level quantities to the surface, much in the same manner that Sec. 3.2.2 introduced a smoothed form of Lagrange augmentation. It is important to note that the nodal fluid pressures and contact tractions at integration points calculated by the finite element solver remain unchanged; the current discussion is restricted only to the local fluid load support ratio used in finding $μloc$, which was specified by Eq. (2.4) as $−p/tn$. To introduce the necessary smoothing, we average the integration point stresses and fluid pressures in the solid element attached to the current contact face, then evaluate the local fluid load support at each integration point as $−pel/(n(1)·Tel·n(1))$, where $Tel$ is the element-averaged total Cauchy stress and pel is the element-averaged fluid pressure. Investigating this method showed that the frictional behavior and contact force magnitudes did not change in any substantial way beyond presenting a cleaner response. This is to be expected, as best practices for simulating biphasic materials call for resolving boundary layers at surfaces where fluid pressures may change rapidly using mesh biasing; the thin elements at the surface prevent the element-level fluid pressure from differing significantly from that at the surface. In FEBio, we have made this internal calculation dependent on the contact flag “smooth_fls.” As this is an internal change, for simplicity we will continue to refer to –p and tn throughout this paper, with the understanding that $μloc$ has been calculated by turning on this flag.

### 4.3 Guidelines for Selecting Constitutive Models for Articular Cartilage.

Several classical cartilage experiments from the literature, demonstrating characteristic biphasic frictional behavior, are illustrated below to validate the accuracy of the implementation. Due to the nature of porous-permeable friction, the development of IFP in the cartilage matrix is critical to obtain an accurate friction response. As noted by both experimental [37,39] and theoretical [4042] studies, the magnitude of IFP in porous soft tissues is directly coupled to the tension-compression nonlinearity arising from the fibrous nature of the solid matrix. Therefore, in order to conclusively validate the frictional algorithm, we must utilize material models for cartilage which are able to develop physiologically realistic IFP, or else the magnitude of the computational results would not match the experimental frictional data. For example, a linear isotropic solid matrix would not be an appropriate choice for modeling the porous matrix of cartilage, since it can develop a maximum IFP of only 33% in unconfined compression [43], whereas the extremely low friction seen in cartilage under MCA conditions is due to an initial IFP that approaches 100% as verified experimentally [8,37].

When adding fiber reinforcement to cartilage models, a potential issue may arise in calculating the fluid load support. The normal component $tn=n(1)·T·n(1)$ of the contact traction may be expressed as $tn=−p+tne$ when using Eq. (2.1), where $tne=n(1)·Te·n(1)$ is the contribution to the total traction resulting from solid matrix stresses. The fibrillar solid matrix of cartilage has a complex ultrastructure, and continuous fiber distributions have been shown to provide increased accuracy in cartilage modeling, even when using a simplified isotropic fiber distribution [44]. However, for finite elements representing the articular surface zone, where actual collagen fibrils are oriented tangential to the articular surface, it is important to reproduce this tangential ultrastructure in the model. Otherwise, a continuum representation of nontangential fibrils, which can only sustain tension, could produce a tensile traction $tne>0$ in the solid matrix, which combines with –p to produce an unacceptable $−p/tn>(1−φ)−1$.3 This nonsensical construct would then produce an interstitial fluid load support $−p/tn$ which is greater than its theoretical upper bound due to the artifact of modeling nontangential fibers at the surface, compromising the validity of the friction model in Eq. (2.5). Therefore, it is important to adopt fibril-reinforced constitutive models for articular cartilage that satisfy this constraint of tangential fibril ultrastructure in the superficial zone. This constraint was not adopted in a previous study from our group, where we utilized multilayered models of knee articular cartilage to account for depth-dependent inhomogeneity, while using an isotropic fibril distribution for all layers [48], though the finite element mesh in that study was sufficiently coarse to exhibit an IFP smaller than 100% at the articular surface.

There are distinctions between modeling SCA and MCA contact, with SCA contact much more forgiving of nontangential fiber distributions. In standard plug-on-flat SCA sliding, the entire cartilage surface is in contact with the flat counterface at all times, and the tangential normal strains are much smaller in magnitude than the axial normal compressive strain, producing a very narrow range of angles out of the tangential plane over which fibers are in tension [45]; those circumstances are less likely to prevail over the entire contact interface for nonconforming MCA contact, thus requiring a strict usage of only tangential fiber distributions in the element layer underlying each contact surface. Even so, in MCA, one may still encounter fluid load support fractions $>(1−φ)−1$ due to the presence of nontangential fibers at greater depths from the contact surface, which propagate these elevated values of $−p/tn$ to the element layer underlying the contact surface. Fortunately, we were able to eliminate this issue by adopting a physiologically realistic modeling assumption [49]: we let the fiber toe region in the deeper zones be longer than that in the superficial zone, even if only slightly. Since the magnitude of $−p/tn$ depends critically on the ratio of elastic moduli of fibers and ground matrix [37,41], a longer toe region prevents the tensile modulus in these deeper zone fibers from contributing as significantly as the superficial zone fibers to $−p/tn$. In this study, we assumed a piecewise linear modulus transition in the superficial zone, and then modeled a toe region throughout the remaining depth.

The constitutive models selected in the following validations are chosen to develop appropriate levels of IFP and, when possible, to match experimental displacement data such as creep deformation. The results show several valid methods of constructing constitutive models that avoid the pitfalls described above. However, the primary purpose of this work is to verify and validate a computational model for biphasic friction and not to demonstrate the state of the art in modeling articular cartilage.

### 4.4 Verification of Stationary Contact Areas Friction.

Verification was performed against the work by Krishnan et al. [9], who conducted a theoretical analysis to determine $μeff$ when an impermeable glass slide was reciprocated against cartilage in SCA conditions under the influence of both constant and cyclical applied loads. In Sec. 4.5, we will also validate our implementation against experimental results from that same study, and defer more detailed discussion of the meaning of the observed behaviors until then.

In these two verification examples, dimensions, constitutive models, and parameters were adopted from the original analysis [9]. The cartilage sample was a cylindrical plug with a diameter of 4.8 mm and a thickness of 2.1 mm. The constitutive model for the cartilage was a biphasic material with constant permeability and a solid matrix represented by a conewise linear elastic (CLE) model with cubic symmetry [41]. The constant permeability was set at k =0.0025 mm$4/$N$·$s, with CLE parameters4$λ+1=12.86$ MPa, $λ−1=0.3$ MPa, $λ2=0.48$ MPa, and $μ=0.17$ MPa. Here, $λ+1, λ−1$, and λ2 represent tensile, compressive, and off-diagonal first Lamé coefficients, and μ is the second Lamé coefficient. The friction parameters were given by $μeq=0.153$ and $φ=0.1$. In each example, loading was applied to the glass slide while it reciprocated against cartilage over a 9 mm range (±4.5 mm in each half-cycle) at a speed of 1 mm/s.

In the SCA model with a constant load, a load of 0.00134 N was ramped up linearly in 1 s and then held constant until 2500 s. The penalty method was used for contact enforcement, with penalty scale factor $α=0.5$. No pressure penalty was needed, as contact with an impermeable interface precluded the need to enforce fluid pressure continuity across an interface according to Eq. (3.12). For the case of SCA under the influence of a cyclical load, the load was ramped up to 0.00134 N linearly in 0.1 s, and then a sinusoidal load with amplitude 0.00089 N was applied at a frequency of 0.05 Hz, such that the average load magnitude was equal to that in the constant load case. Here, we note that these load magnitudes are smaller than those used in the Krishnan et al. [9] analysis by a factor of 1000, to match their study's theoretical assumption of infinitesimal strains given that FEBio only solves for finite deformations. Finally, due to symmetry considerations, only half of the geometry was modeled, thus all loads were similarly reduced to half their given value.

The undeformed mesh used in both verification problems is shown in Fig. 1(a). Zero fluid pressure was prescribed on the cylindrical boundary of the tissue, which is exposed to ambient conditions. The bottom surface was fixed in all directions and symmetry boundary conditions were applied at the cut. As anticipated from Eq. (2.5), an exact linear relationship between the friction coefficient and fluid load support fraction $Wp/W$ emerges for both the theory and FEBio model (Fig. 1(b)); the displayed equation is a linear regression to the FEBio results with $R2=1$. The regression also verifies the implementation of Eq. (2.5), which specifies a linear relationship between the effective friction coefficient and fluid load support. It is of further interest to note the $y−$intercept given by the regression is 0.153, exactly in agreement with the prescribed value of $μeq$, which is the value of $μeff$ theoretically achieved when fluid load support becomes zero. Similarly, $μ=0.015$ when the fluid load support fraction is unity, according to the linear regression, which is the value predicted from both Eqs. (2.3) and (2.5) given the parameters $μeq$ and $φ$ employed in this problem. Plots of both the friction coefficient (Fig. 1(c)) and the fluid load support fraction (Fig. 1(d)) show exact agreement between the theoretical predictions and the FEBio model, demonstrating the ability of the algorithm to capture the link between IFP and friction.

Fig. 1
Fig. 1
Close modal

Similarly exact agreement can be seen for the case of SCA friction under the influence of a cyclical applied load, both in the friction coefficient (Figs. 2(a) and 2(c)) and the fluid load support fraction (Figs. 2(b) and 2(d)). The unusual results here, where μ can far exceed the maximum $μeq$ obtained under a constant load, are due to negative fluid pressure, as explained in Sec. 4.5.2. The relationship between μ and fluid load support remains virtually the same as in Fig. 1(b) for the constant loading case; linear regression (not shown) produced $y=−0.1376x+0.153$ with $R2=0.9999$. In both cases, the ability of the proposed algorithm to capture friction that depends on coupled fluid–solid mechanics in a manner identical to theory serves to verify our implementation.

Fig. 2
Fig. 2
Close modal

### 4.5 Validation of Stationary Contact Areas Friction.

The experimental study by Krishnan et al. [9] offers an ideal choice for validation of our finite element implementation, as the authors performed SCA sliding with both constant and cyclical applied loads on the same samples, under a challenging finite deformation protocol. For the most rigorous validation, we chose to model the geometry and material properties of an individual sample (sample 3C) and predict the resulting friction coefficient. Though two finite deformation tests were performed on sample 3C on separate days, we tried to match compressive displacement data from both tests simultaneously to obtain a single set of material parameters. To reduce computational time while seeking material parameters, 3 deg wedge models of the sample-specific height and diameter (see below) were created to take advantage of axisymmetry, with appropriately scaled applied loads. No sliding was applied, and only frictionless contact was used for fitting material properties; these models were able to run in a few minutes and provide creep displacement data which was remarkably predictive of the final three-dimensional model displacement. As seen in the following sections, the combined fitting of creep displacements produced by both constant (Fig. 4(a)) and cyclical (Figs. 5(a) and 5(b)) applied loads was able to generate material parameters which provided a suitable prediction of μ in both cases.

The cartilage plug was modeled as a cylinder with a diameter of 4.8 mm and a thickness of 2.83 mm according to the measured dimensions of sample 3C. Due to symmetry, only half of the geometry was modeled; hence, all loads reported in the following sections were halved during simulation. As in Sec. 4.4, the bottom surface of the cartilage plug was fixed, and the symmetry boundary was prevented from displacement normal to the symmetry plane. Zero fluid pressure was prescribed on the lateral cylindrical surface. Extrapolating the experimental friction data from the constant load SCA experiment (Fig. 4(c) indicates equilibrium was not attained in 2500 s) produced $μeq=0.2565$. The solid–solid contact fraction was prescribed to be $φ=0.05$ such that μ could better match experimental data when $Wp/W→1$.

The constitutive model for the cartilage itself was separated into two layers: a superficial zone comprised only of the topmost element layer and a middle/deep zone making up the remainder of the tissue. For both layers, the solid matrix was modeled as a mixture of a Holmes–Mow solid [50] and a Newtonian viscous solid, reinforced by a continuous fiber distribution. In the superficial zone, a circular distribution of fibers in a plane parallel to the articular surface was used, with individual fibers having an exponential power law strain energy density. A spherical fiber distribution was chosen for the deep zone (see the discussion in Sec. 4.3 regarding fibers not parallel to the surface). The Holmes-Mow solid parameters were E =0.65 MPa, $ν=0.3$, and $βhm=0.55$, where E is the Young's modulus, ν is the Poisson's ratio, and βhm is an exponential stiffening coefficient [51]. The fiber parameters are αf, βf, and ξ, where αf is a power-law coefficient, βf is an exponential coefficient, and ξ is related to the tensile modulus. For both layers, $αf=0$ and ξ = 40 MPa were selected. The parameter βf was chosen to be 2 in the superficial zone and 2.5 in the deep zone, where $βf=2$ ensures a piecewise linear curve when the modulus transitions from compressive to tensile and $βf>2$ ensures a continuous transition (a toe region). The inclusion of Newtonian viscosity to the solid, representing the viscoelasticity of the proteoglycans [52], provided a better fit to the amplitude of the creep displacement during cyclical loading. The parameters for this model were $κv=0 MPa·s$ and $μv=1 MPa·s$, where κv is the bulk viscosity and μv the shear viscosity of the solid. Finally, strain-dependent Holmes-Mow permeability [50] was prescribed, with $αp=0, k0=0.0045$ mm$4/$N$·$s, and m =8. Here, k0 represents the zero-strain permeability, m is an exponential stiffening coefficient, and αp controls an additional power-law dependence we chose to neglect. For convenience, these constitutive models are provided in Appendix  B and may also be found in the FEBio manual [32].

#### 4.5.1 SCA Friction Under Constant Load.

During this component of the SCA study, a constant load of 13.4 N was applied to the tissue surface. Following experimental data, the load was linearly ramped up over 17.5 s and then maintained until 2500 s elapsed. This load was applied to a glass slide, modeled as a rigid, impermeable loading platen in contact with the cartilage surface; reciprocation of the slide was prescribed over a 9 mm range (±4.5 mm in each half-cycle) at a speed of 1 mm/s. The impermeable platen was selected as the secondary contact surface, with the cartilage chosen to be the primary. The penalty method was used for contact enforcement, with α = 5. A pressure penalty was not needed in this problem due to the impermeable counterface.

Figure 3 shows the undeformed and final deformed state of the mesh. The large deformation induced by the SCA configuration arises due to complete interstitial fluid depressurization. As fluid exits the tissue, the pores collapse and the tissue volume shrinks. Without the inclusion of solid matrix viscosity, this volume reduction can lead to extreme mesh distortions that potentially prematurely terminate the finite element analysis, as briefly mentioned below and covered in Sec. 5 in more detail. The material parameters chosen were able to provide a good fit to the experimental creep displacement (Fig. 4(a), $R2=0.982$). Simulating the experiment with these material properties showed the friction coefficient μ rising toward an asymptotic limit as the fluid load support decreases, and the increase in friction mirrors the decrease in fluid load support (Figs. 4(b) and 4(c)), as seen in Sec. 4.4. Though Krishnan et al. did not measure fluid load support (Fig. 4(b)), the predicted friction coefficient μ is in good agreement with the experimental values (Fig. 4(c), $R2=0.981$). Furthermore, a closer look at the early-time response of μ in Fig. 4(d) shows a distinctive similarity between the behavior of the experimental results and computational model, where the friction coefficient initially decreases over the first 17.5 s before beginning to rise sharply, a consequence of the somewhat lengthy load application. As the load steadily increases up until t =17.5 s, interstitial fluid continues to pressurize, thus lowering the friction coefficient according to Eq. (2.5). Once the applied load reaches a steady value, relaxation begins and the fluid load support decreases, leading to the subsequent monotonic increase in friction coefficient. The ability to capture such experimental effects suggests that our contact algorithm substantially recapitulates the physics of biphasic contact.

Fig. 3
Fig. 3
Close modal
Fig. 4
Fig. 4
Close modal

In this model, the fluid load support peaked at 93.8% immediately following load application but decreased to a low of 9.2% in the final sliding cycle at t =2500 s, producing a corresponding rise in the friction coefficient from $μmin=0.028$ to $μ=0.236$. As likely occurred in the experiment, the fluid pressure has not entirely subsided by 2500 s and thus $μeq$ was not attained. Though a good prediction of μ was achieved for this sliding configuration, Krishnan et al. [9] report $μmin=0.0089$ for sample 3C, smaller than the minimum value obtained from the FEBio simulation (Fig. 4(d)). This minor discrepancy is likely due to our adoption of a fixed value for $φ$: the theoretical minimum friction coefficient $μmin$ can be calculated from Eq. (2.5) when $−p/tn=1$ as $μmin=φμeq=0.0128$, using our prescribed values for $φ$ and the experimental $μeq$; thus our choice of $φ$ is unable to reproduce the experimental $μmin=0.0089$. A smaller value of $φ$ would resolve this discrepancy, especially when accounting for the fact that $φ$ may increase temporally as pressurized fluid trapped between the glass and cartilage asperities eventually gets depleted as it flows into the cartilage, as explained in an earlier theoretical study [29].

As a modeling detail, note that the right-side corner of the cartilage in Fig. 3(b) is curled upward due to stick behavior in the stick-slip response. This elevated corner is penetrating the glass slide as the friction from the platen moving to the left drags it upward. The contact violation is minor and only occurs as the cartilage becomes depressurized and softer, due to the lack of IFP, along with the concomitant increase in friction associated with decreased IFP. Note that it is numerically advantageous not to force this corner back down. If a higher penalty or augmented Lagrange method is used to strictly enforce near-perfect contact, the hex element at the corner becomes degenerate as both the top and right sides become flush with the platen, generally leading to numerical convergence difficulties. That this happens is not an indictment of the algorithm or modeling process, but rather a consequence of assuming a plane-ended cylindrical sample for the geometry, combined with a constitutive model which may not provide sufficient shear stiffness or fiber-matrix interaction to prevent shearing at the surface when significant depressurization occurs [53,54].

#### 4.5.2 SCA Friction Under Cyclical Load.

In the cyclical SCA study, the normal loading waveform was composed of a 13.4 N steady load, superimposed with a cyclical load of amplitude 8.9 N. To apply this load in FEBio in a way that matched the experimental procedure, the normal load was ramped up to 4.5 N in 1 s, where it intersected with the trough of a sinusoidal waveform of amplitude 8.9 N centered at 13.4 N. The resulting waveform had a frequency of 0.05 Hz and peak and trough force values of 22.3 and 4.5 N, respectively. Due to symmetry, only half the geometry was modeled; thus, all referenced loads were halved for the simulation. This load was applied to a rigid body representing the impermeable glass slide. As before, the slide reciprocated over a 9 mm range, traversing 4.5 mm on each stroke before reversing direction. The rigid body was taken to be the secondary surface, with the cartilage surface the primary. Contact conditions were enforced with the penalty method and α = 5.

The simulated creep displacement was successfully fitted to the experimental response (Figs. 5(a) and 5(b), $R2=0.923$), given the constraint that a single set of material parameters had to fit both constant and cyclical loading-induced SCA creep displacements from experiments performed on the same sample on different days. From these fitted properties, the friction coefficient predicted by the present algorithm suitably replicates the experimentally-recorded μ (Figs. 5(c) and 5(d), $R2=0.977$), with both responses displaying large oscillations growing in magnitude. These results were predicted using the same $μeq=0.2565$ as obtained from the constant load experiment in Sec. 4.5.1; in that experiment, μ rose asymptotically toward $μeq$ as depressurization occurred. Here, as evident from Fig. 5(d), cyclical loading in the FEBio model was able to produce fluctuating values for μ as high as $μ=0.51$, in close agreement with experimental results (Fig. 5(c)). As explained in the study by Krishnan et al. [9], such high friction results from negative (subambient) interstitial fluid pressures, as seen in Fig. 6(a); such negative pressures cause suction between the cartilage and glass slide. The correct prediction of a friction coefficient μ greater than the constant loading maximum value $μeq$ validates our implementation against an effect unique to hydrated tissues whose frictional response is coupled to the IFP.

Fig. 5
Fig. 5
Close modal
Fig. 6
Fig. 6
Close modal

Results of Fig. 6(b) serve to further verify the implementation of the theoretical model described in Sec. 2 by comparing the computational friction coefficient μ (calculated using Eq. (4.2)) to the theoretical value $μeff$ obtained directly from Eq. (2.3) by using the computational values of the fluid load support $Wp/W$. As noted previously, data from either end of the sliding range is excluded, producing the gaps in the computational data points for μ in Fig. 6(b). The agreement between μ and $μeff$ further verifies the successful implementation of the model and demonstrates that the production of such elevated friction coefficients directly follows from Eq. (2.3).

To illustrate the development of negative fluid load support, Fig. 7 shows four snapshots of the finite element mesh during one 20 s sinusoidal loading period, displaying interstitial fluid pressure contours. When the cyclical normal load is greatest in magnitude (Fig. 7(a)), the cartilage plug achieves its smallest height and highest fluid pressure of the cycle. The load magnitude then begins to reduce (Fig. 7(b)), allowing the hyperelastic solid matrix to recover and achieve its maximum thickness for the cycle in the 10 s required for the load magnitude to reach its minimum (Fig. 7(c)). The increase in height causes collapsed pores to partially expand. However, the extremely low permeability of the cartilage prevents the expanding pores from imbibing fluid on such a rapid time scale. As a result, the fluid pressure begins to dip below ambient levels until it is strongly negative when minimum load magnitude is reached (Fig. 7(c)). The load magnitude subsequently increases, the sample repressurizes, and the next cycle begins (Fig. 7(d)). As more fluid is lost throughout the testing duration, the plug becomes further compacted. The pores then have progressively less fluid volume and hence the magnitude of the developed suction grows, as the elastic recovery of the solid matrix further deforms the pores. This explains why the magnitude of the oscillations in the creep displacement, friction coefficient, and fluid load support continue to increase as the test progresses.

Fig. 7
Fig. 7
Close modal

### 4.6 Spherical Counterface Migrating on Strip in MCA Contact.

In an earlier theoretical study [14], it was predicted that elevated IFP could be sustained indefinitely if the contact configuration produced a continuously migrating contact area on a biphasic layer, when the convective velocity of sliding was much greater than the diffusive velocity of interstitial fluid flow. These theoretical predictions were later confirmed by an experimental friction study demonstrating virtually no increase in μ under a constant load and MCA reciprocal sliding contact for an hour [10]. In this section, we validate our finite element implementation of frictional sliding biphasic contact against the MCA experiments of Caligaris and Ateshian [10], by examining the temporal dependence of the friction coefficient on loading conditions and confirming that low friction may be maintained indefinitely under a migrating contact configuration.

Reciprocal sliding contact was modeled between a cartilage strip and a cartilage plug, or a cartilage strip and a glass lens. The modeled dimensions of the cartilage strip were taken to be 20 × 5.5 × 2.5 mm. Both the cartilage plug and glass lens were modeled identically, as a 16 mm diameter core of a hollow sphere with an outer radius of 18 mm. The thickness was taken to be 2 mm for both the plug and lens. Reciprocal sliding with a range of 20 mm (±10 mm per half-cycle) was prescribed, under a normal compressive load of 6.3 N. Equilibrium friction $μeq=0.201$ was chosen based on the average of the SCA equilibrium values reported for experiments E1 and E2 in Ref. [10]. Based on Ref. [6], the solid–solid contact fractions were prescribed as $φ=0.12$ for glass–cartilage (GC) contact and $φ=0.0144$ for cartilage–cartilage (CC) contact. Augmented Lagrangian contact enforcement was used with augmentation smoothing (Sec. 3.2.2) and a convergence tolerance on the contact traction set at $Ptol=0.1$ for CC and $Ptol=0.2$ for GC (see Ref. [27] for a discussion of augmented Lagrange parameters). Penalty parameters were α = 1 and $αp=1$ for CC, while GC only required specifying α = 5. In addition, the “update_penalty” contact flag was turned on for both cases (Sec. 4.1). The glass lens was modeled as an impermeable rigid body and designated the secondary contact surface, with the biphasic strip the primary. In CC contact, a two-pass algorithm was used to improve model fidelity, thus each surface played the role of primary and secondary once. As before, symmetry considerations resulting in the applied loads being halved, as only half the problem was modeled. For both the plug and the strip, the top layer of elements underlying the contact surface was designated the superficial zone, with the remainder of the cartilage comprising the middle/deep zone. Figure 8(a) shows the mesh for this validation and indicates the cartilage zones.

Fig. 8
Fig. 8
Close modal

For both zones, the cartilage was modeled as a biphasic material with a constant isotropic permeability and a neo-Hookean solid matrix reinforced by discrete fiber bundles. Fiber reinforcement in the superficial zone consisted of two orthogonal fiber bundles parallel to the articular surface with an exponential power law strain energy density (Appendix  B), and in the deep zone a third bundle was added normal to the articular surface. Material parameters were taken from a prior study [55] which reported E =0.306 MPa and k =0.0011 mm$4/$N$·$s, while assuming ν = 0. Fiber parameters were found by choosing $αf=0$ and $βf=2.5$ and using a simple uniaxial tension model in FEBio with two orthogonal fibers to show [55]'s tensile modulus of 14.1 MPa corresponded to $ξ=17.5$ MPa. The choice of $βf=2.5$ allows a continuous compressive-to-tensile modulus transition as seen experimentally [49]; however, per the discussion in Sec. 4.3, superficial zone fibers had $βf=2$ prescribed instead.

The results for both CC and GC sliding show very high fluid load support sustained over the full 3600 s testing duration (Fig. 8(b)), with the fluid load support only decreasing from 97.8% to 93.4% in CC and from 97.9% to 95.6% in GC during this time. These elevated values of the fluid load support contributed to calculated friction coefficients which remained low throughout the analysis duration (Fig. 8(c)), achieving final values of $μ=0.018$ for CC5 and $μ=0.032$ for GC. These values are an order of magnitude below $μeq=0.201$ and in close agreement with the experiment (Fig. 8(d)), which reported final values of $μ=0.022$ and $μ=0.038$ for CC and GC, respectively. As the experimental data shown here represent only a single sample, we may also compare the equilibrium friction to the average of experiments E1 and E2 in Ref. [10], who found final values of $μ=0.016$ for CC and $μ=0.027$ for GC. The match obtained between friction measured experimentally and predicted by the present algorithm by using reasonable properties from the literature clearly shows the influence of the contact modality on sustained IFP and hence friction, also indicating that the biphasic friction algorithm correctly captures these physics.

Contour plots of the fluid pressure after 155 s and 3595 s for both CC and GC models (Fig. 9) show the mechanisms by which such low friction is accomplished. Due to the presence of two deformable surfaces, the CC model has lower fluid pressure by 155 s compared to CC and exhibits a greater drop by 3595 s (Figs. 9(a) and 9(b)). The model is still able to sustain elevated fluid load support because the contact area is initially greater and further grows over time such that contact pressures are also correspondingly smaller. In contrast, the GC case exhibits higher fluid pressure and a smaller contact area compared to the CC model due to the rigid impermeable glass lens. As the contact area barely grows during 3600 s of sliding, the fluid pressure only decreases slightly over time (Figs. 9(c) and 9(d)).

Fig. 9
Fig. 9
Close modal

## 5 Discussion

This study has presented the formulation, implementation, and validation of a finite element algorithm for large deformation frictional contact between porous-permeable hydrated materials. The algorithm builds on our recently developed elastic frictional contact algorithm [27] and our earlier work on frictionless biphasic contact [19,32], resulting in a conceptually simple framework which is robust enough to solve complex large deformation problems which often arise in biomechanics. To the best of our knowledge, this study provides the first formulation and implementation of a frictional algorithm for hydrated soft tissues, valid under finite deformation and large sliding.

The numerical examples presented in Sec. 4 were selected to verify and validate the contact framework by replicating classical experiments in the cartilage friction literature, demonstrating the ability of the presented algorithm to faithfully reproduce complex frictional behavior unique to hydrated porous media. Verification was achieved by replicating theoretical results for infinitesimal strain SCA contact under both constant and cyclical loads (Figs. 1 and 2); these illustrations also served to verify the linearity of the relationship between $μloc$ and $−p/tn$ embodied in Eq. (2.5) (Fig. 1(b)). Validation was performed by replicating experimental results for three canonical experiments in the cartilage friction literature. It is well known that SCA contact leads to depressurization of cartilage and a large increase in the friction coefficient over time [1,4,10,42]. When complete depressurization occurs, the tissue achieves its equilibrium friction coefficient $μ=μeq$. Additionally, applying a cyclical load to cartilage during SCA experiments results in pumping the interstitial fluid, leading to subambient fluid pressures producing suction against the glass slide during a fraction of each loading cycle [9]. This suction may increase the friction coefficient to values significantly above the depressurized equilibrium value $μeq$. Validation of these unique behaviors was accomplished by fitting creep data for a single plug used in constant and cyclical loading experiments, and then predicting the resulting effective friction coefficients for those scenarios using the same material parameters (Sec. 4.5). Our friction algorithm was able to produce μ as high as 0.51 under cyclical loading in a model where $μeq=0.2565$ was prescribed (Fig. 5), a unique effect that can only occur in a hydrated material where load-sharing is split between fluid and solid constituents. Finally, the influence of contact modality on tissue IFP was explored under MCA contact, computationally replicating the familiar sustained elevated fluid load support and low friction seen experimentally (Fig. 8).

Though the developed algorithms have been extensively validated by comparison with classical experiments in the cartilage friction literature, it is important to remember that cartilage friction is complex and arises from a multitude of factors. Elevated IFP is required for low friction, and can only be achieved with anisotropic material models that include fiber reinforcement and a strong tension-compression nonlinearity [37,39,41]. The transient response of the measured friction coefficient is also reliant upon the tissue hydraulic permeability, which controls the rate of depressurization, and cartilage permeability is known to be strain-dependent [56]. If these complex effects are not captured accurately with an appropriate material model, the finite element algorithm may not be able to reproduce experimental results because a contact algorithm that relies on IFP is limited by the behavior of the bulk material, not just the properties of the bearing surfaces. Beyond physical considerations, there are finite element limitations as well. Due to the high water content of cartilage, the depressurization which accompanies elevated friction can induce very large deformations in cartilage and hence the finite element mesh. Since the stiffness of the tissue is partly due to IFP, the large mesh deformation occurs as the biphasic material depressurizes and becomes softer. These effects often combine with the elevated friction coefficient produced by depressurization to create distorted meshes or a failure to converge as equilibrium approaches, e.g., in SCA contact. Again, these issues are not strictly due to the contact algorithm, which only exists to apply the correct boundary conditions due to contact, yet they generally will not arise in frictionless contact since the mesh does not experience frictional tractions in that case. When attempting to perform or debug a frictional biphasic contact analysis, it is important for the user to be cognizant of both physical modeling considerations and finite element limitations, which may be affecting the results independent of the contact algorithm. In regular geometries of the type which do not strictly exist in reality (e.g., rectangular or cylindrical specimens with perfect right-angled corners), performing mesh refinement often leads to smaller elements near the edge of the contact producing negative Jacobians, as the depressurized elements do not have the stiffness to resist shear deformation and the node on the corner collapses inward. This is a finite element meshing difficulty which is revealed by using frictional algorithms applied to porous media. In this study, a mesh convergence analysis was performed, showing that all variables of interest (including IFP and friction coefficients) changed by less than 3% with further mesh refinement.

Mesh limitations also can affect the calculated friction coefficient μ during MCA contact, Eq. (4.2), and produce values inconsistent with $μeff$ calculated from the fluid load support as per Eqs. (2.3) or (2.5). This occurs because experimental cartilage friction is so exceedingly low during MCA, on the order of 0.02, that mesh discretization error can easily introduce errors of 50% or more in the finite element value of μ in this range of low values. To understand this, consider a frictionless model of a lens on strip with a coarse mesh in the travel direction. As the lens migrates, the contact might be exactly symmetric on either side, but it is more likely that one side will have more integration points in contact due to the discretization of the contact surfaces. In this scenario, the tangential contact tractions would not cancel out exactly (producing $F≠0$), even during frictionless contact ($μeq=μeff=0$), and thus would produce an apparent nonzero finite element friction coefficient μ according to Eq. (4.2). Since experimental MCA friction is so low, the discretization-induced finite element μ may not be necessarily negligible compared to experimental values, given a poor enough mesh. Fortunately, this issue may be easily identified by running frictionless sliding contact analyses with a given mesh ($μeq=0$); if the finite element μ is found to be non-negligible compared to experimental values, the mesh should be refined.

The guidelines provided in Sec. 4.3 for selecting constitutive models for the purposes of biphasic friction simulation may appear restrictive. It is important to remember, however, that users must define finite element models with a specific purpose in mind, and this purpose guides the necessary level of detail. Accordingly, the Sec. 4.3 guidelines are formulated to model cartilage friction in great detail, also producing accurate contact tractions and stress distributions. The suggestion of modeling at least two distinct tissue layers, with different fiber orientations and properties in each layer, is especially important in MCA contact to ensure fluid load support does not produce negative values of $μloc$. In contrast, if friction is deemed to have a negligible effect on a particular analysis, it may be sufficient to use a single layer model with isotropic continuous fiber reinforcement and $μeq=0$, as done in some of our earlier studies [48,57]. Similarly, though continuous fiber distributions represent physiological fibril distributions, they tend to require more computational time. If the effective behavior of cartilage is required, but not necessarily the most demanding physiological structural model, one may choose to use two or three orthogonal fiber bundles as in Sec. 4.6. These are determinations the user must make.

An additional novelty in the present work is the development of smoothing methods for Lagrange augmentation (Sec. 3.2.2) and fluid load support calculations (Sec. 4.2); these two methods represent options users must decide whether to use, along with the choice of penalty versus augmented Lagrangian enforcement of contact constraints. The latter must generally be decided on a case-by-case basis, as both methods (penalty and augmented Lagrange) work well depending upon the application. In our experience, it is most efficient to begin by selecting the penalty method, and then evaluating the contact gap and fluid pressures in the output file. The purpose of contact enforcement is to minimize the contact gap to a degree de-emed appropriate, while ensuring continuity of fluid pressures across a contact interface. If increasing the penalty scale factor produces poor numerical convergence, augmented Lagrangian enforcement is adopted. If Lagrange augmentations are chosen, we recommend applying augmentation smoothing (“smooth_aug”) as well, particularly in cases where the spatial distribution of contact variables is of interest, as this method provides cleaner results. For similar reasons, we suggest the fluid load support smoothing option (“smooth_fls”) always be used as well.

While the verification and validation examples presented in this study relied on some of our own prior experiments, we emphasize that multiple other researchers have made significant advancements to the field of cartilage lubrication [1,58]. Here, we benefited from the availability of raw experimental data, as well as our thorough familiarity with the experimental conditions employed in those studies. The formulation developed in this work, though only verified and validated for articular cartilage, remains valid for frictional boundary contact between porous materials whose frictional response is coupled to IFP in the manner described by Eq. (2.3). Another potential application is in hydrogel mechanics, where recent experimental work has shown a lubrication mechanism strikingly similar to that in cartilage [59]. However, as frictional relations are constitutive in nature, users who wish to apply this work to materials other than cartilage must perform their own validation to evaluate the appropriateness of this framework for their application.

The formulation presented here can account for boundary lubrication as well. The effects of different surface chemistries may be simulated by absorbing them into the user-defined equilibrium friction term $μeq$. For example, boundary lubrication by hyaluronic acid or lubricin, which is known to influence the value of $μeq$ [10,60], may be accounted for by looking up literature values of $μeq$ for the desired surface chemistry effect. Mixed lubrication may also be accounted for by specifying a value of $φ$ which takes into account separation between asperities by pressurized fluid pools on the contacting surfaces. The current formulation considers $φ$ to be a user-defined constant. It would be straightforward to modify this treatment to let $φ=φs(1)φs(2)$, where $φs(i)$ are the solid volume fractions of primary and secondary surfaces i at the integration points in the current configuration, or to introduce a constitutive relation accounting for pore deformation during contact; these modifications would allow $φ$ to evolve with contact, completing the description of mixed lubrication. Users who desire a constitutive relation for $φ$ have the ability to implement such a model as a plug-in [61]. In the present work, a constant $φ$ was sufficient for validation purposes.

Given that this work employs a boundary contact algorithm, there are certain lubrication phenomena that it cannot simulate in a straightforward manner. For example, Burris, Moore, and coworkers recently described a hydrodynamic rehydration mechanism in cartilage, wherein pressurized fluid may be entrained into the cartilage layers at the leading edge of contact for sufficiently high sliding speeds (∼60 mm/s), possibly bolstering IFP [11,12,62]. Fortunately, this type of mechanism may be simulated using our recently developed biphasic fluid–structure interaction module in FEBio [30], which complements the current boundary contact algorithm, offering a way to resolve fluid-film lubrication problems.

This study has addressed an important need in the biomechanics community by formulating and implementing a finite element algorithm for large deformation frictional contact of hydrated materials. Importantly, this algorithm has been implemented into the free, open-source finite element software FEBio, which is available to the general public and specifically tailored to the biomechanics and biophysics research communities [26]. In contrast to the straightforward application of Coulomb's law as used for modeling finite element frictional contact between elastic solids [22,27], friction laws for hydrated materials must account for the porous nature of the contacting materials and resulting interstitial fluid pressurization. This algorithm makes it possible to model frictional responses that were previously unobtainable by computational methods.

## Footnotes

2

Not to be confused with $Te$, introduced in Eq. (2.1) as the effective stress arising from solid matrix deformation.

3

We have previously shown that fibers are in tension when they fall inside or outside the elliptical cone of normal strain [45,46], which occurs when the three principal normal strains at a point exhibit a mix of tensile and compressive values. This mix of compressive and tensile principal normal strains is typically significant at non-conforming biphasic contact interfaces [14,47].

4

These values have been converted from those in Ref. [9], who reported aggregate tensile and compressive moduli of $H+A=13.2$ MPa and $H−A=0.64$ MPa, respectively, with the off-diagonal modulus set at $λ2=0.48$ MPa.

5

Note that inserting $μeq=0.201$ and $φ=0.0144$ into Eq. (2.5) and using the equilibrium CC fluid load support of 93.4% produces $μeff=0.016$, slightly below this value. This difference is addressed in Sec. 5.

## Acknowledgment

This study was supported with funds from the National Institute of General Medical Sciences (Grant No. GM083925) and by the National Science Foundation Graduate Research Fellowship Program (No. DGE-11-44155).

### Appendix A: Contact Kinematics

The contact kinematics developed in this Appendix are summarized from the detailed treatment of Zimmerman and Ateshian [27], which should be consulted for further detail.

##### A.1 Slip Kinematics.
The kinematics of slip are developed by mapping points between the surfaces $γ(i)$ as they move relative to one another. The relationship between spatial points $x(i)(η(i)α,t)$ on each surface is given by
$x(2)=x(1)+gn(1)$
(A1)
where $n(1)$ is the unit outward normal to $γ(1)$ given by
$n(1)=g1(1)×g2(1)|g1(1)×g2(1)|$
(A2)
and the gap function g is defined as
$g=(x(2)−x(1))·n(1)$
(A3)

The relationship of Eq. (A2) is graphically shown in Fig. 10. Here, we note that $x(1)(η(1)α,t)$ is the spatial position of a material point $X(1)$ on the primary surface, and $x(2)(η(2)α,t)$ is the corresponding spatial intersection point on the secondary surface, through which different material points $X(2)$ identified by parametric coordinates $η(2)α$ may convect. At any given instant, $η(2)α$ are termed the parametric coordinates of intersection.

Fig. 10
Fig. 10
Close modal
##### A.2 Stick Kinematics.
Implicit in the concept of stick is the assumption that the contact projection was previously resolved; thus, points are not mapped between surfaces during stick. Rather, the current contact point on $γ(2)$ is given by the parametric coordinates of intersection $η(2)α$ from the previous time point, now denoted as $η(2)pα$ with the subscripted p referring to the previous time. The spatial position of the material point $X(2)$ identified by $η(2)pα$ is then given by
$xs(2)=x(2)(η(2)pα,t)=x(1)(η(1)α,t)+gs$
(A4)
where the vector gap $gs$ in stick is defined to be
$gs=xs(2)−x(1)$
(A5)

Here, $gs$ 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 $gs=0$.

##### A.3 Velocities.

Coulomb's law of kinetic friction requires that the friction force be aligned with the relative slip velocity between the two surfaces. Despite the name, Coulomb's law is a constitutive relation, and hence must obey the Principle of Material-Frame Indifference; this requires a frame-invariant relative velocity. As points in stick do not experience relative motion, the development of velocities below is only concerned with opposing contact points in slip.

As parametric coordinates $η(1)α$ of integration points on the primary surface $γ(1)$ represent material coordinates, the velocity $v(1)$ of these points is evaluated from the material time derivative in the material frame
$v(1)(η(1)α,t)=∂x(1)(η(1)α,t)∂t$
(A6)
In contrast, different material points $X(2)$ may convect through the intersection point $x(2)$, and so the velocity $v(2)$ at the intersection point on $γ(2)$ is evaluated from the material time derivative in the spatial frame
$v(2)(η(2)α,t)=∂x(2)(η(2)α,t)∂t+η˙(2)αgα(2)$
(A7)
where $∂x(2)/∂t$ represents the velocity of the intersection point on $γ(2)$, while $η˙(2)α$ are the contravariant components of the convective velocity of material passing through the intersection point $x(2)$. The quantity $η˙(2)αgα(2)$ represents the relative slip velocity between the material on $γ(2)$ and that on $γ(1)$. Importantly, by definition $∂x(2)/∂t$ is evaluated while holding $η(2)α$ constant. From these relations, a more practical formulation of the slip velocity can be achieved [27]. Taking the material time derivative of Eq. (A1) and recalling the contact persistency condition $g˙=0$ [24] produces $v(2)=v(1)+gn˙(1)$, then substituting Eqs. (A6) and (A7) into this expression yields the desired frame-invariant measure of relative velocity between $γ(1)$ and $γ(2)$ [63],
$vr≡η˙(2)αgα(2)=gn˙(1)+∂x(1)(η(1)α,t)∂t−∂x(2)(η(2)α,t)∂t$
(A8)
where $n˙(1)$ is evaluated from Eq. (A2) as
$n˙(1)=PN·g˙1(1)×g2(1)+g1(1)×g˙2(1)Jη(1)$
(A9)
Here, $g˙α(1)$ is the material time derivative of $gα(1)$ in the material frame, evaluated from Eq. (3.3) as
$g˙α(1) (η(1)β,t)=∂gα(1)(η(1)β,t)∂t$
(A10)
and $PN$ is the tangential plane projection tensor,
$PN=I−n(1)⊗n(1)$
(A11)
A unit vector in the slip direction can then be found by projecting the relative velocity $vr$ onto the tangent plane of $γ(1)$, yielding
$s(1)=PN(1)·vr|PN(1)·vr|$
(A12)

### Appendix B: Constitutive Models

The strain energy density of the compressible Holmes–Mow solid is given by [50]
$ΨHM=14λ+2μβhm(eQ−1)$
(B1)
where
$Q=βhmλ+2μ[(2μ−λ)(I1−3)+λ(I2−3)−(λ+2μ)lnJ2]$
(B2)
and I1 and I2 are the first and second invariants of the right Cauchy-Green deformation tensor C. Here, $J=detF$ is the Jacobian of the deformation and λ and μ are Lamé coefficients related to Young's modulus E and Poisson's ratio ν through
$λ=E(1+ν)(1−2ν)μ=E2(1+ν)$
The Newtonian viscous stress in a Kelvin–Voigt viscoelastic solid is given by
$TNVv=(κv−23μv)(trD)I+2μvD$
(B3)

where κv and μv are the bulk and shear viscosities and D is the rate-of-deformation tensor.

A single fiber bundle with exponential power law behavior has the strain energy density
$Ψf=ξαfβf(exp [αf(In−1)βf]−1)$
(B4)
and the associated stress
$Tf=H(In−1)2InJ∂Ψf∂Inn⊗n$
(B5)

where $In=λn2$ is the square of the fiber stretch, $H(·)$ is the Heaviside step function, and $n=F·N/λn$ is a unit vector in the fiber direction in the current configuration, whereas N is the unit vector in the fiber direction in the reference frame.

The strain energy density for the neo-Hookean solid is given by Bonet and Wood [33]
$ΨNH=μ(12(I1−3)−lnJ+ν1−2ν(lnJ)2)$
(B6)
Strain-dependent Holmes-Mow permeability makes use of the constitutive relation
$k(J)=k0(J−φ01−φ0)αpe12m(J2−1)$

where $φ0$ is the solidity and the corresponding permeability tensor is $k=k(J)I$.

## References

1.
Ateshian
,
G. A.
,
2009
, “
The Role of Interstitial Fluid Pressurization in Articular Cartilage Lubrication
,”
J. Biomech.
,
42
(
9
), pp.
1163
1176
.10.1016/j.jbiomech.2009.04.040
2.
McCutchen
,
C.
,
1959
, “
Mechanism of Animal Joints: Sponge-Hydrostatic and Weeping Bearings
,”
Nature
,
184
(
4695
), pp.
1284
1285
.10.1038/1841284a0
3.
McCutchen
,
C. W.
,
1962
, “
The Frictional Properties of Animal Joints
,”
Wear
,
5
(
1
), pp.
1
17
.10.1016/0043-1648(62)90176-X
4.
Forster
,
H.
, and
Fisher
,
J.
,
1996
, “
The Influence of Loading Time and Lubricant on the Friction of Articular Cartilage
,”
Proc. Inst. Mech. Eng., Part H
,
210
(
2
), pp.
109
119
.10.1243/PIME_PROC_1996_210_399_02
5.
Ateshian
,
G.
,
1997
, “
A Theoretical Formulation for Boundary Friction in Articular Cartilage
,”
ASME J. Biomech. Eng.
,
119
(
1
), pp.
81
86
.10.1115/1.2796069
6.
Ateshian
,
G.
,
Wang
,
H.
, and
Lai
,
W.
,
1998
, “
The Role of Interstitial Fluid Pressurization and Surface Porosities on the Boundary Friction of Articular Cartilage
,”
ASME J. Tribol.
,
120
(
2
), pp.
241
248
.10.1115/1.2834416
7.
Soltz
,
M. A.
, and
Ateshian
,
G. A.
,
1998
, “
Experimental Verification and Theoretical Prediction of Cartilage Interstitial Fluid Pressurization at an Impermeable Contact Interface in Confined Compression
,”
J. Biomech.
,
31
(
10
), pp.
927
934
.10.1016/S0021-9290(98)00105-5
8.
Krishnan
,
R.
,
Kopacz
,
M.
, and
Ateshian
,
G. A.
,
2004
, “
Experimental Verification of the Role of Interstitial Fluid Pressurization in Cartilage Lubrication
,”
J. Orthop. Res.
,
22
(
3
), pp.
565
570
.10.1016/j.orthres.2003.07.002
9.
Krishnan
,
R.
,
Mariner
,
E. N.
, and
Ateshian
,
G. A.
,
2005
, “
Effect of Dynamic Loading on the Frictional Response of Bovine Articular Cartilage
,”
J. Biomech.
,
38
(
8
), pp.
1665
1673
.10.1016/j.jbiomech.2004.07.025
10.
Caligaris
,
M.
, and
Ateshian
,
G. A.
,
2008
, “
Effects of Sustained Interstitial Fluid Pressurization Under Migrating Contact Area, and Boundary Lubrication by Synovial Fluid, on Cartilage Friction
,”
Osteoarthritis Cartilage
,
16
(
10
), pp.
1220
1227
.10.1016/j.joca.2008.02.020
11.
Burris
,
D. L.
, and
Moore
,
A. C.
,
2017
, “
Cartilage and Joint Lubrication: New Insights Into the Role of Hydrodynamics
,”
Biotribology
,
12
, pp.
8
14
.10.1016/j.biotri.2017.09.001
12.
Moore
,
A. C.
, and
Burris
,
D. L.
,
2017
, “
Tribological Rehydration of Cartilage and Its Potential Role in Preserving Joint Health
,”
Osteoarthritis Cartilage
,
25
(
1
), pp.
99
107
.10.1016/j.joca.2016.09.018
13.
Northwood
,
E.
,
Fisher
,
J.
, and
Kowalski
,
R.
,
2007
, “
Investigation of the Friction and Surface Degradation of Innovative Chondroplasty Materials Against Articular Cartilage
,”
Proc. Inst. Mech. Eng., Part H
,
221
(
3
), pp.
263
79
.10.1243/09544119JEIM178
14.
Ateshian
,
G. A.
, and
Wang
,
H.
,
1995
, “
A Theoretical Solution for the Frictionless Rolling Contact of Cylindrical Biphasic Articular Cartilage Layers
,”
J. Biomech.
,
28
(
11
), pp.
1341
1355
.10.1016/0021-9290(95)00008-6
15.
Donzelli
,
P. S.
, and
Spilker
,
R. L.
,
1998
, “
A Contact Finite Element Formulation for Biological Soft Hydrated Tissues
,”
Comput. Methods Appl. Mech. Eng.
,
153
(
1–2
), pp.
63
79
.10.1016/S0045-7825(97)00065-0
16.
Chen
,
X.
,
Chen
,
Y.
, and
Hisada
,
T.
,
2005
, “
Development of a Finite Element Procedure of Contact Analysis for Articular Cartilage With Large Deformation Based on the Biphasic Theory
,”
JSME Int. J., Ser. C
,
48
(
4
), pp.
537
546
.10.1299/jsmec.48.537
17.
Yang
,
T.
, and
Spilker
,
R. L.
,
2007
, “
A Lagrange Multiplier Mixed Finite Element Formulation for Three-Dimensional Contact of Biphasic Tissues
,”
ASME Biomech. Eng.
,
129
(
3
), pp.
457
471
.10.1115/1.2737056
18.
Guo
,
H.
, and
Spilker
,
R. L.
,
2011
, “
Biphasic Finite Element Modeling of Hydrated Soft Tissue Contact Using an Augmented Lagrangian Method
,”
ASME J. Biomech. Eng.
,
133
(
11
), p.
111001
.10.1115/1.4005378
19.
Ateshian
,
G. A.
,
Maas
,
S.
, and
Weiss
,
J. A.
,
2010
, “
Finite Element Algorithm for Frictionless Contact of Porous Permeable Media Under Finite Deformation and Sliding
,”
ASME J. Biomech.
,
132
(
6
), p.
061006
.10.1115/1.4001034
20.
Federico
,
S.
,
La Rosa
,
G.
,
Herzog
,
W.
, and
Wu
,
J. Z.
,
2004
, “
Effect of Fluid Boundary Conditions on Joint Contact Mechanics and Applications to the Modeling of Osteoarthritic Joints
,”
ASME J. Biomech. Eng.
,
126
(
2
), pp.
220
225
.10.1115/1.1691445
21.
Pawaskar
,
S. S.
,
Fisher
,
J.
, and
Jin
,
Z.
,
2010
, “
Robust and General Method for Determining Surface Fluid Flow Boundary Conditions in Articular Cartilage Contact Mechanics Modeling
,”
ASME J. Biomech. Eng.
,
132
(
3
), p.
031001
.10.1115/1.4000869
22.
Laursen
,
T.
, and
Simo
,
J.
,
1993
, “
A Continuum-Based Finite Element Formulation for the Implicit Solution of Multibody, Large Deformation-Frictional Contact Problems
,”
Int. J. Numer. Methods Eng.
,
36
(
20
), pp.
3451
3485
.10.1002/nme.1620362005
23.
Wriggers
,
P.
,
2002
,
Computational Contact Mechanics
,
Wiley & Sons
,
Hoboken, NJ
.
24.
Wriggers
,
P.
, and
Laursen
,
T. A.
,
2007
,
Computational Contact Mechanics, CISM Courses and Lectures
,
Springer
,
Wien
.
25.
Freutel
,
M.
,
Schmidt
,
H.
,
Dürselen
,
L.
,
Ignatius
,
A.
, and
Galbusera
,
F.
,
2014
, “
Finite Element Modeling of Soft Tissues: Material Models, Tissue Interaction and Challenges
,”
Clin. Biomech.
,
29
(
4
), pp.
363
372
.10.1016/j.clinbiomech.2014.01.006
26.
Maas
,
S. A.
,
Ellis
,
B. J.
,
Ateshian
,
G. A.
, and
Weiss
,
J. A.
,
2012
, “
FEBio: Finite Elements for Biomechanics
,”
ASME J. Biomech. Eng.
,
134
(
1
), p.
011005
.10.1115/1.4005694
27.
Zimmerman
,
B. K.
, and
Ateshian
,
G. A.
,
2018
, “
A Surface-to-Surface Finite Element Algorithm for Large Deformation Frictional Contact in FEBio
,”
ASME J. Biomech. Eng.
,
140
(
8
), p.
081013
.10.1115/1.4040497
28.
Mow
,
V. C.
,
Kuei
,
S.
,
Lai
,
W. M.
, and
Armstrong
,
C. G.
,
1980
, “
Biphasic Creep and Stress Relaxation of Articular Cartilage in Compression: Theory and Experiments
,”
ASME J. Biomech. Eng.
,
102
(
1
), pp.
73
84
.10.1115/1.3138202
29.
Soltz
,
M. A.
,
Basalo
,
I. M.
, and
Ateshian
,
G. A.
,
2003
, “
Hydrostatic Pressurization and Depletion of Trapped Lubricant Pool During Creep Contact of a Rippled Indenter Against a Biphasic Articular Cartilage Layer
,”
ASME J. Biomech. Eng.
,
125
(
5
), pp.
585
93
.10.1115/1.1610020
30.
Shim
,
J. J.
,
Maas
,
S. A.
,
Weiss
,
J. A.
, and
Ateshian
,
G. A.
,
2021
, “
Finite Element Implementation of Biphasic-Fluid Structure Interactions in FEBio
,”
ASME J. Biomech. Eng.
,
143
(
9
), p.
09100
.10.1115/1.4050646
31.
Ateshian
,
G. A.
,
Maas
,
S.
, and
Weiss
,
J. A.
,
2012
, “
Solute Transport Across a Contact Interface in Deformable Porous Media
,”
J. Biomech.
,
45
(
6
), pp.
1023
1027
.10.1016/j.jbiomech.2012.01.003
32.
Maas
,
S.
, Herron, M.,
Weiss
,
J.
, and
Ateshian
,
G. A.
,
2021
, “
FEBio 3.4 Theory Manual
,” FEBio, accessed Sept. 2, 2021, https://help.febio.org/FEBioTheory/FEBio_tm_3-4.html
33.
Bonet
,
J.
, and
Wood
,
R. D.
,
2008
,
Nonlinear Continuum Mechanics for Finite Element Analysis
, 2nd ed.,
Cambridge University Press
,
Cambridge, UK
.
34.
Simo
,
J.
, and
Laursen
,
T.
,
1992
, “
An Augmented Lagrangian Treatment of Contact Problems Involving Friction
,”
Comput. Struct.
,
42
(
1
), pp.
97
116
.10.1016/0045-7949(92)90540-G
35.
Bertsekas
,
D. P.
,
1982
,
Constrained Optimization and Lagrange Multiplier Methods
,
Academic Press
,
New York
.
36.
Soltz
,
M. A.
, and
Ateshian
,
G. A.
,
2000
, “
Interstitial Fluid Pressurization During Confined Compression Cyclical Loading of Articular Cartilage
,”
Ann. Biomed. Eng.
,
28
(
2
), pp.
150
159
.10.1114/1.239
37.
Park
,
S.
,
Krishnan
,
R.
,
Nicoll
,
S. B.
, and
Ateshian
,
G. A.
,
2003
, “
Cartilage Interstitial Fluid Load Support in Unconfined Compression
,”
J. Biomech.
,
36
(
12
), pp.
1785
1796
.10.1016/S0021-9290(03)00231-8
38.
Burris
,
D.
, and
Sawyer
,
W.
,
2009
, “
Addressing Practical Challenges of Low Friction Coefficient Measurements
,”
Tribol. Lett.
,
35
(
1
), pp.
17
23
.10.1007/s11249-009-9438-2
39.
Huang
,
C.-Y.
,
Soltz
,
M. A.
,
Kopacz
,
M.
,
Mow
,
V. C.
, and
Ateshian
,
G. A.
,
2003
, “
Experimental Verification of the Roles of Intrinsic Matrix Viscoelasticity and Tension-Compression Nonlinearity in the Biphasic Response of Cartilage
,”
ASME J. Biomech. Eng.
,
125
(
1
), pp.
84
93
.10.1115/1.1531656
40.
Cohen
,
B.
,
Lai
,
W.
, and
Mow
,
V.
,
1998
, “
A Transversely Isotropic Biphasic Model for Unconfined Compression of Growth Plate and Chondroepiphysis
,”
ASME J. Biomech. Eng.
,
120
(
4
), pp.
491
496
.10.1115/1.2798019
41.
Soltz
,
M. A.
, and
Ateshian
,
G. A.
,
2000
, “
A Conewise Linear Elasticity Mixture Model for the Analysis of Tension-Compression Nonlinearity in Articular Cartilage
,”
ASME J. Biomech. Eng.
,
122
(
6
), pp.
576
586
.10.1115/1.1324669
42.
Ateshian
,
G. A.
,
Soltz
,
M. A.
,
Mauck
,
R. L.
,
Basalo
,
I. M.
,
Hung
,
C. T.
, and
Lai
,
W. M.
,
2003
, “
The Role of Osmotic Pressure and Tension-Compression Nonlinearity in the Frictional Response of Articular Cartilage
,”
Transp. Porous Media
,
50
(
1/2
), pp.
5
33
.10.1023/A:1020618514874
43.
Armstrong
,
C.
,
Lai
,
W.
, and
Mow
,
V.
,
1984
, “
An Analysis of the Unconfined Compression of Articular Cartilage
,”
ASME J. Biomech. Eng.
,
106
(
2
), pp.
165
173
.10.1115/1.3138475
44.
Ateshian
,
G. A.
,
Rajan
,
V.
,
Chahine
,
N. O.
,
Canal
,
C. E.
, and
Hung
,
C. T.
,
2009
, “
Modeling the Matrix of Articular Cartilage Using a Continuous Fiber Angular Distribution Predicts Many Observed Phenomena
,”
ASME J. Biomech. Eng.
,
131
(
6
), p.
061003
.10.1115/1.3118773
45.
Ateshian
,
G. A.
,
2007
, “
Anisotropy of Fibrous Tissues in Relation to the Distribution of Tensed and Buckled Fibers
,”
ASME J. Biomech. Eng.
,
129
(
2
), pp.
240
249
.10.1115/1.2486179
46.
Hou
,
C.
, and
Ateshian
,
G. A.
,
2016
, “
A Gauss-Kronrod-Trapezoidal Integration Scheme for Modeling Biological Tissues With Continuous Fiber Distributions
,”
Comput. Methods Biomech. Biomed Eng.
,
19
(
8
), pp.
883
93
.10.1080/10255842.2015.1075518
47.
Ateshian
,
G.
,
Lai
,
W.
,
Zhu
,
W.
, and
Mow
,
V.
,
1994
, “
An Asymptotic Solution for the Contact of Two Biphasic Cartilage Layers
,”
J. Biomech.
,
27
(
11
), pp.
1347
1360
.10.1016/0021-9290(94)90044-2
48.
Jones
,
B.
,
Hung
,
C. T.
, and
Ateshian
,
G.
,
2015
, “
Biphasic Analysis of Cartilage Stresses in the Patellofemoral Joint
,”
J. Knee Surg.
,
29
(
02
), pp.
092
098
.10.1055/s-0035-1568989
49.
Chahine
,
N. O.
,
Wang
,
C. C.
,
Hung
,
C. T.
, and
Ateshian
,
G. A.
,
2004
, “
Anisotropic Strain-Dependent Material Properties of Bovine Articular Cartilage in the Transitional Range From Tension to Compression
,”
J. Biomech.
,
37
(
8
), pp.
1251
1261
.10.1016/j.jbiomech.2003.12.008
50.
Holmes
,
M.
, and
Mow
,
V. C.
,
1990
, “
The Nonlinear Characteristics of Soft Gels and Hydrated Connective Tissues in Ultrafiltration
,”
J. Biomech.
,
23
(
11
), pp.
1145
1156
.10.1016/0021-9290(90)90007-P
51.
Ateshian
,
G.
,
Warden
,
W.
,
Kim
,
J.
,
Grelsamer
,
R.
, and
Mow
,
V.
,
1997
, “
Finite Deformation Biphasic Material Properties of Bovine Articular Cartilage From Confined Compression Experiments
,”
J. Biomech.
,
30
(
11–12
), pp.
1157
1164
.10.1016/S0021-9290(97)85606-0
52.
Zhu
,
W.
,
Mow
,
V. C.
,
Koob
,
T. J.
, and
Eyre
,
D. R.
,
1993
, “
Viscoelastic Shear Properties of Articular Cartilage and the Effects of Glycosidase Treatments
,”
J. Orthop. Res.
,
11
(
6
), pp.
771
781
.10.1002/jor.1100110602
53.
Buckley
,
M. R.
,
Gleghorn
,
J. P.
,
Bonassar
,
L. J.
, and
Cohen
,
I.
,
2008
, “
Mapping the Depth Dependence of Shear Properties in Articular Cartilage
,”
J. Biomech.
,
41
(
11
), pp.
2430
2437
.10.1016/j.jbiomech.2008.05.021
54.
Guerin
,
H. A.
, and
Elliott
,
D. M.
,
2005
, “
The Role of Fiber-Matrix Interactions in a Nonlinear Fiber-Reinforced Strain Energy Model of Tendon
,”
ASME J. Biomech. Eng.
,
127
(
2
), pp.
345
350
.10.1115/1.1865212
55.
Oungoulian
,
S. R.
,
Hehir
,
K. E.
,
Zhu
,
K.
,
Willis
,
C. E.
,
Marinescu
,
A. G.
,
Merali
,
N.
,
Ahmad
,
C. S.
,
Hung
,
C. T.
, and
Ateshian
,
G. A.
,
2014
, “
Effect of Glutaraldehyde Fixation on the Frictional Response of Immature Bovine Articular Cartilage Explants
,”
J. Biomech.
,
47
(
3
), pp.
694
701
.10.1016/j.jbiomech.2013.11.043
56.
Lai
,
W.
,
Mow
,
V. C.
, and
Roth
,
V.
,
1981
, “
Effects of Nonlinear Strain-Dependent Permeability and Rate of Compression on the Stress Behavior of Articular Cartilage
,”
ASME J. Biomech. Eng.
,
103
(
2
), pp.
61
66
.10.1115/1.3138261
57.
Todd
,
J. N.
,
Allan
,
A. N.
,
Maak
,
T. G.
, and
Weiss
,
J. A.
,
2021
, “
Characterization and Finite Element Validation of Transchondral Strain in the Human Hip During Static and Dynamic Loading
,”
J. Biomech.
,
114
, p.
110143
.10.1016/j.jbiomech.2020.110143
58.
Ateshian
,
G. A.
,
Mow
,
V.
, and
Huiskes
,
R.
,
2005
, “
Friction, Lubrication, and Wear of Articular Cartilage and Diarthrodial Joints
,”
Basic Orthopaedic Biomechanics & Mechano-Biology
, Lippincott Williams & Wilkins, Philadelphia, PA, Vol.
3
, pp.
447
494
.
59.
Delavoipière
,
J.
,
Tran
,
Y.
,
Verneuil
,
E.
,
Heurtefeu
,
B.
,
Hui
,
C. Y.
, and
Chateauminois
,
A.
,
2018
, “
Friction of Poroelastic Contacts With Thin Hydrogel Films
,”
Langmuir
,
34
(
33
), pp.
9617
9626
.10.1021/acs.langmuir.8b01466
60.
Schmidt
,
T. A.
,
Gastelum
,
N. S.
,
Nguyen
,
Q. T.
,
Schumacher
,
B. L.
, and
Sah
,
R. L.
,
2007
, “
Boundary Lubrication of Articular Cartilage: Role of Synovial Fluid Constituents
,”
Arthritis Rheumatism
,
56
(
3
), pp.
882
891
.10.1002/art.22446
61.
Maas
,
S. A.
,
LaBelle
,
S. A.
,
Ateshian
,
G. A.
, and
Weiss
,
J. A.
,
2018
, “
A Plugin Framework for Extending the Simulation Capabilities of FEBio
,”
Biophys. J.
,
115
(
9
), pp.
1630
1637
.10.1016/j.bpj.2018.09.016
62.
Burris
,
D.
,
Ramsey
,
L.
,
Graham
,
B.
,
Price
,
C.
, and
Moore
,
A.
,
2019
, “
How Sliding and Hydrodynamics Contribute to Articular Cartilage Fluid and Lubrication Recovery
,”
Tribol. Lett.
,
67
(
2
), pp.
1
10
.10.1007/s11249-019-1158-7
63.
Curnier
,
A.
,
He
,
Q.-C.
, and
Klarbring
,
A.
,
1995
, “
Continuum Mechanics Modelling of Large Deformation Contact With Friction
,”
Contact Mechanics
,
Springer
, Boston, MA, pp.
145
158
.