7.2 Biphasic Contact

Contact Integral

See Section Biphasic Material for a review of biphasic materials, reference ¹ for additional details on biphasic frictionless contact, and reference ² for biphasic frictional contact. The presentation here follows that of ². The contact interface is defined between surfaces \(\gamma^{\left(1\right)}\) and \(\gamma^{\left(2\right)}\). Due to continuity requirements on the traction and fluxes, the external virtual work resulting from contact tractions \(\mathbf{t}^{\left(i\right)}\) and solvent fluxes \(w_{n}^{\left(i\right)}\) (\(i=1,2)\), may be combined into the contact integral

\[ \begin{equation} \delta G_{c}=\int_{\gamma^{(1)}}\left(\left(\delta\mathbf{v}^{(1)}-\delta\mathbf{v}^{(2)}\right)\cdot\mathbf{t}^{(1)}+\left(\delta p^{(1)}-\delta p^{(2)}\right)w_{n}\right)da^{(1)}\label{eq:contact-int} \end{equation} \]

where \(\mathbf{t}^{(1)}\) is the contact traction on the primary surface, which is equal and opposite to that on the secondary surface, \(\mathbf{t}^{(1)}=-\mathbf{t}^{(2)}\); \(w_{n}\equiv w_{n}^{(1)}\) is the outward normal component of the fluid flux \(\mathbf{w}^{(1)}\) from the primary surface, \(\delta\mathbf{v}^{(i)}\) are virtual velocities, \(\delta p^{(i)}\) are virtual fluid pressures, and \(da^{(1)}\) is an elemental area on the primary surface \(\gamma^{(1)}\). The contact integral is then written over the invariant parametric space of \(\gamma^{(1)}\), which is denoted by \(\Gamma_{\eta}^{(1)}\) ³, facilitating its linearization as required for use with an iterative solution method such as Newton's method. The elemental area is given by

\[ \begin{equation} \begin{aligned}da^{(1)} & =J_{\eta}^{(1)}d\eta_{(1)}^{1}d\eta_{(1)}^{2}, & J_{\eta}^{(1)} & =\left|\mathbf{g}_{1}^{(1)}\times\mathbf{g}_{2}^{(1)}\right|\end{aligned} \label{eq:da-jeta} \end{equation} \]

using the covariant basis vectors

\[ \begin{equation} \mathbf{g}_{\alpha}^{(i)}=\frac{\partial\mathbf{x}^{(i)}}{\partial\eta_{(i)}^{\alpha}}\label{eq:basis-vectors} \end{equation} \]

where \(\mathbf{x}^{(i)}\left(\eta_{(i)}^{\alpha},t\right)\) is the spatial representation of surface \(\gamma^{(i)}\) as it deforms in time \(t\), in terms of contravariant surface parametric coordinates \(\eta_{(i)}^{\alpha}\). Casting eq.\eqref{eq:contact-int} into convenient matrix notation and switching the domain of integration to the parametric frame yields the invariant biphasic contact integral

\[ \begin{equation} \delta G_{c}=\int_{\Gamma_{\eta}^{(1)}}\left(\begin{bmatrix}\delta\mathbf{v}^{(1)} & \delta p^{(1)}\end{bmatrix}\cdot\begin{bmatrix}\mathbf{t}^{(1)}\\ w_{n} \end{bmatrix}+\begin{bmatrix}\delta\mathbf{v}^{(2)} & \delta p^{(2)}\end{bmatrix}\cdot\begin{bmatrix}-\mathbf{t}^{(1)}\\ -w_{n} \end{bmatrix}\right)J_{\eta}^{(1)}d\eta_{(1)}^{1}d\eta_{(1)}^{2}\label{eq:contact-int-invariant-1} \end{equation} \]

Since \(\Gamma_{\eta}^{(1)}\) represents an invariant material frame, the linearization of \(\delta G_{c}\) can be accomplished by applying the directional derivative operator directly to the integrand,

\[ \begin{equation} D\delta G_{c}=\int_{\Gamma_{\eta}^{(1)}}D\left(\begin{bmatrix}\delta\mathbf{v}^{(1)} & \delta p^{(1)}\end{bmatrix}\cdot\begin{bmatrix}\mathbf{t}^{(1)}\\ w_{n} \end{bmatrix}J_{\eta}^{(1)}+\begin{bmatrix}\delta\mathbf{v}^{(2)} & \delta p^{(2)}\end{bmatrix}\cdot\begin{bmatrix}-\mathbf{t}^{(1)}\\ -w_{n} \end{bmatrix}J_{\eta}^{(1)}\right)d\eta_{(1)}^{1}d\eta_{(1)}^{2}\label{eq:D-contact-int} \end{equation} \]

where it is understood that for any function \(f\) in this biphasic analysis,

\[ \begin{equation} Df\equiv\sum_{i=1}^{2}Df\left[\Delta\mathbf{u}^{(i)}\right]+Df\left[\Delta p^{(i)}\right] \end{equation} \]

Biphasic Friction Formulation

An examination of eq.\eqref{eq:contact-int-invariant-1} reveals that a contact framework must provide expressions for both the normal fluid flux \(w_{n}\) and the contact traction \(\mathbf{t}^{(1)}\). For the contact traction, this work implements a Coulomb-like framework for frictional contact ⁴. In this framework the relationship between sticking and slipping, and thus the contact traction \(\mathbf{t}^{(i)}\) on the opposing surfaces, is described by a slip criterion \(\Psi\), where on the primary surface

\[ \begin{equation} \Psi=\left|\mathbf{t}_{T}^{(1)}\right|-\mu_{\text{eff}}\left|t_{n}\right| \end{equation} \]

where

\[ \begin{equation} \mu_{\text{eff}}=\mu_{\text{eq}}\left[1+\left(1-\varphi\right)\frac{p}{t_{n}}\right]\label{eq:mueff-biphasic} \end{equation} \]

Here, \(t_{n}=\mathbf{t}^{(1)}\cdot\mathbf{n}^{(1)}\) is the normal component of the contact traction (negative in compression), \(\mathbf{t}_{T}^{(1)}=\mathbf{t}^{(1)}-t_{n}\mathbf{n}^{(1)}\) is the tangential component of the contact traction, and the effective friction coefficient \(\mu_{\text{eff}}\) is given by eq.\eqref{eq:mueff-biphasic} with no distinction made between static and kinetic coefficients of friction. The user-defined parameter \(\mu_{\text{eq}}\) represents the friction coefficient when the fluid pressure has subsided (\(p=0\)), whereas the user-defined parameter \(\varphi\) represents the fraction of the apparent contact area where solid-on-solid contact takes place ⁵. The ratio \(-p/t_{n}\) represents the local fluid load support at each point on the contact surface. Since the friction coefficient \(\mu_{\text{eff}}\) cannot be negative, the theoretical upper bound on the local fluid load support is \(\left.-p/t_{n}\right|_{\max}=\left(1-\varphi\right)^{-1}\). This upper bound is enforced whenever numerical errors produce a greater value of the local fluid load support.

The value of the slip criterion determines the stick-slip status via

\[ \begin{equation} \Psi\begin{cases} <0 & \text{sticking}\\ =0 & \text{slipping} \end{cases}\label{eq:slip-criterion-1} \end{equation} \]

Following our prior study ⁴, 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 (Section Tied Biphasic Contact), 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.

Contact Kinematics

Slip Kinematics

The kinematics of slip are developed by mapping points between the surfaces \(\gamma^{(i)}\) as they move relative to one another. The relationship between spatial points \(\mathbf{x}^{(i)}\left(\eta_{(i)}^{\alpha},t\right)\) on each surface is given by

\[ \begin{equation} \mathbf{x}^{(2)}=\mathbf{x}^{(1)}+g\mathbf{n}^{(1)}\label{eq:slip-x2} \end{equation} \]

where \(\mathbf{n}^{(1)}\) is the unit outward normal to \(\gamma^{(1)}\) given by

\[ \begin{equation} \mathbf{n}^{(1)}=\frac{\mathbf{g}_{1}^{(1)}\times\mathbf{g}_{2}^{(1)}}{\left|\mathbf{g}_{1}^{(1)}\times\mathbf{g}_{2}^{(1)}\right|}\label{eq:n-1} \end{equation} \]

and the gap function \(g\) is defined as

\[ \begin{equation} g=\left(\mathbf{x}^{(2)}-\mathbf{x}^{(1)}\right)\cdot\mathbf{n}^{(1)}\label{eq:slip-gap} \end{equation} \]

Here we note that \(\mathbf{x}^{(1)}\left(\eta_{(1)}^{\alpha},t\right)\) is the spatial position of a material point \(X^{(1)}\) on the primary surface, and \(\mathbf{x}^{(2)}\left(\eta_{(2)}^{\alpha},t\right)\) 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.

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 \(\gamma^{(2)}\) is given by the parametric coordinates of intersection \(\eta_{(2)}^{\alpha}\) from the previous time point, now denoted as \(\eta_{(2)p}^{\alpha}\) with the subscripted \(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

\[ \begin{equation} \mathbf{x}_{s}^{(2)}=\mathbf{x}^{(2)}\left(\eta_{(2)p}^{\alpha},t\right)=\mathbf{x}^{(1)}\left(\eta_{(1)}^{\alpha},t\right)+\mathbf{g}_{s}\label{eq:stick-x2} \end{equation} \]

where the vector gap \(\mathbf{g}_{s}\) in stick is defined to be

\[ \begin{equation} \mathbf{g}_{s}=\mathbf{x}_{s}^{(2)}-\mathbf{x}^{(1)}\label{eq:stick-gap} \end{equation} \]

Here \(\mathbf{g}_{s}\) is the vectorial distance, at the current time \(t\), between material points which were in contact at the previous time step; for perfect stick we must have \(\mathbf{g}_{s}=\mathbf{0}\).

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 \(\eta_{(1)}^{\alpha}\) of integration points on the primary surface \(\gamma^{(1)}\) represent material points, the velocity \(\mathbf{v}^{(1)}\) of these points is evaluated from the material time derivative in the material frame,

\[ \begin{equation} \mathbf{v}^{(1)}\left(\eta_{(1)}^{\alpha},t\right)=\frac{\partial\mathbf{x}^{(1)}\left(\eta_{(1)}^{\alpha},t\right)}{\partial t}\label{eq:v1-1} \end{equation} \]

In contrast, different material points \(X^{(2)}\) may convect through the intersection point \(\mathbf{x}^{(2)}\), and so the velocity \(\mathbf{v}^{(2)}\) at the intersection point on \(\gamma^{(2)}\) is evaluated from the material time derivative in the spatial frame,

\[ \begin{equation} \mathbf{v}^{(2)}\left(\eta_{(2)}^{\alpha},t\right)=\frac{\partial\mathbf{x}^{(2)}\left(\eta_{(2)}^{\alpha},t\right)}{\partial t}+\dot{\eta}_{(2)}^{\alpha}\mathbf{g}_{\alpha}^{(2)}\label{eq:v2} \end{equation} \]

where \(\partial\mathbf{x}^{(2)}/\partial t\) represents the velocity of the intersection point on \(\gamma^{(2)}\), while \(\dot{\eta}_{(2)}^{\alpha}\) are the contravariant components of the convective velocity of material passing through the intersection point \(\mathbf{x}^{(2)}\). The quantity \(\dot{\eta}_{(2)}^{\alpha}\mathbf{g}_{\alpha}^{(2)}\) represents the relative slip velocity between the material on \(\gamma^{(2)}\) and that on \(\gamma^{(1)}\). Importantly, by definition \(\partial\mathbf{x}^{(2)}/\partial t\) is evaluated while holding \(\eta_{(2)}^{\alpha}\) constant. From these relations, a more practical formulation of the slip velocity can be achieved ⁴. Taking the material time derivative of eq.\eqref{eq:slip-x2} and recalling the contact persistency condition \(\dot{g}=0\) ⁶ produces \(\mathbf{v}^{(2)}=\mathbf{v}^{(1)}+g\dot{\mathbf{n}}^{(1)}.\) Substituting Eqs.\eqref{eq:v1-1}-\eqref{eq:v2} into this expression yields the desired frame-invariant measure of relative velocity between \(\gamma^{(1)}\) and \(\gamma^{(2)}\) ⁷,

\[ \begin{equation} \mathbf{v}^{r}\equiv\dot{\eta}_{(2)}^{\alpha}\mathbf{g}_{\alpha}^{(2)}=g\dot{\mathbf{n}}^{(1)}+\frac{\partial\mathbf{x}^{(1)}\left(\eta_{(1)}^{\alpha},t\right)}{\partial t}-\frac{\partial\mathbf{x}^{(2)}\left(\eta_{(2)}^{\alpha},t\right)}{\partial t}, \end{equation} \]

where \(\dot{\mathbf{n}}^{(1)}\) is evaluated from eq.\eqref{eq:n-1} as

\[ \begin{equation} \dot{\mathbf{n}}^{(1)}=\mathbf{P}_{N}\cdot\frac{\dot{\mathbf{g}}_{1}^{(1)}\times\mathbf{g}_{2}^{(1)}+\mathbf{g}_{1}^{(1)}\times\dot{\mathbf{g}}_{2}^{(1)}}{J_{\eta}^{(1)}}\,. \end{equation} \]

Here, \(\dot{\mathbf{g}}_{\alpha}^{(1)}\) is the material time derivative of \(\mathbf{g}_{\alpha}^{(1)}\) in the material frame, evaluated from eq.\eqref{eq:basis-vectors} as

\[ \begin{equation} \dot{\mathbf{g}}_{\alpha}^{(1)}\left(\eta_{\left(1\right)}^{\beta},t\right)=\frac{\partial\mathbf{g}_{\alpha}^{(1)}\left(\eta_{\left(1\right)}^{\beta},t\right)}{\partial t}\,, \end{equation} \]

and \(\mathbf{P}_{N}\) is the tangential plane projection tensor,

\[ \begin{equation} \mathbf{P}_{N}=\mathbf{I}-\mathbf{n}^{(1)}\otimes\mathbf{n}^{(1)}\,. \end{equation} \]

A unit vector in the slip direction can then be found by projecting the relative velocity \(\mathbf{v}^{r}\) onto the tangent plane of \(\gamma^{(1)}\), yielding

\[ \begin{equation} \mathbf{s}^{(1)}=\frac{\mathbf{P}_{N}^{(1)}\cdot\mathbf{v}^{r}}{\left|\mathbf{P}_{N}^{(1)}\cdot\mathbf{v}^{r}\right|}\,.\label{eq:s-def} \end{equation} \]

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,

\[ \begin{equation} \mathbf{t}^{(1)}=\varepsilon\mathbf{g}_{s}\,,\label{eq:stick-traction} \end{equation} \]

where \(\varepsilon\) is the contact penalty parameter and we have utilized eq.\eqref{eq:stick-gap} to define the gap vector \(\mathbf{g}_{s}\) 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.\eqref{eq:slip-gap},

\[ \begin{equation} t_{n}=\varepsilon g\,.\label{eq:slip-tn} \end{equation} \]

The total traction vector in slip is then directly prescribed as

\[ \begin{equation} \mathbf{t}^{(1)}=t_{n}\left(\mathbf{n}^{(1)}+\mu_{\text{eff}}\mathbf{s}^{(1)}\right)\,,\label{eq:slip-traction} \end{equation} \]

where \(\mathbf{s}^{(1)}\) is the unit vector in the slip direction, given by eq.\eqref{eq:s-def}. We note that eq.\eqref{eq:slip-traction} has the same form as the classical Coulomb friction considered in Section Sliding-Elastic and ⁴, with the standard Coulomb friction coefficient replaced by \(\mu_{\text{eff}}\), as evaluated from \eqref{eq:mueff-biphasic}. A trial state and return map, adapted from Section Stick-Slip Algorithm ⁴ and presented in Section Stick-Slip Algorithm, is employed to differentiate between stick and slip.

The jump conditions for a biphasic mixture require continuity of the fluid pressure across an interface and therefore an expression for the normal fluid flux \(w_{n}\) can be obtained by penalizing the fluid pressure gap \(\pi\) between contacting points ⁸,

\[ \begin{equation} \begin{aligned}w_{n} & =\varepsilon_{p}\pi\,, & t_{n}<0\\ p^{(i)} & =0\,, & t_{n}=0 \end{aligned} \label{eq:wn} \end{equation} \]

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

\[ \begin{equation} \pi=p^{(1)}-p^{(2)}\label{eq:pi-biphasic-penalty} \end{equation} \]

and \(\varepsilon_{p}\) is the pressure penalty parameter which has units of hydraulic permeability per unit length (e.g. m\(^{3}\)/N\(\cdot\)s, similar to the hydraulic permeability of a membrane). Equation \eqref{eq:wn} is valid for both stick and slip. Note that \(p^{\left(1\right)}=p^{\left(1\right)}\left(\eta_{\left(1\right)}^{\alpha}\right)\) is evaluated at a point with parametric coordinates \(\eta_{\left(1\right)}^{\alpha}\) on the primary surface, whereas the pressure on the secondary surface, \(p^{(2)}=p^{(2)}\left(\eta_{(2)}^{\alpha},t\right)\), is evaluated at the parametric coordinates of the intersection of a ray issued from that primary surface point and normal to \(\gamma^{\left(1\right)}\). As detailed previously ⁴ and summarized in SectionContact Kinematics, the parametric coordinates of intersection are dependent on the stick-slip status, so care must be taken to evaluate \(w_{n}\) from eq.\eqref{eq:wn} using the contact kinematics determined by eq.\eqref{eq:slip-criterion-1}. In practice this is accomplished by calculating \(w_{n}\) once the stick-slip status has been resolved for each iteration.

Augmented Lagrangian Scheme

The augmented Lagrangian scheme presented herein was developed in Section Augmented Lagrangian Scheme ⁴ as a modification of the approach proposed by Simo and Laursen ⁹. Briefly, this is a first-order augmentation scheme that utilizes Uzawa's algorithm ¹⁰, where the multipliers are updated outside of the Newton step, producing a double loop algorithm ⁶ and preserving quadratic convergence of Newton's method near solution points.

During stick, the traction is calculated by augmenting the vector gap \(\mathbf{g}_{s}\),

\[ \begin{equation} \mathbf{t}^{(1)}=\boldsymbol{\lambda}_{s}+\varepsilon\mathbf{g}_{s}\,,\label{eq:AL-stick-traction} \end{equation} \]

where \(\boldsymbol{\lambda}_{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\),

\[ \begin{equation} t_{n}=\lambda_{n}+\varepsilon g\,,\label{eq:AL-tn} \end{equation} \]

where \(\lambda_{n}\) is the normal Lagrange multiplier. The total traction vector in slip is then directly prescribed as

\[ \begin{equation} \mathbf{t}^{(1)}=t_{n}\left(\mathbf{n}^{(1)}+\mu_{\text{eff}}\mathbf{s}^{(1)}\right)\label{eq:augmented-t} \end{equation} \]

where eq.\eqref{eq:AL-tn} has been used. The update formulas for the Lagrange multipliers can be found in Section Augmented Lagrangian Scheme ⁴. Here it suffices to note that by augmenting only the normal gap and employing eq.\eqref{eq:augmented-t}, 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 Section Stick-Slip Algorithm 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

\[ \begin{equation} \begin{aligned}w_{n} & =\lambda_{p}+\varepsilon_{p}\pi\,, & t_{n}<0\\ p^{(i)} & =0\,, & t_{n}=0 \end{aligned} \label{eq:AL-wn} \end{equation} \]

where \(\lambda_{p}\) is the fluid pressure Lagrange multiplier, and the update formula for \(\lambda_{p}\) has been detailed in our prior work on frictionless biphasic contact ⁸ and is given by

\[ \begin{equation} \lambda_{p}\leftarrow\lambda_{p}+\varepsilon_{p}\pi \end{equation} \]

Unlike \(\boldsymbol{\lambda}_{s}\) and \(\lambda_{n}\) ⁴, \(\lambda_{p}\) is not dependent on the stick-slip status. However, as noted before, eq.\eqref{eq:AL-wn} must be evaluated after the stick-slip status has been resolved.

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 ⁴. We initially assume stick and calculate a trial traction \(\tilde{\mathbf{t}}^{(1)}\), utilizing either eq.\eqref{eq:stick-traction} or eq.\eqref{eq:AL-stick-traction}. The trial normal and tangential components \(\tilde{t}_{n}\) and \(\tilde{\mathbf{t}}_{T}^{(1)}\) are evaluated from \(\tilde{\mathbf{t}}^{(1)}\) and inserted into the slip criterion \(\Psi\),

\[ \begin{equation} \Psi=\left|\tilde{\mathbf{t}}_{T}^{(1)}\right|-\mu_{\text{eff}}\left|\tilde{t}_{n}\right|. \end{equation} \]

Based on the slip criterion and trial traction vector, the traction vector is calculated from the return mapping as

\[ \begin{equation} \mathbf{t}^{(1)}=\begin{cases} \tilde{\mathbf{t}}^{(1)} & \Psi<0\,,\quad\text{sticking}\\ t_{n}(\mathbf{n}^{(1)}+\mu_{\text{eff}}\mathbf{s}^{(1)}) & \Psi=0\,,\quad\text{slipping} \end{cases} \end{equation} \]

where \(t_{n}\) is given by either eq.\eqref{eq:slip-tn} or eq.\eqref{eq:AL-tn}.

Linearization and Discretization Outline

The following sections provide the linearization and discretization of stick and slip for biphasic frictional contact, culminating in the final forms of the residual vectors and stiffness matrices which have been implemented into FEBio. Section Definitions and Notation begins with an important discussion of kinematic definitions and notation which is assumed to hold without repetition. Sections Biphasic Stick and Biphasic Slip linearize and discretize biphasic stick and slip, respectively. Importantly, Section Biphasic Stick sets up and uses more formal notation and definitions, before we dispense of this cumbersome notation for the remainder of the treatment.

Definitions and Notation

Evaluating the linearizations of eq.\eqref{eq:contact-int-invariant-1} requires directional derivatives of kinematic quantities, some of which depend on the stick-slip status. To simplify the presentation, the continuum linearization is presented only for a few select quantities. The remainder of the linearization is deferred until after discretization, as many expressions are much easier to manipulate in discretized form. To keep the equations more manageable, the discretization of slip is split into frictionless and frictional terms, and the final form of the stiffness matrices follows this split. The full model is easily obtained by summing the frictionless and frictional contributions. Here we emphasize that, as a consequence of the double-loop Uzawa algorithm discussed in Section Augmented Lagrangian Scheme (reprised from ⁴), all Lagrange multipliers are updated outside of each Newton step, thus \(D\lambda=0\), where \(\lambda\) is any Lagrange multiplier, i.e. \(\lambda=\lambda_{n},\boldsymbol{\lambda}_{s},\lambda_{p}\), etc. As a consequence, Lagrange multipliers do not appear in any of the linearized or discretized equations. In what follows, Greek indices are associated with covariant basis vectors on the contacting surfaces, and thus vary from 1 to 2. Repeated Greek indices indicate implicit summation over their range.

Many kinematic quantities remain unchanged from Section Sliding-Elastic on elastic frictional contact ⁴. We thus accept without repetition the definitions of \(\No\), \(\Nbo\), \(\mc\), \(\mb\), \(\Nt\), \(\Mc\), \(\Mb\), \(\mathbf{c}^{(1)}\), \(\mathbf{m}^{(1)}\), \(\mathbf{Q}^{(1)}\), \(\bar{\mathbf{A}}_{c}^{(1)}\), \(\Ac\), \(\mathbf{P}_{N}\), \(\mathbf{P}_{S},\) \(\mathbf{R}^{(1)}\), \(\mathbf{B}^{(1)}\), \(\mathbf{L}^{(1)}\), \(\mathbf{J}_{c}^{(1)}\), \(\So\), and \(\tilde{\mathbf{S}}^{(1)}\). These terms are all defined in Sections Frictionless Terms-Frictional Terms ⁴. Directional derivatives which will not be duplicated here include \(D\no\), \(D\dot{\mathbf{n}}^{(1)}\), \(D\so\), \(Dt_{n}=\varepsilon Dg\), \(Dg\), \(D\jeta\), \(D\eta_{(2)}^{\alpha}\), and \(D\mathbf{P}_{N}\) (see Sections Linearization, Frictionless Terms and Frictional Terms).

The final residual vectors and stiffness matrices presented below are written in integral form. In the FEBio implementation, a Gaussian quadrature scheme is adopted to perform numerical integration. A detailed treatment of Gaussian quadrature, and equations for numerically integrating the contact integrals and stiffnesses, may be found in Section Integration Scheme.

As a final implementation detail, we note that multiplying all biphasic entries in the residual vector and stiffness matrices by the time step leads to better convergence. In this context, biphasic refers to all terms which are not purely related to solid-solid contact. In the residual, this is all terms which are not \(\mathbf{f}_{m}^{(i)}\) (see below). Similarly, this includes all stiffness entries except the solid-solid \(\mathbf{K}_{mn}^{(i,j)}\) terms (see below).

Biphasic Stick

Linearization

In stick parametric coordinates are fixed on both surfaces, hence directional derivatives of \(\mathbf{x}^{(i)}\), \(p^{(i)}\), \(\delta\mathbf{v}^{(i)}\), and \(\delta p^{(i)}\) are given by

\[ \begin{equation} \begin{aligned}D\mathbf{x}^{(1)} & =\Delta\mathbf{u}^{(1)} & D\mathbf{x}^{(2)} & =\Delta\mathbf{u}^{(2)}\\ Dp^{(1)} & =\Delta p^{(1)} & Dp^{(2)} & =\Delta p^{(2)}\\ D\delta\mathbf{v}^{(1)} & =\mathbf{0} & D\delta\mathbf{v}^{(2)} & =\mathbf{0}\\ D\delta p^{(1)} & =0 & D\delta p^{(2)} & =0 \end{aligned} \label{eq:biphasic-stick-linearization} \end{equation} \]

From the definitions of eq.\eqref{eq:stick-gap} and eq.\eqref{eq:stick-traction}, and utilizing eq.\eqref{eq:biphasic-stick-linearization}, it follows that

\[ \begin{equation} D\mathbf{t}^{(1)}=\varepsilon\left(\Delta\mathbf{u}^{(2)}-\Delta\mathbf{u}^{(1)}\right)\,.\label{eq:biph-stick-Dt} \end{equation} \]

Similarly, application of eq.\eqref{eq:biphasic-stick-linearization} to Eqs.\eqref{eq:wn}-\eqref{eq:pi-biphasic-penalty} yields

\[ \begin{equation} Dw_{n}=\varepsilon_{p}\left(\Delta p^{(1)}-\Delta p^{(2)}\right)\label{eq:biph-stick-Dwn} \end{equation} \]

The biphasic contact integral of eq.\eqref{eq:contact-int}, written over an invariant domain, can then be linearized directly to find

\[ \begin{equation} \begin{aligned}D\delta G_{c}= & \int_{\Gamma_{\eta}^{(1)}}\left(\begin{bmatrix}\delta\mathbf{v}^{(1)} & \delta\mathbf{v}^{(2)}\end{bmatrix}\cdot\left(\begin{bmatrix}1\\ -1 \end{bmatrix}D\mathbf{t}^{(1)}J_{\eta}^{(1)}+\begin{bmatrix}\mathbf{t}^{(1)}\\ -\mathbf{t}^{(1)} \end{bmatrix}DJ_{\eta}^{(1)}\right)\right.\\ & \left.+\begin{bmatrix}\delta p^{(1)} & \delta p^{(2)}\end{bmatrix}\cdot\left(\begin{bmatrix}1\\ -1 \end{bmatrix}Dw_{n}J_{\eta}^{(1)}+\begin{bmatrix}w_{n}\\ -w_{n} \end{bmatrix}DJ_{\eta}^{(1)}\right)\right)d\eta_{(1)}^{1}d\eta_{(1)}^{2} \end{aligned} \label{eq:biph-cont-int-linearized} \end{equation} \]

Note that for convenience, terms in eq.\eqref{eq:biph-cont-int-linearized} have been grouped differently than in eq.\eqref{eq:D-contact-int}.

Discretization

Let the continuous variables on the primary and secondary surfaces be interpolated over element faces according to ⁸⁴

\[ \begin{equation} \begin{aligned}\delta\mathbf{v}^{(1)} & =\sum_{a=1}^{m^{(1)}}N_{a}^{(1)}\delta\mathbf{v}_{a}^{(1)} & \delta\mathbf{v}^{(2)} & =\sum_{b=1}^{m_{k}^{(2)}}N_{b}^{(2)}\delta\mathbf{v}_{b}^{(2)}\\ \Delta\mathbf{u}^{(1)} & =\sum_{c=1}^{m^{(1)}}N_{c}^{(1)}\Delta\mathbf{u}_{c}^{(1)} & \Delta\mathbf{u}^{(2)} & =\sum_{d=1}^{m_{k}^{(2)}}N_{d}^{(2)}\Delta\mathbf{u}_{d}^{(2)}\\ \delta p^{(1)} & =\sum_{a=1}^{m^{(1)}}N_{a}^{(1)}\delta p_{a}^{(1)} & \delta p^{(2)} & =\sum_{b=1}^{m_{k}^{(2)}}N_{b}^{(2)}\delta p_{b}^{(2)}\\ \Delta p^{(1)} & =\sum_{c=1}^{m^{(1)}}N_{c}^{(1)}\Delta p_{c}^{(1)} & \Delta p^{(2)} & =\sum_{d=1}^{m_{k}^{(2)}}N_{d}^{(2)}\Delta p_{d}^{(2)} \end{aligned} \label{eq:biphasic-stick-var-discretization} \end{equation} \]

where \(N_{a}^{(i)}\) represent interpolation functions on the element faces of \(\gamma^{(i)}\), \(m^{(1)}\) is the number of nodes and interpolation functions on each primary element face, \(m_{k}^{(2)}\) is the number of nodes and interpolation functions on the secondary element face which is intersected by the ray issued from the \(k\text{th}\) integration point on the primary element face, and \(\delta\mathbf{v}_{a}^{(i)}\), \(\Delta\mathbf{u}_{a}^{(i)}\),\(\delta p_{a}^{(i)}\), and \(\Delta p_{a}^{(i)}\) represent respective nodal values of \(\delta\mathbf{v}^{(i)}\), \(\Delta\mathbf{u}^{(i)}\), \(\delta p^{(i)}\), and \(\Delta p^{(i)}\). From this point forward the summation signs will be written simply as \(\sum_{a}\), where it is assumed they have the same meaning described above.

Inserting the discretization of eq.\eqref{eq:biphasic-stick-var-discretization} into eq.\eqref{eq:contact-int-invariant-1} yields

\[ \begin{equation} \begin{aligned}\delta G_{c} & =\sum_{a}\begin{bmatrix}\delta\mathbf{v}_{a}^{(1)} & \delta p_{a}^{(1)}\end{bmatrix}\cdot\int_{\Gamma_{\eta}^{(1)}}\begin{bmatrix}\mathbf{f}_{a}^{(1)}\\ w_{a}^{(1)} \end{bmatrix}J_{\eta}^{(1)}d\eta_{(1)}^{1}d\eta_{(1)}^{2}\\ & +\sum_{b}\begin{bmatrix}\delta\mathbf{v}_{b}^{(2)} & \delta p_{b}^{(2)}\end{bmatrix}\cdot\int_{\Gamma_{\eta}^{(1)}}\begin{bmatrix}\mathbf{f}_{b}^{(2)}\\ w_{b}^{(2)} \end{bmatrix}J_{\eta}^{(1)}d\eta_{(1)}^{1}d\eta_{(1)}^{2} \end{aligned} \label{eq:biph-stick-int-disc} \end{equation} \]

where the residuals are

\[ \begin{equation} \begin{aligned}\mathbf{f}_{a}^{(1)}= & N_{a}^{(1)}\mathbf{t}^{(1)}\,, & \mathbf{f}_{b}^{(2)}= & -N_{b}^{(2)}\mathbf{t}^{(1)}\\ w_{a}^{(1)}= & N_{a}^{(1)}w_{n}\,\,, & w_{b}^{(2)}= & -N_{b}^{(2)}w_{n} \end{aligned} \end{equation} \]

and \(\mathbf{t}^{(1)}\) is obtained from eq.\eqref{eq:stick-traction} in the penalty case and from eq.\eqref{eq:AL-stick-traction} if augmented Lagrangian regularization is employed. Similarly, \(w_{n}\) is calculated from either eq.\eqref{eq:wn} or eq.\eqref{eq:AL-wn} for penalty and augmented Lagrange methods, respectively.

We may now discretize individual terms and place them into matrix notation, anticipating their substitution into eq.\eqref{eq:biph-cont-int-linearized}. By placing eq.\eqref{eq:biph-stick-int-disc} into Eqs.\eqref{eq:biph-stick-Dt}-\eqref{eq:biph-stick-Dwn} and inserting the resulting linearizations into eq.\eqref{eq:biph-cont-int-linearized}, we obtain the stiffness matrix for biphasic stick as

\[ \begin{equation} \begin{aligned}D\delta G_{c} & =\sum_{a}\begin{bmatrix}\delta\mathbf{v}_{a}^{(1)} & \delta p_{a}^{(1)}\end{bmatrix}\cdot\int_{\Gamma_{\eta}^{(1)}}\left(\sum_{c}\begin{bmatrix}\mathbf{K}_{ac}^{(1,1)} & \mathbf{0}\\ \mathbf{g}_{ac}^{(1,1)} & g_{ac}^{(1,1)} \end{bmatrix}\cdot\begin{bmatrix}\Delta\mathbf{u}_{c}^{(1)}\\ \Delta p_{c}^{(1)} \end{bmatrix}\right.\\ & \qquad\qquad\qquad\left.+\sum_{d}\begin{bmatrix}\mathbf{K}_{ad}^{(1,2)} & \mathbf{0}\\ \mathbf{g}_{ad}^{(1,2)} & g_{ad}^{(1,2)} \end{bmatrix}\cdot\begin{bmatrix}\Delta\mathbf{u}_{d}^{(2)}\\ \Delta p_{d}^{(2)} \end{bmatrix}\right)J_{\eta}^{(1)}d\eta_{(1)}^{1}d\eta_{(1)}^{2}\\ & +\sum_{b}\begin{bmatrix}\delta\mathbf{v}_{b}^{(2)} & \delta p_{b}^{(2)}\end{bmatrix}\cdot\int_{\Gamma_{\eta}^{(1)}}\left(\sum_{c}\begin{bmatrix}\mathbf{K}_{bc}^{(2,1)} & \mathbf{0}\\ \mathbf{g}_{bc}^{(2,1)} & g_{bc}^{(2,1)} \end{bmatrix}\cdot\begin{bmatrix}\Delta\mathbf{u}_{c}^{(1)}\\ \Delta p_{c}^{(1)} \end{bmatrix}\right.\\ & \qquad\qquad\qquad\left.+\sum_{d}\begin{bmatrix}\mathbf{K}_{bd}^{(2,2)} & \mathbf{0}\\ \mathbf{g}_{bd}^{(2,2)} & g_{bd}^{(2,2)} \end{bmatrix}\cdot\begin{bmatrix}\Delta\mathbf{u}_{d}^{(2)}\\ \Delta p_{d}^{(2)} \end{bmatrix}\right)J_{\eta}^{(1)}d\eta_{(1)}^{1}d\eta_{(1)}^{2} \end{aligned} \label{eq:biphasic-stick-stiffness} \end{equation} \]

where

\[ \begin{equation} \begin{aligned}\mathbf{K}_{ac}^{(1,1)} & =-\varepsilon_{n}N_{a}^{(1)}N_{c}^{(1)}\mathbf{I}+N_{a}^{(1)}\mathbf{t}^{(1)}\otimes\mathbf{A}_{c}^{(1)}\cdot\mathbf{n}^{(1)}\\ \mathbf{K}_{ad}^{(1,2)} & =\varepsilon_{n}N_{a}^{(1)}N_{d}^{(2)}\mathbf{I}\\ \mathbf{K}_{bc}^{(2,1)} & =\varepsilon_{n}N_{b}^{(2)}N_{c}^{(1)}\mathbf{I}-N_{b}^{(2)}\mathbf{t}^{(1)}\otimes\mathbf{A}_{c}^{(1)}\cdot\mathbf{n}^{(1)}\\ \mathbf{K}_{bd}^{(2,2)} & =-\varepsilon_{n}N_{b}^{(2)}N_{d}^{(2)}\mathbf{I} \end{aligned} \end{equation} \]

and

\[ \begin{equation} \begin{aligned}\mathbf{g}_{ac}^{(1,1)} & =N_{a}^{(1)}w_{n}\mathbf{A}_{c}^{(1)}\cdot\mathbf{n}^{(1)} & g_{ac}^{(1,1)} & =\varepsilon_{p}N_{a}^{(1)}N_{c}^{(1)}\\ \mathbf{g}_{ad}^{(1,2)} & =\mathbf{0} & g_{ad}^{(1,2)} & =-\varepsilon_{p}N_{a}^{(1)}N_{d}^{(2)}\\ \mathbf{g}_{bc}^{(2,1)} & =-N_{b}^{(2)}w_{n}\mathbf{A}_{c}^{(1)}\cdot\mathbf{n}^{(1)} & g_{bc}^{(2,1)} & =-\varepsilon_{p}N_{b}^{(2)}N_{c}^{(1)}\\ \mathbf{g}_{bd}^{(2,2)} & =\mathbf{0} & g_{bd}^{(2,2)} & =\varepsilon_{p}N_{b}^{(2)}N_{d}^{(2)} \end{aligned} \end{equation} \]

In eq.\eqref{eq:biphasic-stick-stiffness}, \(\begin{bmatrix}\Delta\mathbf{u}_{c}^{(1)} & \Delta p_{c}^{(1)}\end{bmatrix}^{T}\) is the vector of incremental changes in the degrees of freedom of the \(c\)th node of the current element face on \(\gamma^{(1)}\). Similarly, \(\begin{bmatrix}\Delta\mathbf{u}_{d}^{(2)} & \Delta p_{d}^{(2)}\end{bmatrix}^{T}\) represents the incremental changes in the degrees of freedom of the \(d\)th node of the element face on \(\gamma^{(2)}\) which contains the intersection point \(X^{(2)}\) associated with the \(k\)th integration point on the current element face on \(\gamma^{(1)}\) ¹⁴. See Section Integration Scheme for a description of the Gaussian quadrature integration scheme used to evaluate residuals and stiffness matrices (e.g. Eqs.\eqref{eq:biph-stick-int-disc} and \eqref{eq:biphasic-stick-stiffness}). Briefly, it should be noted that for all terms associated with \(\gamma^{(2)}\) (e.g. \(\mathbf{K}_{ad}^{(1,2)}\), \(\Delta\mathbf{u}_{d}^{(2)}\), \(\delta p_{b}^{(2)}\)) there may be \(k\in\left[1,\,n_{int}^{(e)}\right]\) distinct element faces on \(\gamma^{(2)}\) associated with all the integration points \(X^{(1)}\) on the \(e\)th element face of \(\gamma^{(1)}\), based on the location of \(X^{(2)}\) obtained from eq.\eqref{eq:stick-x2}.

In Section Sliding-Elastic on sliding-elastic frictional contact, we split the contact stiffness matrices in a different way, which allowed us to clearly separate like terms. Here, due to the complexity of the biphasic contact formulation, it was determined that the form of eq.\eqref{eq:biphasic-stick-stiffness} provides more clarity.

Biphasic Slip

Linearization

During slip, the contact integral over \(\gamma^{(1)}\) is performed over integration points \(X^{(1)}\) with prescribed parametric coordinates \(\eta_{(1)}^{\alpha}\). However, the point on \(\gamma^{(2)}\) in contact with \(\gamma^{(1)}\) has parametric coordinates \(\eta_{(2)}^{\alpha}\) which change with variations in \(\mathbf{x}^{(1)}\) and \(\mathbf{n}^{(1)}\), in accordance with eq.\eqref{eq:slip-gap}. Consequently, directional derivatives of \(\mathbf{x}^{(i)}\), \(p^{(i)}\), \(\delta\mathbf{v}^{(i)}\), and \(\delta p^{(i)}\) are given by

\[ \begin{equation} \begin{aligned}D\mathbf{x}^{(1)} & =\Delta\mathbf{u}^{(1)} & D\mathbf{x}^{(2)} & =\Delta\mathbf{u}^{(2)}+\mathbf{g}_{\alpha}^{(2)}D\eta_{(2)}^{\alpha}\\ Dp^{(1)} & =\Delta p^{(1)} & Dp^{(2)} & =\Delta p^{(2)}+\frac{\partial p^{(2)}}{\partial\eta_{(2)}^{\alpha}}D\eta_{(2)}^{\alpha}\\ D\delta\mathbf{v}^{(1)} & =\mathbf{0} & D\delta\mathbf{v}^{(2)} & =\frac{\partial\delta\mathbf{v}^{(2)}}{\partial\eta_{(2)}^{\alpha}}D\eta_{(2)}^{\alpha}\\ D\delta p^{(1)} & =0 & D\delta p^{(2)} & =\frac{\partial\delta p^{(2)}}{\partial\eta_{(2)}^{\alpha}}D\eta_{(2)}^{\alpha} \end{aligned} \label{eq:biphasic-slip-linearization} \end{equation} \]

where \(D\eta_{(2)}^{\alpha}\) is evaluated by our modification ¹ of a method proposed by Laursen and Simo ¹¹, and may be found in Section Frictionless Terms. The linearization of the slip traction proceeds as before ⁴, with the addition of a term involving the linearization of \(\mueff\). From eq.\eqref{eq:mueff-biphasic} it follows that

\[ \begin{equation} D\mueff=\mueq\left(1-\varphi\right)\frac{\Delta p^{(1)}}{t_{n}}-\mueq\left(1-\varphi\right)\frac{p^{(1)}}{t_{n}^{2}}Dt_{n}\label{eq:Dmueff} \end{equation} \]

where \(Dt_{n}=\varepsilon Dg\) according to eq.\eqref{eq:slip-tn}; this term has been provided previously in eq.(7.1-77) ⁴. Equation \eqref{eq:Dmueff} was derived by recalling that \(\mueq\) and \(\varphi\) are constants. We then obtain

\[ \begin{equation} D\mathbf{t}^{(1)}=Dt_{n}\left(\no+\mueff\so\right)+t_{n}\left(D\no+\left(D\mueff\right)\so+\mueff D\so\right) \end{equation} \]

Finally, from eq.\eqref{eq:wn} (evaluated at \(\mathbf{x}^{\left(2\right)}\) determined by eq.\eqref{eq:slip-x2}) and eq.\eqref{eq:biphasic-slip-linearization} it follows that

\[ \begin{equation} Dw_{n}=\varepsilon_{p}\left(\Delta p^{(1)}-\Delta p^{(2)}-\frac{\partial p^{(2)}}{\partial\eta_{(2)}^{\alpha}}D\eta_{(2)}^{\alpha}\right)\label{eq:Dwn} \end{equation} \]

Note that the form of \(Dw_{n}\) given in eq.\eqref{eq:Dwn} contains additional terms not present in eq.\eqref{eq:biph-stick-Dwn}. This is because the parametric coordinates of intersection may vary in slip, but are invariant in stick; as a consequence, whether or not the linearization of \(p^{(2)}\) depends on \(\eta_{(2)}^{\alpha}\) is determined by the stick-slip status. Per the discussion following eq.\eqref{eq:pi-biphasic-penalty} in the main text, the fluid flux must be calculated after the stick-slip status has been resolved.

Discretization

Let the continuous variables on the primary and secondary surfaces be interpolated over each element face according to

\[ \begin{equation} \begin{aligned}\delta\mathbf{v}^{(1)} & =\sum_{a}N_{a}^{(1)}\delta\mathbf{v}_{a}^{(1)} & \delta\mathbf{v}^{(2)} & =\sum_{b}N_{b}^{(2)}\delta\mathbf{v}_{b}^{(2)}\\ \Delta\mathbf{u}^{(1)} & =\sum_{c}N_{c}^{(1)}\Delta\mathbf{u}_{c}^{(1)} & \Delta\mathbf{u}^{(2)} & =\sum_{d}N_{d}^{(2)}\Delta\mathbf{u}_{d}^{(2)}\\ \delta p^{(1)} & =\sum_{a}N_{a}^{(1)}\delta p_{a}^{(1)} & \delta p^{(2)} & =\sum_{b}N_{b}^{(2)}\delta p_{b}^{(2)}\\ \Delta p^{(1)} & =\sum_{c}N_{c}^{(1)}\Delta p_{c}^{(1)} & \Delta p^{(2)} & =\sum_{d}N_{d}^{(2)}\Delta p_{d}^{(2)} \end{aligned} \label{eq:biphasic-slip-var-discretization} \end{equation} \]

In the case of slip, the contact integral of eq.\eqref{eq:contact-int-invariant-1} can be split into normal and tangential parts, \(\delta G_{c}=\delta G_{c}^{n}+\delta G_{c}^{t}\), such that

\[ \begin{equation} \begin{aligned}\delta G_{c}= & \int_{\Gamma_{\eta}^{(1)}}\left(\begin{bmatrix}\delta\mathbf{v}^{(1)} & \delta\mathbf{v}^{(2)}\end{bmatrix}\cdot\begin{bmatrix}t_{n}\mathbf{n}^{(1)}\\ -t_{n}\mathbf{n}^{(1)} \end{bmatrix}+\begin{bmatrix}\delta p^{(1)} & \delta p^{(2)}\end{bmatrix}\cdot\begin{bmatrix}w_{n}\\ -w_{n} \end{bmatrix}\right)J_{\eta}^{(1)}d\eta_{(1)}^{1}d\eta_{(1)}^{2}\\ & +\int_{\Gamma_{\eta}^{(1)}}\begin{bmatrix}\delta\mathbf{v}^{(1)} & \delta\mathbf{v}^{(2)}\end{bmatrix}\cdot\begin{bmatrix}\mueff t_{n}\mathbf{s}^{(1)}\\ -\mueff t_{n}\mathbf{s}^{(1)} \end{bmatrix}\,J_{\eta}^{(1)}d\eta_{(1)}^{1}d\eta_{(1)}^{2} \end{aligned} \label{eq:biphasic-slip-discretized} \end{equation} \]

This split will be useful for the full linearization. Discretizing this expression yields

\[ \begin{equation} \begin{aligned}\delta G_{c} & =\sum_{a}\begin{bmatrix}\delta\mathbf{v}_{a}^{(1)} & \delta p_{a}^{(1)}\end{bmatrix}\cdot\int_{\Gamma_{\eta}^{(1)}}\begin{bmatrix}\mathbf{f}_{a}^{(1)}\\ w_{a}^{(1)} \end{bmatrix}J_{\eta}^{(1)}d\eta_{(1)}^{1}d\eta_{(1)}^{2}\\ & +\sum_{b}\begin{bmatrix}\delta\mathbf{v}_{b}^{(2)} & \delta p_{b}^{(2)}\end{bmatrix}\cdot\int_{\Gamma_{\eta}^{(1)}}\begin{bmatrix}\mathbf{f}_{b}^{(2)}\\ w_{b}^{(2)} \end{bmatrix}J_{\eta}^{(1)}d\eta_{(1)}^{1}d\eta_{(1)}^{2} \end{aligned} \end{equation} \]

where

\[ \begin{equation} \begin{aligned}\mathbf{f}_{a}^{(1)}= & N_{a}^{(1)}\mathbf{t}^{(1)}\,, & \mathbf{f}_{b}^{(2)}= & -N_{b}^{(2)}\mathbf{t}^{(1)}\\ w_{a}^{(1)}= & N_{a}^{(1)}w_{n}\,\,, & w_{b}^{(2)}= & -N_{b}^{(2)}w_{n} \end{aligned} \end{equation} \]

and \(\mathbf{t}^{(1)}\) is given by Eqs.\eqref{eq:slip-tn}-\eqref{eq:slip-traction} for the penalty method and Eqs.\eqref{eq:AL-tn}-\eqref{eq:augmented-t} for augmented Lagrangian regularization; similarly, \(w_{n}\) is given by eq.\eqref{eq:wn} or eq.\eqref{eq:AL-wn} for penalty and augmented Lagrangian schemes, respectively. For the following linearization and discretization the normal and tangential components, representing frictionless and frictional contributions to the contact integral, will be treated separately and may then be added together as in eq.\eqref{eq:biphasic-slip-discretized}.

Frictionless Terms

The linearization of the frictional part of eq.\eqref{eq:biphasic-slip-discretized} makes use of eq.\eqref{eq:biphasic-slip-linearization} to find

\[ \begin{equation} \begin{aligned}D\delta G_{c} & =\int_{\Gamma_{\eta}^{(1)}}\begin{bmatrix}\delta\mathbf{v}^{(1)} & \delta p^{(1)}\end{bmatrix}\cdot\left(\begin{bmatrix}t_{n}\mathbf{n}^{(1)}\\ w_{n} \end{bmatrix}\frac{1}{J_{\eta}^{(1)}}DJ_{\eta}^{(1)}+\begin{bmatrix}t_{n}\\ 0 \end{bmatrix}D\mathbf{n}^{(1)}+\begin{bmatrix}\mathbf{n}^{(1)}\\ 0 \end{bmatrix}Dt_{n}\right.\\ & \qquad\qquad\left.+\begin{bmatrix}0\\ 1 \end{bmatrix}Dw_{n}\right)J_{\eta}^{(1)}d\eta_{(1)}^{1}d\eta_{(1)}^{2}\\ & +\int_{\Gamma_{\eta}^{(1)}}\begin{bmatrix}\frac{\partial\delta\mathbf{v}^{(2)}}{\partial\eta_{(2)}^{\alpha}} & \frac{\partial\delta p^{(2)}}{\partial\eta_{(2)}^{\alpha}}\end{bmatrix}\cdot\begin{bmatrix}-t_{n}\mathbf{n}^{(1)}\\ -w_{n} \end{bmatrix}D\eta_{(2)}^{\alpha}\,J_{\eta}^{(1)}d\eta_{(1)}^{1}d\eta_{(1)}^{2}\\ & +\int_{\Gamma_{\eta}^{(1)}}\begin{bmatrix}\delta\mathbf{v}^{(2)} & \delta p^{(2)}\end{bmatrix}\cdot\left(\begin{bmatrix}-t_{n}\mathbf{n}^{(1)}\\ -w_{n} \end{bmatrix}\frac{1}{J_{\eta}^{(1)}}DJ_{\eta}^{(1)}+\begin{bmatrix}-t_{n}\\ 0 \end{bmatrix}D\mathbf{n}^{(1)}+\begin{bmatrix}-\mathbf{n}^{(1)}\\ 0 \end{bmatrix}Dt_{n}\right.\\ & \qquad\qquad\left.\begin{bmatrix}0\\ -1 \end{bmatrix}Dw_{n}\right)J_{\eta}^{(1)}d\eta_{(1)}^{1}d\eta_{(1)}^{2} \end{aligned} \label{eq:Dgc-biph-fricless-slip} \end{equation} \]

where we note that the virtual variables on the secondary surface now enter the linearization, as parametric coordinates on the secondary surface vary during slip according to eq.\eqref{eq:slip-gap}. The only linearization which has not been previously discretized is \(Dw_{n}\), and it follows from placing eq.\eqref{eq:biphasic-slip-var-discretization} into eq.\eqref{eq:Dwn} that

\[ \begin{equation} Dw_{n}=\sum_{c}\begin{bmatrix}-\varepsilon_{p}\Nc\mathbf{p}^{(2)}+\varepsilon_{p}P_{c}\no & \varepsilon_{p}\Nc\end{bmatrix}\cdot\begin{bmatrix}\Delta\mathbf{u}_{c}^{(1)}\\ \Delta p_{c}^{(1)} \end{bmatrix}+\sum_{d}\begin{bmatrix}\varepsilon_{p}\Nd\mathbf{p}^{(2)} & -\varepsilon_{p}\Nd\end{bmatrix}\cdot\begin{bmatrix}\Delta\mathbf{u}_{d}^{(2)}\\ \Delta p_{d}^{(2)} \end{bmatrix}\label{eq:Dwn-discrete} \end{equation} \]

where we define

\[ \begin{equation} \begin{aligned}\mathbf{p}^{(2)} & =\frac{\partial p^{(2)}}{\partial\eta_{(2)}^{\alpha}}\bar{\mathbf{g}}_{(2)}^{\alpha}\\ P_{c} & =ga^{\alpha\beta}\frac{\partial p^{(2)}}{\partial\eta_{(2)}^{\alpha}}\frac{\partial\Nc}{\partial\eta_{(1)}^{\beta}} \end{aligned} \label{eq:def-p-pc} \end{equation} \]

and for convenience we note that

\[ \begin{equation} \frac{\partial p^{(2)}}{\partial\eta_{(2)}^{\alpha}}D\eta_{(2)}^{\alpha}=\sum_{c}\begin{bmatrix}\Nc\mathbf{p}^{(2)}-P_{c}\no & 0\end{bmatrix}\cdot\begin{bmatrix}\Delta\mathbf{u}_{c}^{(1)}\\ \Delta p_{c}^{(1)} \end{bmatrix}+\sum_{d}\begin{bmatrix}-\Nd\mathbf{p}^{(2)} & 0\end{bmatrix}\cdot\begin{bmatrix}\Delta\mathbf{u}_{d}^{(2)}\\ \Delta p_{d}^{(2)} \end{bmatrix}\label{eq:dp-deta} \end{equation} \]

Directional derivatives of virtual variables on the secondary surface will lead to expressions of the form \(\left(\partial\Nb/\partial\eta_{(2)}^{\lambda}\right)D\eta_{(2)}^{\lambda}\), where

\[ \begin{equation} \frac{\partial N_{b}^{(2)}}{\partial\eta_{(2)}^{\lambda}}D\eta_{(2)}^{\lambda}=\begin{bmatrix}\sum_{c}N_{c}^{(1)}\bar{\mathbf{m}}_{b}^{(2)}-G_{bc}\mathbf{n}^{(1)} & \sum_{d}-N_{d}^{(2)}\bar{\mathbf{m}}_{b}^{(2)}\end{bmatrix}\cdot\begin{bmatrix}\Delta\mathbf{u}_{c}^{(1)}\\ \Delta\mathbf{u}_{d}^{(2)} \end{bmatrix}\label{eq:dNb-deta} \end{equation} \]

In eq.\eqref{eq:dNb-deta}, the pressure degrees of freedom do not enter into the linearization. In an effort to keep the expression more compact, we have not included these variables. However, this expression could easily be cast into the form of e.g. eq.\eqref{eq:dp-deta} by adding zeros where necessary. Discretizing eq.\eqref{eq:Dgc-biph-fricless-slip} and making use of linearizations found previously ⁴, along with Eqs.\eqref{eq:Dwn-discrete} and \eqref{eq:dp-deta}, allows the resulting stiffness matrix to be expressed as

\[ \begin{equation} \begin{aligned}D\delta G_{c} & =\sum_{a}\begin{bmatrix}\delta\mathbf{v}_{a}^{(1)} & \delta p_{a}^{(1)}\end{bmatrix}\cdot\int_{\Gamma_{\eta}^{(1)}}\left(\sum_{c}\begin{bmatrix}\mathbf{K}_{ac}^{(1,1)} & \mathbf{0}\\ \mathbf{g}_{ac}^{(1,1)} & g_{ac}^{(1,1)} \end{bmatrix}\cdot\begin{bmatrix}\Delta\mathbf{u}_{c}^{(1)}\\ \Delta p_{c}^{(1)} \end{bmatrix}\right.\\ & \qquad\qquad\qquad\left.+\sum_{d}\begin{bmatrix}\mathbf{K}_{ad}^{(1,2)} & \mathbf{0}\\ \mathbf{g}_{ad}^{(1,2)} & g_{ad}^{(1,2)} \end{bmatrix}\cdot\begin{bmatrix}\Delta\mathbf{u}_{d}^{(2)}\\ \Delta p_{d}^{(2)} \end{bmatrix}\right)J_{\eta}^{(1)}d\eta_{(1)}^{1}d\eta_{(1)}^{2}\\ & +\sum_{b}\begin{bmatrix}\delta\mathbf{v}_{b}^{(2)} & \delta p_{b}^{(2)}\end{bmatrix}\cdot\int_{\Gamma_{\eta}^{(1)}}\left(\sum_{c}\begin{bmatrix}\mathbf{K}_{bc}^{(2,1)} & \mathbf{0}\\ \mathbf{g}_{bc}^{(2,1)} & g_{bc}^{(2,1)} \end{bmatrix}\cdot\begin{bmatrix}\Delta\mathbf{u}_{c}^{(1)}\\ \Delta p_{c}^{(1)} \end{bmatrix}\right.\\ & \qquad\qquad\qquad\left.+\sum_{d}\begin{bmatrix}\mathbf{K}_{bd}^{(2,2)} & \mathbf{0}\\ \mathbf{g}_{bd}^{(2,2)} & g_{bd}^{(2,2)} \end{bmatrix}\cdot\begin{bmatrix}\Delta\mathbf{u}_{d}^{(2)}\\ \Delta p_{d}^{(2)} \end{bmatrix}\right)J_{\eta}^{(1)}d\eta_{(1)}^{1}d\eta_{(1)}^{2} \end{aligned} \label{eq:biphasic-fricless-slip-stiffness} \end{equation} \]

where

\[ \begin{equation} \begin{aligned}\mathbf{K}_{ac}^{(1,1)} & =-N_{a}^{(1)}N_{c}^{(1)}\left(\varepsilon\Nt\right)-\Na\left(t_{n}\Ac+t_{n}\Mc\cdot\No\right)\\ \mathbf{K}_{ad}^{(1,2)} & =N_{a}^{(1)}N_{d}^{(2)}\left(\varepsilon\Nt\right)\\ \mathbf{K}_{bc}^{(2,1)} & =N_{b}^{(2)}N_{c}^{(1)}\left(\varepsilon\Nt\right)+N_{b}^{(2)}\left(t_{n}\Ac+t_{n}\Mc\cdot\No\right)\\ & \qquad\qquad+\Nc\left(t_{n}\Mb\right)+G_{bc}\left(t_{n}\No\right)\\ \mathbf{K}_{bd}^{(2,2)} & =-N_{b}^{(2)}N_{d}^{(2)}\left(\varepsilon\Nt\right)-\Nd\left(t_{n}\Mb\right) \end{aligned} \end{equation} \]

and

\[ \begin{equation} \begin{aligned}\mathbf{g}_{ac}^{(1,1)} & =-N_{a}^{(1)}N_{c}^{(1)}\left(\varepsilon_{p}\mathbf{p}^{(2)}\right)+N_{a}^{(1)}\left(w_{n}\Ac\cdot\no+\varepsilon_{p}P_{c}\no\right)\\ \mathbf{g}_{ad}^{(1,2)} & =N_{a}^{(1)}N_{d}^{(2)}\left(\varepsilon_{p}\mathbf{p}^{(2)}\right)\\ \mathbf{g}_{bc}^{(2,1)} & =N_{b}^{(2)}N_{c}^{(1)}\left(\varepsilon_{p}\mathbf{p}^{(2)}\right)-N_{b}^{(2)}\left(w_{n}\Ac\cdot\no+\varepsilon_{p}P_{c}\no\right)-N_{c}^{(1)}\left(w_{n}\mb\right)+G_{bc}\left(w_{n}\no\right)\\ \mathbf{g}_{bd}^{(2,2)} & =-N_{b}^{(2)}N_{d}^{(2)}\left(\varepsilon_{p}\mathbf{p}^{(2)}\right)+N_{d}^{(2)}\left(w_{n}\mb\right) \end{aligned} \end{equation} \]

and

\[ \begin{equation} \begin{aligned}g_{ac}^{(1,1)} & =-N_{a}^{(1)}N_{c}^{(1)}\left(-\varepsilon_{p}\right)\\ g_{ad}^{(1,2)} & =N_{a}^{(1)}N_{d}^{(2)}\left(-\varepsilon_{p}\right)\\ g_{bc}^{(2,1)} & =N_{b}^{(2)}N_{c}^{(1)}\left(-\varepsilon_{p}\right)\\ g_{bd}^{(2,2)} & =-N_{b}^{(2)}N_{d}^{(2)}\left(-\varepsilon_{p}\right) \end{aligned} \end{equation} \]

The expressions above are very similar to those which can be found in our frictionless biphasic contact paper ⁸, although the present framework is more general, as that previous study evaluated expressions in the limit as \(g\to0\) (see Zimmerman and Ateshian ⁴ for a detailed discussion of this assumption and the benefits of relaxing it).

Frictional Terms

Linearizing the frictional contribution follows from the second term of eq.\eqref{eq:biphasic-slip-discretized} as

\[ \begin{equation} \begin{aligned}D\delta G_{c} & =\int_{\Gamma_{\eta}^{(1)}}\begin{bmatrix}\delta\mathbf{v}^{(1)} & \delta p^{(1)}\end{bmatrix}\cdot\left(\begin{bmatrix}t_{n}\so\\ 0 \end{bmatrix}D\mueff+\begin{bmatrix}\mueff\so\\ 0 \end{bmatrix}Dt_{n}\right.\\ & \left.+\begin{bmatrix}\mueff t_{n}\\ 0 \end{bmatrix}D\so+\begin{bmatrix}\mueff t_{n}\so\\ 0 \end{bmatrix}\frac{1}{\jeta}D\jeta\right)\jeta d\eta_{(1)}^{1}d\eta_{(1)}^{2}\\ & +\int_{\Gamma_{\eta}^{(1)}}\begin{bmatrix}\delta\mathbf{v}^{(2)} & \delta p^{(2)}\end{bmatrix}\cdot\left(\begin{bmatrix}-t_{n}\so\\ 0 \end{bmatrix}D\mueff+\begin{bmatrix}-\mueff\so\\ 0 \end{bmatrix}Dt_{n}\right.\\ & \left.+\begin{bmatrix}-\mueff t_{n}\\ 0 \end{bmatrix}D\so+\begin{bmatrix}-\mueff t_{n}\so\\ 0 \end{bmatrix}\frac{1}{\jeta}D\jeta\right)\jeta d\eta_{(1)}^{1}d\eta_{(1)}^{2}\\ & +\int_{\Gamma_{\eta}^{(1)}}\begin{bmatrix}\frac{\partial\delta\mathbf{v}^{(2)}}{\partial\eta_{(2)}^{\alpha}} & \frac{\partial\delta p^{(2)}}{\partial\eta_{(2)}^{\alpha}}\end{bmatrix}\cdot\begin{bmatrix}-\mueff t_{n}\so\\ 0 \end{bmatrix}D\eta_{(2)}^{\alpha}\jeta d\eta_{(1)}^{1}d\eta_{(1)}^{2} \end{aligned} \label{eq:Dgc-biph-fric-slip} \end{equation} \]

The remaining quantity in this expression to be determined is the discretization of \(D\mueff\); from eq.\eqref{eq:Dmueff} and eq.\eqref{eq:biphasic-slip-var-discretization} it follows that

\[ \begin{equation} \begin{aligned}D\mueff & =\sum_{c}\begin{bmatrix}\varepsilon\Nc\mueq\Gamma_{p}\frac{1}{t_{n}^{2}}\Nbo\cdot\no+\mueq\Gamma_{p}\frac{1}{t_{n}}\No\cdot\mc & \Nc\mueq\left(1-\varphi\right)\frac{1}{t_{n}}\end{bmatrix}\cdot\begin{bmatrix}\Delta\mathbf{\mathbf{u}}_{c}^{(1)}\\ \Delta p_{c}^{(1)} \end{bmatrix}\\ & +\sum_{d}\begin{bmatrix}-\varepsilon_{n}\Nd\mueq\Gamma_{p}\frac{1}{t_{n}^{2}}\Nbo\cdot\no & 0\end{bmatrix}\cdot\begin{bmatrix}\Delta\mathbf{\mathbf{u}}_{d}^{(2)}\\ \Delta p_{d}^{(2)} \end{bmatrix} \end{aligned} \label{eq:Dmueff-discrete} \end{equation} \]

where we made the definition

\[ \begin{equation} \Gamma_{p}=\left(1-\varphi\right)p^{(1)} \end{equation} \]

just to be used in eq.\eqref{eq:Dmueff-discrete} for space considerations. Finally, inserting Eqs.\eqref{eq:biphasic-slip-var-discretization} and \eqref{eq:Dmueff-discrete} into eq.\eqref{eq:Dgc-biph-fric-slip} yields the stiffness matrix for the frictional terms

\[ \begin{equation} \begin{aligned}D\delta G_{c} & =\sum_{a}\begin{bmatrix}\delta\mathbf{v}_{a}^{(1)} & \delta p_{a}^{(1)}\end{bmatrix}\cdot\int_{\Gamma_{\eta}^{(1)}}\left(\sum_{c}\begin{bmatrix}\mathbf{K}_{ac}^{(1,1)} & \mathbf{k}_{ac}^{(1,1)}\\ \mathbf{0} & 0 \end{bmatrix}\cdot\begin{bmatrix}\Delta\mathbf{u}_{c}^{(1)}\\ \Delta p_{c}^{(1)} \end{bmatrix}\right.\\ & \qquad\qquad\qquad\left.+\sum_{d}\begin{bmatrix}\mathbf{K}_{ad}^{(1,2)} & \mathbf{k}_{ad}^{(1,2)}\\ \mathbf{0} & 0 \end{bmatrix}\cdot\begin{bmatrix}\Delta\mathbf{u}_{d}^{(2)}\\ \Delta p_{d}^{(2)} \end{bmatrix}\right)J_{\eta}^{(1)}d\eta_{(1)}^{1}d\eta_{(1)}^{2}\\ & +\sum_{b}\begin{bmatrix}\delta\mathbf{v}_{b}^{(2)} & \delta p_{b}^{(2)}\end{bmatrix}\cdot\int_{\Gamma_{\eta}^{(1)}}\left(\sum_{c}\begin{bmatrix}\mathbf{K}_{bc}^{(2,1)} & \mathbf{k}_{bc}^{(2,1)}\\ \mathbf{0} & 0 \end{bmatrix}\cdot\begin{bmatrix}\Delta\mathbf{u}_{c}^{(1)}\\ \Delta p_{c}^{(1)} \end{bmatrix}\right.\\ & \qquad\qquad\qquad\left.+\sum_{d}\begin{bmatrix}\mathbf{K}_{bd}^{(2,2)} & \mathbf{k}_{bd}^{(2,2)}\\ \mathbf{0} & 0 \end{bmatrix}\cdot\begin{bmatrix}\Delta\mathbf{u}_{d}^{(2)}\\ \Delta p_{d}^{(2)} \end{bmatrix}\right)J_{\eta}^{(1)}d\eta_{(1)}^{1}d\eta_{(1)}^{2} \end{aligned} \label{eq:biphasic-fric-slip-stiffness} \end{equation} \]

where

\[ \begin{equation} \begin{aligned}\mathbf{K}_{ac}^{(1,1)} & =-N_{a}^{(1)}N_{c}^{(1)}\left(\varepsilon\mueq\tilde{\mathbf{S}}^{(1)}+\mueff t_{n}\mathbf{B}^{(1)}\right)-N_{a}^{(1)}\left(t_{n}\so\otimes\bar{\mathbf{h}}_{c-}^{(1)}\right.\\ & \qquad\qquad\left.+\mueff t_{n}g\mathbf{P}_{S}\cdot\mathbf{c}^{(1)}\otimes\mathbf{h}_{c+}^{(1)}-\mueff t_{n}\mathbf{J}_{c}^{(1)}\right)\\ \mathbf{K}_{ad}^{(1,2)} & =N_{a}^{(1)}N_{d}^{(2)}\left(\varepsilon\mueq\tilde{\mathbf{S}}^{(1)}+\mueff t_{n}\mathbf{B}^{(1)}\right)\\ \mathbf{K}_{bc}^{(2,1)} & =N_{b}^{(2)}N_{c}^{(1)}\left(\varepsilon\mueq\tilde{\mathbf{S}}^{(1)}+\mueff t_{n}\mathbf{B}^{(1)}\right)+N_{b}^{(2)}\left(t_{n}\so\otimes\bar{\mathbf{h}}_{c-}^{(1)}\right.\\ & \qquad\qquad\left.+\mueff t_{n}g\mathbf{P}_{S}\cdot\mathbf{c}^{(1)}\otimes\mathbf{h}_{c+}^{(1)}-\mueff t_{n}\mathbf{J}_{c}^{(1)}\right)\\ & \qquad\qquad+\Nc\left(-\mueff t_{n}\so\otimes\mb\right)+G_{bc}\left(\mueff t_{n}\So\right)\\ \mathbf{K}_{bd}^{(2,2)} & =-N_{b}^{(2)}N_{d}^{(2)}\left(\varepsilon\mueq\tilde{\mathbf{S}}^{(1)}+\mueff t_{n}\mathbf{B}^{(1)}\right)-\Nd\left(-\mueff t_{n}\so\otimes\mb\right) \end{aligned} \label{eq:biphasic-fric-slip-solid-stiffness} \end{equation} \]

and

\[ \begin{equation} \begin{aligned}\mathbf{k}_{ac}^{(1,1)} & =-N_{a}^{(1)}N_{c}^{(1)}\left(-\mueq\left(1-\varphi\right)\so\right)\\ \mathbf{k}_{ad}^{(1,2)} & =\mathbf{0}\\ \mathbf{k}_{bc}^{(2,1)} & =N_{b}^{(2)}N_{c}^{(1)}\left(-\mueq\left(1-\varphi\right)\so\right)\\ \mathbf{k}_{bd}^{(2,2)} & =\mathbf{0} \end{aligned} \label{eq:biphasic-fric-slip-pressure-stiffness} \end{equation} \]

In eq.\eqref{eq:biphasic-fric-slip-solid-stiffness}, we have defined

\[ \begin{equation} \begin{aligned}\bar{\mathbf{h}}_{c-}^{(1)} & =\mueq\No\cdot\mc-\mueff\Ac\cdot\no\end{aligned} \label{eq:h-bar-c-minus} \end{equation} \]

by analogy with our definitions for \(\mathbf{h}_{c+}^{(1)}\) and \(\mathbf{h}_{c-}^{(1)}\) in Section Frictional Terms. The stiffness matrix of the frictional contribution to the virtual work, given by eq.\eqref{eq:biphasic-fric-slip-stiffness}, is nonsymmetric. Summing Eqs.\eqref{eq:biphasic-fricless-slip-stiffness} and \eqref{eq:biphasic-fric-slip-stiffness} produces the total stiffness matrix for the case of frictional biphasic slip; this total stiffness matrix is also nonsymmetric. However, as noted before ⁸, the stiffness matrix for biphasic materials is nonsymmetric by construction, so there is no expectation that the contact stiffness matrices be symmetric.

Equation \eqref{eq:biphasic-fric-slip-solid-stiffness} contains the solid-solid stiffness terms. Comparing eq.\eqref{eq:biphasic-fric-slip-solid-stiffness} with eq.(7.1-112) shows that the same form of the equations is recovered. Interestingly, however, some of the terms in sliding-elastic (Section Sliding-Elastic) which were multiplied by the friction coefficient \(\mu\) are now multiplied by \(\mueff\), which indicates certain terms which are important in biphasic contact. Of course, in the absence of pressurized fluid, \(\mueff\to\mu\) and we recover the sliding-elastic formulation.

---

1. Ateshian, GA; Maas, SA; Weiss, J.A.. "Finite element algorithm for frictionless contact of porous permeable media under finite deformation and sliding." J. Biomech. Engn., vol. 132, pp. 1006-1019 (2010). ↩↩↩

2. Zimmerman, Brandon K; Maas, Steve A; Weiss, Jeffrey A; Ateshian, Gerard A. "A Finite Element Algorithm for Large Deformation Biphasic Frictional Contact Between Porous-Permeable Hydrated Soft Tissues." J Biomech Eng, vol. 144 (2022). ↩↩

3. Bonet, Javier; Wood, Richard D.. "Nonlinear continuum mechanics for finite element analysis." Cambridge University Press (1997). ↩

4. Zimmerman, Brandon K; Ateshian, Gerard A. "A Surface-to-Surface Finite Element Algorithm for Large Deformation Frictional Contact in febio." J Biomech Eng, vol. 140 (2018). ↩↩↩↩↩↩↩↩↩↩↩↩↩↩↩↩↩↩↩

5. Ateshian, Gerard A. "The role of interstitial fluid pressurization in articular cartilage lubrication." J Biomech, vol. 42, pp. 1163-76 (2009). ↩

6. Wriggers, P; Laursen, Tod A. "Computational contact mechanics." Springer, vol. no. 498 (2007). ↩↩

7. Curnier, Alain; He, Qi-Chang; Klarbring, Anders. "Continuum mechanics modelling of large deformation contact with friction." Contact mechanics, pp. 145--158 (1995). ↩

8. Ateshian, G. A.; Weiss, J. A.. "Anisotropic hydraulic permeability under finite deformation." Journal of biomechanical engineering, vol. 132, pp. 111004 (2010). ↩↩↩↩↩

9. Simo, J.C.; Armero, F.. "Geometrically Non-linear Enhanced Strain Mixed Methods and the Method of Incompatible Modes." International Journal for Numerical Methods in Engineering, vol. 33, pp. 1413-1419 (1992). ↩

10. Bertsekas, Dimitri P. "Constrained optimization and Lagrange multiplier methods." Academic Press (1982). ↩

11. Laursen, T. A.; Simo, J. C.. "Continuum-based finite element formulation for the implicit solution of multibody, large deformation frictional contact problems." International Journal for Numerical Methods in Engineering, vol. 36, pp. 3451-3485 (1993). ↩