0

### TECHNICAL PAPERS: Joint/Whole Body

J Biomech Eng. 2006;129(3):386-392. doi:10.1115/1.2720915.

The knee joint is partially stabilized by the interaction of multiple ligament structures. This study tested the interdependent functions of the anterior cruciate ligament (ACL) and the medial collateral ligament (MCL) by evaluating the effects of ACL deficiency on local MCL strain while simultaneously measuring joint kinematics under specific loading scenarios. A structural testing machine applied anterior translation and valgus rotation (limits $100N$ and $10Nm$, respectively) to the tibia of ten human cadaveric knees with the ACL intact or severed. A three-dimensional motion analysis system measured joint kinematics and MCL tissue strain in 18 regions of the superficial MCL. ACL deficiency significantly increased MCL strains by 1.8% $(p<0.05)$ during anterior translation, bringing ligament fibers to strain levels characteristic of microtrauma. In contrast, ACL transection had no effect on MCL strains during valgus rotation (increase of only 0.1%). Therefore, isolated valgus rotation in the ACL-deficient knee was nondetrimental to the MCL. The ACL was also found to promote internal tibial rotation during anterior translation, which in turn decreased strains near the femoral insertion of the MCL. These data advance the basic structure-function understanding of the MCL, and may benefit the treatment of ACL injuries by improving the knowledge of ACL function and clarifying motions that are potentially harmful to secondary stabilizers.

Commentary by Dr. Valentin Fuster
J Biomech Eng. 2006;129(3):393-399. doi:10.1115/1.2737432.

Because fall experiments with volunteers can be both challenging and risky, especially with older volunteers, we wished to develop computer simulations of falls to provide a theoretical framework for understanding and extending experimental results. To perform a preliminary validation of the articulated total body (ATB) model for passive falls, we compared the model predictions of fall direction, impact location, and impact velocity as a function of disturbance type (faint, slip, step down, trip) and gait speed (fast, normal, slow) to experimental results with young adult volunteers. The three-dimensional ATB model had 17 segments and 16 joints. Its physical characteristics, environment definitions, contact functions, and initial conditions were representative of our experiment. For each combination of disturbance and gait speed, the ATB model was left to fall passively under gravity once disturbed, i.e., no joint torques were applied, until impact with the floor occurred. Finally, we also determined the sensitivity of the model predictions to changes in the model’s parameters. Our model predictions of fall angles and impact angles were qualitatively in agreement with those observed experimentally for ten and seven of the 12 original simulations, respectively. Quantitatively, the model predictions of fall angles, impact angles, and impact velocities were within one experimental standard deviation for seven, three, and nine of the 12 original simulations, respectively, and within two experimental standard deviations for ten, nine, and 11 of the 12 original simulations, respectively. Finally, the fall angle and impact angle region did not change for 92% and 95% of the 74 input variation simulations, respectively, and the impact velocities were within the experimental standard deviations for 78% of the 74 input variation simulations. Based on our simulations and a sensitivity analysis, we conclude that our preliminary validation of the ATB model for passive falls was successful. In fact, these ATB model simulations represent a significant step forward in fall simulations. We believe that with additional work, the ATB model could be used to accurately simulate a variety of human falls and may be useful in further understanding the etiology and mechanisms of fall injuries such as hip fractures.

Commentary by Dr. Valentin Fuster

### TECHNICAL PAPERS: Fluids/Heat/Transport

J Biomech Eng. 2006;129(3):330-340. doi:10.1115/1.2720910.

Electroporation is an approach used to enhance transdermal transport of large molecules in which the skin is exposed to a series of electric pulses. The structure of the transport inhibiting outer layer, the stratum corneum, is temporarily destabilized due to the development of microscopic pores. Consequently agents that are ordinarily unable to pass into the skin are able to pass through this outer barrier. Of possible concern when exposing biological tissue to an electric field is thermal tissue damage associated with Joule heating. This paper shows the importance of using a composite model in calculating the electrical and thermal effects associated with skin electroporation. A three-dimensional transient finite-volume model of in vivo skin electroporation is developed to emphasize the importance of representing the skin’s composite layers and to illustrate the underlying relationships between the physical parameters of the composite makeup of the skin and resulting thermal damage potential.

Commentary by Dr. Valentin Fuster
J Biomech Eng. 2006;129(3):341-353. doi:10.1115/1.2720911.

Expiratory droplets and droplet nuclei can be pathogen carriers for airborne diseases. Their transport characteristics were studied in detail in two idealized floor-supply-type ventilation flow patterns: Unidirectional–upward and single-side–floor, using a multiphase numerical model. The model was validated by running interferometric Mie imaging experiments using test droplets with nonvolatile content, which formed droplet nuclei, ultimately, in a class-100 clean-room chamber. By comparing the droplet dispersion and removal characteristics with data of two other ceiling-supply ventilation systems collected from a previous work, deviations from the perfectly mixed ventilation condition were found to exist in various cases to different extent. The unidirectional–upward system was found to be more efficient in removing the smallest droplet nuclei (formed from $1.5μm$ droplets) by air extraction, but it became less effective for larger droplets and droplet nuclei. Instead, the single-side–floor system was shown to be more favorable in removing these large droplets and droplet nuclei. In the single-side–floor system, the lateral overall dispersion coefficients for the small droplets and nuclei (initial size $⩽45μm$) were about an order of magnitude higher than those in the unidirectional–upward system. It indicated that bulk lateral airflow transport in the single-side–floor system was much stronger than the lateral dispersion mechanism induced mainly by air turbulence in the unidirectional–upward system. The time required for the droplets and droplet nuclei to be transported to the exhaust vent or deposition surfaces for removal varied with different ventilation flow patterns. Possible underestimation of exposure level existed if the perfectly mixed condition was assumed. For example, the weak lateral dispersion in the unidirectional ventilation systems made expiratory droplets and droplet nuclei stay at close distance to the source leading to highly nonuniform spatial distributions. The distance between the source and susceptible patients became an additional concern in exposure analysis. Relative significance of the air-extraction removal mechanism was studied. This can have impact to the performance evaluation of filtration and disinfection systems installed in the indoor environment. These findings revealed the need for further development in a risk-assessment model incorporating the effect of different ventilation systems on distributing expiratory droplets and droplet nuclei nonuniformly in various indoor spaces, such as buildings, aircraft cabins, trains, etc.

Commentary by Dr. Valentin Fuster
J Biomech Eng. 2006;129(3):354-364. doi:10.1115/1.2720912.

This study uses a reconstructed vascular geometry to evaluate the thermal response of tissue during a three-dimensional radiofrequency (rf) tumor ablation. MRI images of a sectioned liver tissue containing arterial vessels are processed and converted into a finite-element mesh. A rf heat source in the form of a spherically symmetric Gaussian distribution, fit from a previously computed profile, is employed. Convective cooling within large blood vessels is treated using direct physical modeling of the heat and momentum transfer within the vessel. Calculations of temperature rise and thermal dose are performed for transient rf procedures in cases where the tumor is located at three different locations near the bifurcation point of a reconstructed artery. Results demonstrate a significant dependence of tissue temperature profile on the reconstructed vasculature and the tumor location. Heat convection through the arteries reduced the steady-state temperature rise, relative to the no-flow case, by up to 70% in the targeted volume. Blood flow also reduced the thermal dose value, which quantifies the extent of cell damage, from $∼3600min$, for the no-flow condition, to $10min$ for basal flow $(13.8cm∕s)$. Reduction of thermal dose below the threshold value of $240min$ indicates ablation procedures that may inadequately elevate the temperature in some regions, thereby permitting possible tumor recursion. These variations are caused by vasculature tortuosity that are patient specific and can be captured only by the reconstruction of the realistic geometry.

Commentary by Dr. Valentin Fuster
J Biomech Eng. 2006;129(3):365-373. doi:10.1115/1.2720913.

Microchannel bioreactors have applications for manipulating and investigating the fluid microenvironment on cell growth and functions in either single culture or co-culture. This study considers two different types of cells distributed randomly as a co-culture at the base of a microchannel bioreactor: absorption cells, which only consume species based on the Michaelis-Menten process, and release cells, which secrete species, assuming zeroth order reaction, to support the absorption cells. The species concentrations at the co-culture cell base are computed from a three-dimensional numerical flow-model incorporating mass transport. Combined dimensionless parameters are proposed for the co-culture system, developed from a simplified analysis under the condition of decreasing axial-concentration. The numerical results of species concentration at the co-culture cell-base are approximately correlated by the combined parameters under the condition of positive flux-parameter. Based on the correlated results, the critical value of the inlet concentration is determined, which depends on the effective microchannel length. For the flow to develop to the critical inlet concentration, an upstream length consisting only of release cells is needed; this upstream length is determined from an analytical solution. The generalized results may find applications in analyzing the mass transport requirements in a co-culture microchannel bioreactor.

Commentary by Dr. Valentin Fuster
J Biomech Eng. 2006;129(3):374-385. doi:10.1115/1.2720914.

Atherosclerosis localizes at a bend and∕or bifurcation of an artery, and low density lipoproteins (LDL) accumulate in the intima. Hemodynamic factors are known to affect this localization and LDL accumulation, but the details of the process remain unknown. It is thought that the LDL concentration will be affected by the filtration flow, and that the velocity of this flow will be affected by deformation of the arterial wall. Thus, a coupled model of a blood flow and a deformable arterial wall with filtration flow would be invaluable for simulation of the flow field and concentration field in sequence. However, this type of highly coupled interaction analysis has not yet been attempted. Therefore, we performed a coupled analysis of an artery with multiple bends in sequence. First, based on the theory of porous media, we modeled a deformable arterial wall using a porohyperelastic model (PHEM) that was able to express both the filtration flow and the viscoelastic behavior of the living tissue, and simulated a blood flow field in the arterial lumen, a filtration flow field and a displacement field in the arterial wall using a fluid-structure interaction (FSI) program code by the finite element method (FEM). Next, based on the obtained results, we further simulated LDL transport using a mass transfer analysis code by the FEM. We analyzed the PHEM in comparison with a rigid model. For the blood flow, stagnation was observed downward of the bends. The direction of the filtration flow was only from the lumen to the wall for the rigid model, while filtration flows from both the wall to the lumen and the lumen to the wall were observed for the PHEM. The LDL concentration was high at the lumen∕wall interface for both the PHEM and rigid model, and reached its maximum value at the stagnation area. For the PHEM, the maximum LDL concentration in the wall in the radial direction was observed at the position of 3% wall thickness from the lumen∕wall interface, while for the rigid model, it was observed just at the lumen∕wall interface. In addition, the peak LDL accumulation area of the PHEM moved about according to the pulsatile flow. These results demonstrate that the blood flow, arterial wall deformation, and filtration flow all affect the LDL concentration, and that LDL accumulation is due to stagnation and the presence of filtration flow. Thus, FSI analysis is indispensable.

Commentary by Dr. Valentin Fuster

### TECHNICAL PAPERS: Cell

J Biomech Eng. 2006;129(3):315-323. doi:10.1115/1.2720908.

The variations in mechanical properties of cells obtained from experimental and theoretical studies can be overcome only through the development of a sound mathematical framework correlating the derived mechanical property with the cellular structure. Such a formulation accounting for the inhomogeneity of the cytoplasm due to stress fibers and actin cortex is developed in this work. The proposed model is developed using the Mori-Tanaka method of homogenization by treating the cell as a fiber-reinforced composite medium satisfying the continuum hypothesis. The validation of the constitutive model using finite element analysis on atomic force microscopy (AFM) and magnetic twisting cytometry (MTC) has been carried out and is found to yield good correlation with reported experimental results. It is observed from the study that as the volume fraction of the stress fiber increases, the stiffness of the cell increases and it alters the force displacement behavior for the AFM and MTC experiments. Through this model, we have also been able to find the stress fiber as a likely cause of the differences in the derived mechanical property from the AFM and MTC experiments. The correlation of the mechanical behavior of the cell with the cell composition, as obtained through this study, is an important observation in cell mechanics.

Commentary by Dr. Valentin Fuster
J Biomech Eng. 2006;129(3):324-329. doi:10.1115/1.2720909.

The endothelial glycocalyx mediates interactions between the blood flow and the endothelium. This study aims to evaluate, quantitatively, effects of structural change of the glycocalyx on stress distribution and shear rate on endothelial cells. In the study, the endothelial glycocalyx is modeled as a surface layer of fiber matrix and when exposed to laminar shear flow, the matrix deforms. Fluid velocity and stress distribution inside the matrix and on cell membranes are studied based on a binary mixture theory. Parameters, such as the height and porosity of the matrix and the drag coefficient between fluid and matrix fibrils, are based on available data and estimation from experiments. Simple theoretical solutions are achieved for fluid velocity and stress distribution in the surface matrix. Degradation of the matrix, e.g., by enzyme digestion, is represented by reductions in the volume fraction of fibrils, height, and drag coefficient. From a force balance, total stress on endothelial surface remains constant regardless of structural alteration of the glycocalyx. However, the stress that is transmitted to endothelial cells by direct “pulling” of fiber branches of the glycocalyx is reduced significantly. Fluid shear rate at the cell membrane, on the other hand, increases. The study gives quantitative insight into the effect of the structural change of the glycocalyx on the shear rate and pulling stress on the endothelium. Results can be used to interpret experiments on effects of the glycocalyx in shear induced endothelial responses.

Commentary by Dr. Valentin Fuster

### TECHNICAL PAPERS: Bone/Orthopedic

J Biomech Eng. 2006;129(3):297-309. doi:10.1115/1.2720906.

The prediction of patient-specific proximal femur mechanical response to various load conditions is of major clinical importance in orthopaedics. This paper presents a novel, empirically validated high-order finite element method (FEM) for simulating the bone response to loads. A model of the bone geometry was constructed from a quantitative computerized tomography (QCT) scan using smooth surfaces for both the cortical and trabecular regions. Inhomogeneous isotropic elastic properties were assigned to the finite element model using distinct continuous spatial fields for each region. The Young’s modulus was represented as a continuous function computed by a least mean squares method. $p$-FEMs were used to bound the simulation numerical error and to quantify the modeling assumptions. We validated the FE results with in-vitro experiments on a fresh-frozen femur loaded by a quasi-static force of up to $1500N$ at four different angles. We measured the vertical displacement and strains at various locations and investigated the sensitivity of the simulation. Good agreement was found for the displacements, and a fair agreement found in the measured strain in some of the locations. The presented study is a first step toward a reliable $p$-FEM simulation of human femurs based on QCT data for clinical computer aided decision making.

Commentary by Dr. Valentin Fuster
J Biomech Eng. 2006;129(3):310-314. doi:10.1115/1.2720907.

Commentary by Dr. Valentin Fuster

### TECHNICAL PAPERS: Soft Tissue

J Biomech Eng. 2006;129(3):400-404. doi:10.1115/1.2721075.

Rotator cuff tears frequently occur and can lead to pain and decreased shoulder function. Repair of the torn tendon back to bone is often successful in relieving pain, but failure of the repair commonly occurs. Post-operative activity level is an important treatment component that has received minimal attention for the shoulder, but may have the potential to enhance tendon to bone healing. The objective of this study was to investigate the effect of short and long durations of various activity levels on the healing supraspinatus tendon to bone insertion site. Rotator cuff tears were surgically created in Sprague–Dawley rats by detaching the supraspinatus tendon from its insertion on the humerus and these tears were immediately repaired back to the insertion site. The post-operative activity level was controlled through shoulder immobilization (IM), cage activity (CA), or moderate exercise (EX) for durations of 4 or 16 weeks. The healing tissue was evaluated utilizing biomechanical testing and a quantitative polarized light microscopy method. We found that activity level had no effect on the elastic properties (stiffness, modulus) of the insertion site at four weeks post injury and repair, and a decreased activity level had a positive effect on these properties at 16 weeks $(IM>CA=EX)$. Furthermore, a decreased activity level had the greatest positive effect on these properties over time $(IM>CA=EX)$. The angular deviation of the collagen, a measure of disorganization, was decreased with a decrease in activity level at 4 weeks $(IM, but was similar between groups at 16 weeks $(IM=CA=EX)$. It appears from this study that decreasing the activity level by immobilizing the shoulder improves tendon to bone healing, which progresses by first increasing the organization of the collagen and then increasing the mechanical properties. Future studies in this area will investigate the effect of passive motion and remobilization on both tendon to bone healing and shoulder function.

Commentary by Dr. Valentin Fuster
J Biomech Eng. 2006;129(3):405-412. doi:10.1115/1.2720918.

Porous-permeable tissues have often been modeled using porous media theories such as the biphasic theory. This study examines the equivalence of the short-time biphasic and incompressible elastic responses for arbitrary deformations and constitutive relations from first principles. This equivalence is illustrated in problems of unconfined compression of a disk, and of articular contact under finite deformation, using two different constitutive relations for the solid matrix of cartilage, one of which accounts for the large disparity observed between the tensile and compressive moduli in this tissue. Demonstrating this equivalence under general conditions provides a rationale for using available finite element codes for incompressible elastic materials as a practical substitute for biphasic analyses, so long as only the short-time biphasic response is sought. In practice, an incompressible elastic analysis is representative of a biphasic analysis over the short-term response $δt⪡Δ2∕∥C4∥∥K∥$, where $Δ$ is a characteristic dimension, $C4$ is the elasticity tensor, and $K$ is the hydraulic permeability tensor of the solid matrix. Certain notes of caution are provided with regard to implementation issues, particularly when finite element formulations of incompressible elasticity employ an uncoupled strain energy function consisting of additive deviatoric and volumetric components.

Commentary by Dr. Valentin Fuster
J Biomech Eng. 2006;129(3):413-422. doi:10.1115/1.2720919.

Articular cartilage is a biological weight-bearing tissue covering the bony ends of articulating joints. Negatively charged proteoglycan (PG) in articular cartilage is one of the main factors that govern its compressive mechanical behavior and swelling phenomenon. PG is nonuniformly distributed throughout the depth direction, and its amount or distribution may change in the degenerated articular cartilage such as osteoarthritis. In this paper, we used a $50MHz$ ultrasound system to study the depth-dependent strain of articular cartilage under the osmotic loading induced by the decrease of the bathing saline concentration. The swelling-induced strains under the osmotic loading were used to determine the layered material properties of articular cartilage based on a triphasic model of the free-swelling. Fourteen cylindrical cartilage-bone samples prepared from fresh normal bovine patellae were tested in situ in this study. A layered triphasic model was proposed to describe the depth distribution of the swelling strain for the cartilage and to determine its aggregate modulus $Ha$ at two different layers, within which $Ha$ was assumed to be linearly dependent on the depth. The results showed that $Ha$ was $3.0±3.2$, $7.0±7.4$, $24.5±11.1MPa$ at the cartilage surface, layer interface, and deep region, respectively. They are significantly different $(p<0.01)$. The layer interface located at $70%±20%$ of the overall thickness from the uncalcified-calcified cartilage interface. Parametric analysis demonstrated that the depth-dependent distribution of the water fraction had a significant effect on the modeling results but not the fixed charge density. This study showed that high-frequency ultrasound measurement together with triphasic modeling is practical for quantifying the layered mechanical properties of articular cartilage nondestructively and has the potential for providing useful information for the detection of the early signs of osteoarthritis.

Commentary by Dr. Valentin Fuster
J Biomech Eng. 2006;129(3):423-429. doi:10.1115/1.2720920.

Cartilage is a charged hydrated fibrous tissue exhibiting a high degree of tension-compression nonlinearity (i.e., tissue anisotropy). The effect of tension-compression nonlinearity on solute transport has not been investigated in cartilaginous tissue under dynamic loading conditions. In this study, a new model was developed based on the mechano-electrochemical mixture model [Yao and Gu, 2007, J. Biomech. Model Mechanobiol., 6, pp. 63–72, Lai, 1991, J. Biomech. Eng., 113, pp. 245–258], and conewise linear elasticity model [Soltz and Ateshian, 2000, J. Biomech. Eng., 122, pp. 576–586;Curnier, 1995, J. Elasticity, 37, pp. 1–38]. The solute desorption in cartilage under unconfined dynamic compression was investigated numerically using this new model. Analyses and results demonstrated that a high degree of tissue tension-compression nonlinearity could enhance the transport of large solutes considerably in the cartilage sample under dynamic unconfined compression, whereas it had little effect on the transport of small solutes (at 5% dynamic strain level). The loading-induced convection is an important mechanism for enhancing the transport of large solutes in the cartilage sample with tension-compression nonlinearity. The dynamic compression also promoted diffusion of large solutes in both tissues with and without tension-compression nonlinearity. These findings provide a new insight into the mechanisms of solute transport in hydrated, fibrous soft tissues.

Commentary by Dr. Valentin Fuster
J Biomech Eng. 2006;129(3):430-440. doi:10.1115/1.2720924.

The atomic force microscope (AFM) has found wide applicability as a nanoindentation tool to measure local elastic properties of soft materials. An automated approach to the processing of AFM indentation data, namely, the extraction of Young’s modulus, is essential to realizing the high-throughput potential of the instrument as an elasticity probe for typical soft materials that exhibit inhomogeneity at microscopic scales. This paper focuses on Hertzian analysis techniques, which are applicable to linear elastic indentation. We compiled a series of synergistic strategies into an algorithm that overcomes many of the complications that have previously impeded efforts to automate the fitting of contact mechanics models to indentation data. AFM raster data sets containing up to 1024 individual force-displacement curves and macroscopic compression data were obtained from testing polyvinyl alcohol gels of known composition. Local elastic properties of tissue-engineered cartilage were also measured by the AFM. All AFM data sets were processed using customized software based on the algorithm, and the extracted values of Young’s modulus were compared to those obtained by macroscopic testing. Accuracy of the technique was verified by the good agreement between values of Young’s modulus obtained by AFM and by direct compression of the synthetic gels. Validation of robustness was achieved by successfully fitting the vastly different types of force curves generated from the indentation of tissue-engineered cartilage. For AFM indentation data that are amenable to Hertzian analysis, the method presented here minimizes subjectivity in preprocessing and allows for improved consistency and minimized user intervention. Automated, large-scale analysis of indentation data holds tremendous potential in bioengineering applications, such as high-resolution elasticity mapping of natural and artificial tissues.

Commentary by Dr. Valentin Fuster
J Biomech Eng. 2006;129(3):441-449. doi:10.1115/1.2721076.

Early in development, the heart is a single muscle-wrapped tube without formed valves. Yet survival of the embryo depends on the ability of this tube to pump blood at steadily increasing rates and pressures. Developmental biologists historically have speculated that the heart tube pumps via a peristaltic mechanism, with a wave of contraction propagating from the inflow to the outflow end. Physiological measurements, however, have shown that the flow becomes pulsatile in character quite early in development, before the valves form. Here, we use a computational model for flow though the embryonic heart to explore the pumping mechanism. Results from the model show that endocardial cushions, which are valve primordia arising near the ends of the tube, induce a transition from peristaltic to pulsatile flow. Comparison of numerical results with published experimental data shows reasonably good agreement for various pressure and flow parameters. This study illustrates the interrelationship between form and function in the early embryonic heart.

Commentary by Dr. Valentin Fuster
J Biomech Eng. 2006;129(3):450-456. doi:10.1115/1.2720928.

The mechanical response of soft tissue is commonly characterized from unconfined uniaxial compression experiments on cylindrical samples. However, friction between the sample and the compression platens is inevitable and hard to quantify. One alternative is to adhere the sample to the platens, which leads to a known no-slip boundary condition, but the resulting nonuniform state of stress in the sample makes it difficult to determine its material parameters. This paper presents an approach to extract the nonlinear material properties of soft tissue (such as liver) directly from no-slip experiments using a set of computationally determined correction factors. We assume that liver tissue is an isotropic, incompressible hyperelastic material characterized by the exponential form of strain energy function. The proposed approach is applied to data from experiments on bovine liver tissue. Results show that the apparent material properties, i.e., those determined from no-slip experiments ignoring the no-slip conditions, can differ from the true material properties by as much as 50% for the exponential material model. The proposed correction approach allows one to determine the true material parameters directly from no-slip experiments and can be easily extended to other forms of hyperelastic material models.

Commentary by Dr. Valentin Fuster
J Biomech Eng. 2006;129(3):457-471. doi:10.1115/1.2737056.

A three-dimensional (3D) contact finite element formulation has been developed for biological soft tissue-to-tissue contact analysis. The linear biphasic theory of Mow, Holmes, and Lai (1984, J. Biomech., 17(5), pp. 377–394) based on continuum mixture theory, is adopted to describe the hydrated soft tissue as a continuum of solid and fluid phases. Four contact continuity conditions derived for biphasic mixtures by Hou (1989, ASME J. Biomech. Eng., 111(1), pp. 78–87) are introduced on the assumed contact surface, and a weighted residual method has been used to derive a mixed velocity-pressure finite element contact formulation. The Lagrange multiplier method is used to enforce two of the four contact continuity conditions, while the other two conditions are introduced directly into the weighted residual statement. Alternate formulations are possible, which differ in the choice of continuity conditions that are enforced with Lagrange multipliers. Primary attention is focused on a formulation that enforces the normal solid traction and relative fluid flow continuity conditions on the contact surface using Lagrange multipliers. An alternate approach, in which the multipliers enforce normal solid traction and pressure continuity conditions, is also discussed. 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. The resulting first-order system of equations is solved in time using the generalized finite difference scheme. The formulation is validated by a series of increasingly complex canonical problems, including the confined and unconfined compression, 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 is analyzed using an idealized 3D geometry. In the joint, both glenoid and humeral head cartilage experience maximum tensile and compressive stresses are at the cartilage-bone interface, away from the center of the contact area.

Commentary by Dr. Valentin Fuster