## 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 $\mu eq$, which is a boundary friction property of the load-bearing surfaces of the solid matrix alone, and the experimentally recorded effective friction coefficient $\mu 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 $\mu eff$ unless otherwise specified. Due to interactions between fluid and solid phases, time- and load-dependent effects produce values of $\mu 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 $\mu 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 [7–12]. 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 [2–4]. 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 [15–18], 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 [22–24]. 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

**T**in a biphasic material is given by

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=\u222bAn\xb7T\xb7n\u2009dA$, 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=\u2212\u222bAp\u2009dA$, 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\u2212\phi )Wp$, where $1\u2212\phi $ 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\u2212(1\u2212\phi )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 *W ^{ss}* becomes equal to the total applied load

*W*.

*F*at the interface results primarily from solid-on-solid contact, and thus is directly proportional to

*W*,

^{ss}where $\mu eq$ is the equilibrium friction coefficient achieved when $Wp=0$. For two idealized perfectly smooth surfaces in contact, $\phi $ 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., $\phi =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 $\phi $. 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, $\phi $ must be obtained experimentally, along with $\mu eq$ [6]; both parameters will therefore need to be specified by the user in computational models (Sec. 4.1).

*F*/

*W*becomes

and it is apparent that in this formulation, $\mu eff$ is a linear function of the interstitial fluid load support, reaching a minimum value of $\mu eff=\phi \mu 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,\u2009\phi =0)$ the friction coefficient reduces to $\mu eff=0$. Furthermore, in the idealized case of two perfectly smooth surfaces in contact ($\phi =1$), the model of Eq. (2.3) correctly reduces to Coulomb's law and produces $\mu eff=\mu eq$.

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 $\u2212p/tn\u22641$. However, on thermodynamic grounds, it is impossible to have a negative friction coefficient $\mu loc$ in a continuum framework. Therefore, the theoretical upper bound on $\u2212p/tn$ is $(1\u2212\phi )\u22121$. It must be emphasized that $\mu 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 $\phi $ 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.

*t*, in terms of contravariant surface parametric coordinates $\eta (i)\alpha $. Casting Eq. (3.1) into convenient matrix notation and switching the domain of integration to the parametric frame yields the invariant biphasic contact integral

*f*in this biphasic analysis,

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.

*w*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 $\Psi $, where on the primary surface

_{n}*t*and

_{n}*p*. The value of the slip criterion determines the stick-slip status via

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.

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

*g*of the gap, given by Eq. (A3),

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 $\mu 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.

*w*can be obtained by penalizing the fluid pressure gap

_{n}*π*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

and $\epsilon p$ is the pressure penalty parameter, which has units of hydraulic permeability per unit length (e.g., m^{3}/N$\xb7$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)(\eta (1)\alpha )$ is evaluated at a point with parametric coordinates $\eta (1)\alpha $ on the primary surface, whereas the pressure on the secondary surface, $p(2)=p(2)(\eta (2)\alpha ,t)$, is evaluated at the parametric coordinates of the intersection of a ray issued from that primary surface point and normal to $\gamma (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 *w _{n}* from Eq. (3.12) using the contact kinematics determined by Eq. (3.8). In practice, this is accomplished by calculating

*w*once the stick-slip status has been resolved for each iteration.

_{n}#### 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.

*g,*

*λ*is the normal Lagrange multiplier. The total traction vector in slip is then directly prescribed as

_{n}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.

*λ*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}Unlike $\lambda s$ and *λ _{n}* [27],

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

_{p}*λ*according to

_{n}*λ*and derive the passive multiplier $\lambda s$

_{n}In the above formulas, $Tel$ is the total Cauchy stress in the solid element attached to the contact face^{2} 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.

## 4 Verifications and Validations

### 4.1 Finite Element Implementation.

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

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

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

_{i}*k*is a measure of the average hydraulic permeability, in the direction normal to the surface. The parameters

_{i}*α*and

*α*are user-defined scale factors. When the autopenalty feature is not used, the penalties reduce to these user-defined parameters, $\epsilon =\alpha $ and $\epsilon p=\alpha 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.

_{p}The primary user-supplied parameters for biphasic friction analyses are the equilibrium friction coefficient $\mu eq$ and the solid–solid contact area fraction $\phi $. For users unfamiliar with these concepts, we briefly remark upon best practices for determining the correct values to apply. Experimentally, $\mu eq$ may only be measured in the context of SCA friction, where typically one counterface is impermeable; once measured, however, values of $\mu 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\u2264\mu eq\u22640.30$ for articular cartilage, though we emphasize that there are no restrictions upon this parameter beyond $\mu eq>0$. Similarly, $\phi $ 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, $\phi \u223d0.1$ for contact between cartilage and an impermeable surface, and $\phi \u223d0.01$ for cartilage-on-cartilage contact. As it represents an area fraction, the only restriction is $0\u2264\phi \u22641$.

*W*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,

^{p}Location-dependent friction coefficients $\mu loc$ (and fluid load support $\u2212p/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 $\mu 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 $\mu loc$ in Eq. (2.5) in the frictional biphasic contact algorithm. Using the model for $\mu 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 $\u2212p/tn$ at the contact surface integration points. However, this evaluation at the surface is susceptible to numerical noise for the reasons given above for *t _{n}*, especially when using the augmented Lagrange method. This spatial noise becomes problematic when evaluating $\mu loc$, especially when a noisy $\u2212p/tn$ exceeds the theoretical upper bound $(1\u2212\phi )\u22121$; in those cases, $\mu 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 $\mu 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 $\mu loc$, which was specified by Eq. (2.4) as $\u2212p/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 $\u2212pel/(n(1)\xb7Tel\xb7n(1))$, where $Tel$ is the element-averaged total Cauchy stress and

*p*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 –

_{el}*p*and

*t*throughout this paper, with the understanding that $\mu loc$ has been calculated by turning on this flag.

_{n}### 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 [40–42] 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)\xb7T\xb7n(1)$ of the contact traction may be expressed as $tn=\u2212p+tne$ when using Eq. (2.1), where $tne=n(1)\xb7Te\xb7n(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 $\u2212p/tn>(1\u2212\phi )\u22121$.^{3} This nonsensical construct would then produce an interstitial fluid load support $\u2212p/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\u2212\phi )\u22121$ due to the presence of nontangential fibers at greater depths from the contact surface, which propagate these elevated values of $\u2212p/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 $\u2212p/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 $\u2212p/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 $\mu 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$\xb7$s, with CLE parameters^{4}$\lambda +1=12.86$ MPa, $\lambda \u22121=0.3$ MPa, $\lambda 2=0.48$ MPa, and $\mu =0.17$ MPa. Here, $\lambda +1,\u2009\lambda \u22121$, and *λ*_{2} represent tensile, compressive, and off-diagonal first Lamé coefficients, and *μ* is the second Lamé coefficient. The friction parameters were given by $\mu eq=0.153$ and $\phi =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 $\alpha =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\u2212$intercept given by the regression is 0.153, exactly in agreement with the prescribed value of $\mu eq$, which is the value of $\mu eff$ theoretically achieved when fluid load support becomes zero. Similarly, $\mu =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 $\mu eq$ and $\phi $ 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.

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 $\mu 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=\u22120.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.

### 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 $\mu eq=0.2565$. The solid–solid contact fraction was prescribed to be $\phi =0.05$ such that *μ* could better match experimental data when $Wp/W\u21921$.

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, $\nu =0.3$, and $\beta 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}*β*, and

_{f}*ξ*, where

*α*is a power-law coefficient,

_{f}*β*is an exponential coefficient, and

_{f}*ξ*is related to the tensile modulus. For both layers, $\alpha f=0$ and

*ξ*= 40 MPa were selected. The parameter

*β*was chosen to be 2 in the superficial zone and 2.5 in the deep zone, where $\beta f=2$ ensures a piecewise linear curve when the modulus transitions from compressive to tensile and $\beta 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 $\kappa v=0\u2009MPa\xb7s$ and $\mu v=1\u2009MPa\xb7s$, where

_{f}*κ*is the bulk viscosity and

_{v}*μ*the shear viscosity of the solid. Finally, strain-dependent Holmes-Mow permeability [50] was prescribed, with $\alpha p=0,\u2009k0=0.0045$ mm$4/$N$\xb7$s, and

_{v}*m*=

*8. Here,*

*k*

_{0}represents the zero-strain permeability,

*m*is an exponential stiffening coefficient, and

*α*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].

_{p}#### 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.

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 $\mu min=0.028$ to $\mu =0.236$. As likely occurred in the experiment, the fluid pressure has not entirely subsided by 2500 s and thus $\mu eq$ was not attained. Though a good prediction of *μ* was achieved for this sliding configuration, Krishnan et al. [9] report $\mu 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 $\phi $: the theoretical minimum friction coefficient $\mu min$ can be calculated from Eq. (2.5) when $\u2212p/tn=1$ as $\mu min=\phi \mu eq=0.0128$, using our prescribed values for $\phi $ and the experimental $\mu eq$; thus our choice of $\phi $ is unable to reproduce the experimental $\mu min=0.0089$. A smaller value of $\phi $ would resolve this discrepancy, especially when accounting for the fact that $\phi $ 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 $\mu eq=0.2565$ as obtained from the constant load experiment in Sec. 4.5.1; in that experiment, *μ* rose asymptotically toward $\mu 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 $\mu =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 $\mu eq$ validates our implementation against an effect unique to hydrated tissues whose frictional response is coupled to the IFP.

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 $\mu 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 $\mu 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.

### 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 $\mu 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 $\phi =0.12$ for glass–cartilage (GC) contact and $\phi =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 $\alpha 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.

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$\xb7$s, while assuming *ν* = 0. Fiber parameters were found by choosing $\alpha f=0$ and $\beta 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 $\xi =17.5$ MPa. The choice of $\beta 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 $\beta 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 $\mu =0.018$ for CC^{5} and $\mu =0.032$ for GC. These values are an order of magnitude below $\mu eq=0.201$ and in close agreement with the experiment (Fig. 8(d)), which reported final values of $\mu =0.022$ and $\mu =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 $\mu =0.016$ for CC and $\mu =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)).

## 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 $\mu loc$ and $\u2212p/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 $\mu =\mu 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 $\mu 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 $\mu 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 $\mu 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\u22600$), even during frictionless contact ($\mu eq=\mu 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 ($\mu 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 $\mu 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 $\mu 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 $\mu eq$. For example, boundary lubrication by hyaluronic acid or lubricin, which is known to influence the value of $\mu eq$ [10,60], may be accounted for by looking up literature values of $\mu eq$ for the desired surface chemistry effect. Mixed lubrication may also be accounted for by specifying a value of $\phi $ which takes into account separation between asperities by pressurized fluid pools on the contacting surfaces. The current formulation considers $\phi $ to be a user-defined constant. It would be straightforward to modify this treatment to let $\phi =\phi s(1)\phi s(2)$, where $\phi 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 $\phi $ to evolve with contact, completing the description of mixed lubrication. Users who desire a constitutive relation for $\phi $ have the ability to implement such a model as a plug-in [61]. In the present work, a constant $\phi $ 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

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

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

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\u2212A=0.64$ MPa, respectively, with the off-diagonal modulus set at $\lambda 2=0.48$ MPa.

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

*g*is defined as

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

##### A.2 Stick Kinematics.

*p*referring to the previous time. The spatial position of the material point $X(2)$ identified by $\eta (2)p\alpha $ is then given by

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.

### Appendix B: Constitutive Models

*I*

_{1}and

*I*

_{2}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

where *κ _{v}* and

*μ*are the bulk and shear viscosities and

_{v}**D**is the rate-of-deformation tensor.

where $In=\lambda n2$ is the square of the fiber stretch, $H(\xb7)$ is the Heaviside step function, and $n=F\xb7N/\lambda 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.

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