On foundations of discrete element analysis of contact in diarthrodial joints

Information about the stress distribution on contact surfaces of adjacent bones is indispensable for analysis of arthritis, bone fracture and remodeling. Numerical solution of the contact problem based on the classical approaches of solid mechanics is sophisticated and time-consuming. However, the solution can be essentially simplified on the following physical grounds. The bone contact surfaces are covered with a layer of articular cartilage, which is a soft tissue as compared to the hard bone. The latter allows ignoring the bone compliance in analysis of the contact problem, i.e. rigid bones are considered to interact through a compliant cartilage. Moreover, cartilage shear stresses and strains can be ignored because of the negligible friction between contacting cartilage layers. Thus, the cartilage can be approximated by a set of unilateral compressive springs normal to the bone surface. The forces in the springs can be computed from the equilibrium equations iteratively accounting for the changing contact area. This is the essence of the discrete element analysis (DEA). Despite the success in applications of DEA to various bone contact problems, its classical formulation required experimental validation because the springs approximating the cartilage were assumed linear while the real articular cartilage exhibited non-linear mechanical response in reported tests. Recent experimental results of Ateshian and his co-workers allow for revisiting the classical DEA formulation and establishing the limits of its applicability. In the present work, it is shown that the linear spring model is remarkably valid within a wide range of large deformations of the cartilage. It is also shown how to extend the classical DEA to the case of strong nonlinearity if necessary.     

An augmented Lagrangian finite element formulation for 3D contact of biphasic tissues

Biphasic contact analysis is essential to obtain a complete understanding of soft tissue biomechanics, and the importance of physiological structure on the joint biomechanics has long been recognised; however, up to date, there are no successful developments of biphasic finite element contact analysis for three-dimensional (3D) geometries of physiological joints. The aim of this study was to develop a finite element formulation for biphasic contact of 3D physiological joints. The augmented Lagrangian method was used to enforce the continuity of contact traction and fluid pressure across the contact interface. The biphasic contact method was implemented in the commercial software COMSOL Multiphysics 4.2^® (COMSOL, Inc., Burlington, MA). The accuracy of the implementation was verified using 3D biphasic contact problems, including indentation with a flat-ended indenter and contact of glenohumeral cartilage layers. The ability of the method to model multibody biphasic contact of physiological joints was proved by a 3D knee model. The 3D biphasic finite element contact method developed in this study can be used to study the biphasic behaviours of the physiological joints.     

Glenohumeral joint cartilage contact in the healthy adult during scapular plane elevation depression with external humeral rotation

The shoulder (glenohumeral) joint has the greatest range of motion of all human joints; as a result, it is particularly vulnerable to dislocation and injury. The ability to non-invasively quantify in-vivo articular cartilage contact patterns of joints has been and remains a difficult biomechanics problem. As a result, little is known about normal in-vivo glenohumeral joint contact patterns or the consequences that surgery has on altering them. In addition, the effect of quantifying glenohumeral joint contact patterns by means of proximity mapping, both with and without cartilage data, is unknown. Therefore, the objectives of this study are to (1) describe a technique for quantifying in-vivo glenohumeral joint contact patterns during dynamic shoulder motion, (2) quantify normal glenohumeral joint contact patterns in the young healthy adult during scapular plane elevation depression with external humeral rotation, and (3) compare glenohumeral joint contact patterns determined both with and without articular cartilage data. Our results show that the inclusion of articular cartilage data when quantifying in-vivo glenohumeral joint contact patterns has significant effects on the anterior–posterior contact centroid location, the superior–inferior contact centroid range of travel, and the total contact path length. As a result, our technique offers an advantage over glenohumeral joint contact pattern measurement techniques that neglect articular cartilage data. Likewise, this technique may be more sensitive than traditional 6-Degree-of-Freedom (6-DOF) joint kinematics for the assessment of overall glenohumeral joint health. Lastly, for the shoulder motion tested, we found that glenohumeral joint contact was located on the anterior–inferior glenoid surface.     

A geometric approach to study the contact mechanisms in the patellofemoral joint of normal versus patellofemoral pain syndrome subjects

The biomechanics of the patellofemoral (PF) joint is complex in nature, and the aetiology of such manifestations of PF instability as patellofemoral pain syndrome (PFPS) is still unclear. At this point, the particular factors affecting PFPS have not yet been determined. This study has two objectives: (1) The first is to develop an alternative geometric method using a three-dimensional (3D) registration technique and linear mapping to investigate the PF joint contact stress using an indirect measure: the depth of virtual penetration (PD) of the patellar cartilage surface into the femoral cartilage surface. (2) The second is to develop 3D PF joint models using the finite element analysis (FEA) to quantify in vivo cartilage contact stress and to compare the peak contact stress location obtained from the FE models with the location of the maximum PD. Magnetic resonance images of healthy and PFPS subjects at knee flexion angles of 15°, 30° and 45° during isometric loading have been used to develop the geometric models. The results obtained from both approaches demonstrated that the subjects with PFPS show higher PD and contact stresses than the normal subjects. Maximum stress and PD increase with flexion angle, and occur on the lateral side in healthy and on the medial side in PFPS subjects. It has been concluded that the alternative geometric method is reliable in addition to being computationally efficient compared with FEA, and has the potential to assess the mechanics of PFPS with an accuracy similar to the FEA.     

An analytical model of joint contact

The stress distribution in the region of contact between a layered elastic sphere and a layered elastic cavity is determined using an analytical model to stimulate contact of articulating joints. The purpose is to use the solution to analyze the effects of cartilage thickness and stiffness, bone stiffness and joint curvature on the resulting stress field, and investigate the possibility of cracking of the material due to tensile and shear stresses. Vertical cracking of cartilage as well as horizontal splitting at the cartilage-calcified cartilage interface has been observed in osteoarthritic joints. The current results indicate that for a given system (material properties mu and nu constant), the stress distribution is a function of the ratio of contact radius to layer thickness (a/h), and while tensile stresses are seen to occur only when a/h is small, tensile strain is observed for all a/h values. Significant shear stresses are observed at the cartilage-bone interface. Softening of cartilage results in an increase in a/h, and a decrease in maximum normal stress. Cartilage thinning increases a/h and the maximum contact stress, while thickening has the opposite effect. A reduction in the indenting radius reduces a/h and increases the maximum normal stress. Bone softening is seen to have negligible effect on the resulting contact parameters and stress distribution.     

A micromechanical model to predict damage and failure in biological tissues. Application to the ligament-to-bone attachment in the human knee joint

Computational models are developed in injury biomechanics to assess lesions in biological tissues based on mechanical measurements. The linear mechanics of fracture theory (LMFT) is a common approach to establish injuries based on thresholds (such as force or strain thresholds) which are straightforward to implement and computationally efficient. However, LMFT does not apply to non-linear heterogeneous materials and does not have the ability to predict failure onset. This paper proposes the cohesive zone model theory (CZMT) as an alternative. CZMT focuses on the development of behaviour laws for crack initiation and propagation at an interface that apply within a fibrous material or at the interface between materials. With the view of evaluating CZMT for biological tissues, the model developed by Raous et al. [1999. A consistent model coupling adhesion, friction and unilateral contact. Comput. Methods Appl. Mech. Eng., 177, 383–399] was applied to the ligament-to-bone interface in the human knee joint. This model accounts for adhesion, friction and damage at the interface and provides a smooth transition from total adhesion to complete failure through the intensity of adhesion variable. A 2D finite element model was developed to mimic previous experiments, and the model parameters were determined using a dichotomy method. The model showed good results by its ability to predict damage. The extension to a 3D geometry, with an inverse problem approach, is, however, required to better estimate the model parameters values. Although it is computationally costly, CZMT supplements the improvements achieved in microimaging techniques to support the development of micro/macro approaches in biomechanical modelling.     

Hand biomechanics during simulated stone tool use

Human radial digits have derived features compared with apes, with long robust thumbs, relatively larger joint surfaces, and hypertrophic thenar muscles. Here we test the hypothesis that these features evolved in the context of making and using stone tools, specifically for producing large gripping forces and for countering large joint contact stresses. We used portable force plates simulating early stone tools to: 1) document and compare the magnitude of external/internal forces and joint stresses in the radial digits during hardhammer percussion and flake use, and 2) examine how variation in digit morphology affects muscle and joint mechanics during stone tool use. Force and kinematic data were collected from a sample representing normal variation in digit morphology (n = 25). The effects of digit size/shape on digit biomechanics were evaluated using partial correlations, controlling for tool reaction forces and impact velocities. Results show that individuals with longer digits require relatively less muscle force to stabilize digital joints, and are exposed to relatively lower joint contact stresses during stone tool use, due in part to an increase in the robusticity of metacarpals and phalanges in humans relative to chimpanzees. These analyses further suggest that Pan- or australopith-like pollical anatomy presents serious performance challenges to habitual tool use. Our data support the hypothesis that evolutionary increases in thumb length, robusticity, and thenar muscle mass enabled Homo to produce more force and to tolerate higher joint stresses during tool use.     

Scaling of limb joint surface areas in anthropoid primates and other mammals

Recently it was reported that limb joint surface areas scale positively allometrically with body weight in anthropoid primates. This was attributed to the biomechanics of weight bearing: larger animals must require relatively larger joint surfaces to withstand relatively greater weight-related stresses on the joints. Our data on humeral and femoral geometry and joint surface areas in 73 species belonging to six mammalian orders (including primates) demonstrate that positive allometry of joint surface areas is not a general phenomenon for mammals and cannot have its basis among Anthropoidea in the biomechanics of weight bearing. We argue that, to the extent that positive allometric scaling of joint surface areas occurs in anthropoid primates, it is an artifact of differences in positional behaviour among distinct taxonomic groups that also happen to differ widely in body weight. Furthermore, we argue that, among mammals ranging in body size from sportive lemurs to giant pandas, functionally similar groups tend to exhibit: (1) linear dimensions (especially diameters and shaft circumferences) that scale in direct proportion to each other; and (2) joint surface areas that scale in direct proportion to the squares of linear dimensions and to the 2/3 power of body mass. In other words, limb bones of functionally similar animals fit the theoretical model of geometric similarity (or skeletal isometry). Differences in relative sizes of joint surface areas are related to differences in force transmission and movement potential among functionally distinct groups of animals.     

Biomechanical model calculation of muscle contraction forces: a double linear programming method

This paper presents a novel scheme for the use of linear programming to calculate muscle contraction forces in models describing musculoskeletal system biomechanics. Models of this kind are frequently found in the biomechanics literature. In most cases they involve muscle contraction force calculations that are statically indeterminate, and hence use optimization techniques to make those calculations. We present a linear programming optimization technique that solves a two-objective problem with two sequential linear programs. We use the technique here to minimize muscle intensity and joint compression force, since those are commonly used objectives. The two linear program model has the advantages of low computation cost, ready implementation on a micro-computer, and stable solutions. We show how to solve the model analytically in simple cases. We also discuss the use of the dual problem of linear programming to gain understanding of the solution it provides.     

An augmented lagrangian method for sliding contact of soft tissue

Despite the importance of sliding contact in diarthrodial joints, only a limited number of studies have addressed this type of problem, with the result that the mechanical behavior of articular cartilage in daily life remains poorly understood. In this paper, a finite element formulation is developed for the sliding contact of biphasic soft tissues. The augmented Lagrangian method is used to enforce the continuity of contact traction and fluid pressure across the contact interface. The resulting method is implemented in the commercial software COMSOL Multiphysics. The accuracy of the new implementation is verified using an example problem of sliding contact between a rigid, impermeable indenter and a cartilage layer for which analytical solutions have been obtained. The new implementation's capability to handle a complex loading regime is verified by modeling plowing tests of the temporomandibular joint (TMJ) disc.     

A contact algorithm for density-based load estimation

An algorithm, which includes contact interactions within a joint, has been developed to estimate the dominant loading patterns in joints based on the density distribution of bone. The algorithm is applied to the proximal femur of a chimpanzee, gorilla and grizzly bear and is compared to the results obtained in a companion paper that uses a non-contact (linear) version of the density-based load estimation method. Results from the contact algorithm are consistent with those from the linear method. While the contact algorithm is substantially more complex than the linear method, it has some added benefits. First, since contact between the two interacting surfaces is incorporated into the load estimation method, the pressure distributions selected by the method are more likely indicative of those found in vivo. Thus, the pressure distributions predicted by the algorithm are more consistent with the in vivo loads that were responsible for producing the given distribution of bone density. Additionally, the relative positions of the interacting bones are known for each pressure distribution selected by the algorithm. This should allow the pressure distributions to be related to specific types of activities. The ultimate goal is to develop a technique that can predict dominant joint loading patterns and relate these loading patterns to specific types of locomotion and/or activities.     

Development and validation of a kinematically-driven discrete element model of the patellofemoral joint

Quantifying the complex loads at the patellofemoral joint (PFJ) is vital to understanding the development of PFJ pain and osteoarthritis. Discrete element analysis (DEA) is a computationally efficient method to estimate cartilage contact stresses with potential application at the PFJ to better understand PFJ mechanics. The current study validated a DEA modeling framework driven by PFJ kinematics to predict experimentally-measured PFJ contact stress distributions. Two cadaveric knee specimens underwent quadriceps muscle [215 N] and joint compression [350 N] forces at ten discrete knee positions representing PFJ positions during early gait while measured PFJ kinematics were used to drive specimen-specific DEA models. DEA-computed contact stress and area were compared to experimentally-measured data. There was good agreement between computed and measured mean and peak stress across the specimens and positions (r = 0.63–0.85). DEA-computed mean stress was within an average of 12% (range: 1–47%) of the experimentally-measured mean stress while DEA-computed peak stress was within an average of 22% (range: 1–40%). Stress magnitudes were within the ranges measured (0.17–1.26 MPa computationally vs 0.12–1.13 MPa experimentally). DEA-computed areas overestimated measured areas (average error = 60%; range: 4–117%) with magnitudes ranging from 139 to 307 mm² computationally vs 74–194 mm² experimentally. DEA estimates of the ratio of lateral to medial patellofemoral stress distribution predicted the experimental data well (mean error = 15%) with minimal measurement bias. These results indicate that kinematically-driven DEA models can provide good estimates of relative changes in PFJ contact stress.     

Finite element analysis of shear stresses at the implant-bone interface of an acetabular press-fit cup during impingement.

After total hip replacement (THR) impingement of the implant components causes shear stresses at the acetabular implant-bone interface. In the current study the finite element method (FEM) was applied to analyse the shear stresses at a fully bonded implant-bone interface assuming total ingrowth of the cup. The FE model of a press-fit acetabular component and the proximal part of the femoral component incorporates non-linear material and large sliding contact. The model was loaded with a superior-medial joint load of 435 N simulating a two-legged stance. Starting at initial impingement, the femoral component was medially rotated by 20 degrees . The peak tilting shear stress of -2.6 MPa at the impingement site takes effect towards the pole of the cup. The torsional shear stress at the impingement site is zero. On each side of the impingement site, there are extrema of torsional shear stress reaching -1.8 and 1.8 MPa, respectively. The global peak shear stress during impingement may indicate a possible starting point for cup loosening. The pattern of the torsional shear stresses suggests that besides the symmetric lever-out, an additional asymmetrical tilting of the cup occurs that can be explained by the orientation of the applied joint load.     

Effect of meniscus replacement fixation technique on restoration of knee contact mechanics and stability

The menisci are important biomechanical components of the knee. We developed and validated a finite element model of meniscal replacement to assess the effect of surgical fixation technique on contact behavior and knee stability. The geometry of femoral and tibial articular cartilage and menisci was segmented from magnetic resonance images of a normal cadaver knee using MIMICS (Materialise, Leuven, Belgium). A finite element mesh was generated using HyperWorks (Altair Inc, Santa Ana, CA). A finite element solver (Abaqus v6.9, Simulia, Providence, RI) was used to compute contact area and stresses under axial loading and to assess stability (reaction force generated during anteroposterior translation of the femur). The natural and surgical attachments of the meniscal horns and peripheral rim were simulated using springs. After total meniscectomy, femoral contact area decreased by 26% with a concomitant increase in average contact stresses (36%) and peak contact stresses (33%). Replacing the meniscus without suturing the horns did little to restore femoral contact area. Suturing the horns increased contact area and reduced peak contact stresses. Increasing suture stiffness correlated with increased meniscal contact stresses as a greater proportion of tibiofemoral load was transferred to the meniscus. A small incremental benefit was seen of simulated bone plug fixation over the suture construct with the highest stiffness (50 N/mm). Suturing the rim did little to change contact conditions. The nominal anteroposterior stiffness reduced by 3.1 N/mm after meniscectomy. In contrast to contact area and stress, stiffness of the horn fixation sutures had a smaller effect on anteroposterior stability. On the other hand suturing the rim of the meniscus affected anteroposterior stability to a much larger degree. This model emphasizes the importance of the meniscus in knee biomechanics. Appropriate meniscal replacement fixation techniques are likely to be critical to the clinical success of meniscal replacement. While contact conditions are mainly sensitive to meniscus horn fixation, the stability of the knee under anteroposterior shear loads appeared to be more sensitive to meniscal rim fixation. This model may also be useful in predicting the effect of biomaterial mechanical properties and meniscal replacement shape on knee contact conditions.     

Computational hip joint simulator for wear and heat generation

This paper presents a computational simulator for the hip to compute the wear and heat generation on artificial joints. The friction produced on artificial hip joints originates wear rates that can lead to failure of the implant. Furthermore, the frictional heating can increase the wear. The developed computational model calculates the wear in the joint and the temperature in the surrounding zone, allowing the use of different combinations of joint materials, daily activities and different individuals. The pressure distribution on the joint bearing surfaces is obtained with the solution of a contact model. The heat generation by friction and the volumetric wear is computed from the pressure distribution and the sliding distance. The temperature is obtained from the solution of a transient heat conduction problem that includes the time-dependent heat generated by friction. The contact and heat conduction problems are solved numerically with the Finite Element Method. The developed computational model performs a full simulation of the acetabular bearing surface behaviour, which is useful for acetabular cup design and material selection. The results obtained by the present model agree with experimental and clinical data, as well as other numerical studies.     

Effects of inserting a pressensor film into articular joints on the actual contact mechanics.

Fuji film has been widely used in studies aimed at obtaining the contact mechanics of articular joints. Once sealed for practical use in biological joints, Fuji Pressensor film has a total effective thickness of 0.30 mm, which is comparable to the cartilage thickness in the joints of many small animals. The average effective elastic modulus of Fuji film is approximately 100 MPa in compression, which is larger by a factor of 100-300 compared to that of normal articular cartilage. Therefore, inserting a Pressensor film into an articular joint will change the contact mechanics of the joint. The measurement precision of the Pressensor film has been determined systematically; however, the changes in contact mechanics associated with inserting the film into joints have not been investigated. This study was aimed at quantifying the changes in the contact mechanics associated with inserting sealed Fuji Pressensor film into joints. Spherical and cylindrical articular joint contact mechanics with and without Pressensor film and for varying degrees of surface congruency were analyzed and compared by using finite element models. The Pressensor film was taken as linearly elastic and the cartilage was assumed to be biphasic, composed of a linear elastic solid phase and an inviscid fluid phase. The present analyses showed that measurements of the joint contact pressures with Fuji Pressensor film will change the maximum true contact pressures by 10-26 percent depending on the loading, geometry of the joints, and the mechanical properties of cartilage. Considering this effect plus the measurement precision of the film (approximately 10 percent), the measured joint contact pressures in a joint may contain errors as large as 14-28 percent.     

Image based weighted center of proximity versus directly measured knee contact location during simulated gait

To understand the mechanical consequences of knee injury requires a detailed analysis of the effect of that injury on joint contact mechanics during activities of daily living. Three-dimensional (3D) knee joint geometric models have been combined with knee joint kinematics to dynamically estimate the location of joint contact during physiological activities—using a weighted center of proximity (WCoP) method. However, the relationship between the estimated WCoP and the actual location of contact has not been defined. The objective of this study was to assess the relationship between knee joint contact location as estimated using the image-based WCoP method, and a directly measured weighted center of contact (WCoC) method during simulated walking. To achieve this goal, we created knee specific models of six human cadaveric knees from magnetic resonance imaging. All knees were then subjected to physiological loads on a knee simulator intended to mimic gait. Knee joint motion was captured using a motion capture system. Knee joint contact stresses were synchronously recorded using a thin electronic sensor throughout gait, and used to compute WCoC for the medial and lateral plateaus of each knee. WCoP was calculated by combining knee kinematics with the MRI-based knee specific model. Both metrics were compared throughout gait using linear regression. The anteroposterior (AP) location of WCoP was significantly correlated with that of WCoC on both tibial plateaus in all specimens (p<0.01, 95% confidence interval of Pearson?s coefficient r>0), but the correlation was not significant in the mediolateral (ML) direction for 4/6 knees (p>0.05). Our study demonstrates that while the location of joint contact obtained from 3D knee joint contact model, using the WCoP method, is significantly correlated with the location of actual contact stresses in the AP direction, that relationship is less certain in the ML direction.     

Stresses on joints of Opuntia laevis: forces necessary for mechanical failure and role of lignified xylem cells as resistance components

Species of Opuntia exhibit a wide range of morphologies. Understanding these morphologies may require knowledge of the mechanical stresses on joints of stem segments and as well as the internal components in joints that withstand joint failure (separation of the terminal cladode from the sub-terminal cladode after weights were applied perpendicularly to the long axis). Results of stress testing terminal cladodes of Opuntia laevis provided the following conclusions: (1) amounts of applied stress for joint failure were not related to the amounts of stress on joints before stress testing; (2) breaking strength (failure stress) was accurately determined for joints from linear plots of M (bending moment) versus I/c (section modulus) [breaking stress for O. laevis was 2.77 kPa]; (3) bending moments at failure were twice as high for tensile portions than for compressive portions of joints; and (4) bending moments at failure were positively correlated with amounts of lignified xylem cells in joints [for each mm² of lignified xylem cells in joints there was an increase of 0.06 N m of bending moment]. These data support the overall hypothesis that bending stresses are the main stresses at joints of Opuntia laevis and that lignified xylem cells are the main components that resist joint failure. Moreover, since tensile portions have more lignified xylem cells than other stem portions, tensile portions can resist more applied stress.     

Finite element algorithm for frictionless contact of porous permeable media under finite deformation and sliding

This study formulates and implements a finite element contact algorithm for solid-fluid (biphasic) mixtures, accommodating both finite deformation and sliding. The finite element source code is made available to the general public. The algorithm uses a penalty method regularized with an augmented Lagrangian method to enforce the continuity of contact traction and normal component of fluid flux across the contact interface. The formulation addresses the need to automatically enforce free-draining conditions outside of the contact interface. The accuracy of the implementation is verified using contact problems, for which exact solutions are obtained by alternative analyses. Illustrations are also provided that demonstrate large deformations and sliding under configurations relevant to biomechanical applications such as articular contact. This study addresses an important computational need in the biomechanics of porous-permeable soft tissues. Placing the source code in the public domain provides a useful resource to the biomechanics community.