The Galerkin weighted residual method has been applied to this strong form to derive alternate velocity-pressure finite element contact formulations. The Lagrange multiplier method is used, and enforces two of the four contact continuity conditions, while the other two conditions are introduced directly into the weighted residual statement. The alternate formulations differ in the choice of continuity conditions to enforce with Lagrange multipliers. In one they enforce the normal solid traction and relative fluid flow continuity conditions on the contact surface, and in the other the multipliers enforce normal solid traction and pressure continuity conditions. The contact nonlinearity is treated with an iterative algorithm, where the assumed area is either extended or reduced based on the validity of the solution relative to contact conditions such as impenetrability and intensility. The resulting first order system of equations is solved in time using the generalized finite difference scheme.

In an independent study, various combinations of the Krylov iterative solvers have been tested with alternate equation re-ordering schemes and levels of preconditioning. An extensive study has been performed using the 3D linear biphasic equations, without contact, the results of which have provided guidance for efficiently solving the symmetric indefinite system found in the present v-p contact formulations.

A biphasic contact patch test has been used to select interpolation functions for the Lagrange multipliers and to confirm that the formulations are capable of transmitting the constant normal traction as a finite element completeness check. Then, the preferred formulation has been validated by a series of increasingly complex canonical problems including the confined and unconfined compression tests, the Hertz contact problem and two biphasic indentation tests. As a clinical demonstration of the capability of the contact analysis, the gleno-humeral joint contact of human shoulders has been analyzed using an idealized 3D geometry. In the joint, both glenoid and humeral head cartilage experience maximum tensile and compressive stresses at the cartilage-bone interface, away from the center of the contact area.