Multibody dynamic simulation of knee contact mechanics

Abstract

Multibody dynamic musculoskeletal models capable of predicting muscle forces and joint contact pressures simultaneously would be valuable for studying clinical issues related to knee joint degeneration and restoration. Current three-dimensional multibody knee models are either quasi-static with deformable contact or dynamic with rigid contact. This study proposes a computationally efficient methodology for combining multibody dynamic simulation methods with a deformable contact knee model. The methodology requires preparation of the articular surface geometry, development of efficient methods to calculate distances between contact surfaces, implementation of an efficient contact solver that accounts for the unique characteristics of human joints, and specification of an application programming interface for integration with any multibody dynamic simulation environment. The current implementation accommodates natural or artificial tibiofemoral joint models, small or large strain contact models, and linear or nonlinear material models. Applications are presented for static analysis (via dynamic simulation) of a natural knee model created from MRI and CT data and dynamic simulation of an artificial knee model produced from manufacturer's CAD data. Small and large strain natural knee static analyses required 1 min of CPU time and predicted similar contact conditions except for peak pressure, which was higher for the large strain model. Linear and nonlinear artificial knee dynamic simulations required 10 min of CPU time and predicted similar contact force and torque but different contact pressures, which were lower for the nonlinear model due to increased contact area. This methodology provides an important step toward the realization of dynamic musculoskeletal models that can predict in vivo knee joint motion and loading simultaneously.

Fig. 1

Contact surface preparation demonstrated for an artificial knee model. Starting from manufacturer’s CAD data (a), the trimmed NURBS contact surfaces (b) are removed from the model. Commercial software (see text) is used to fit untrimmed NURBS patches to the contact surfaces and then merge them into a single patch (c). Finally, other commercial software (see text) is used to minimize the number of parametric B-spline curves representing the surface (d) while maintaining a specified level of surface accuracy.

Fig. 2

Classification of surface types analyzed during contact calculations. Element surfaces created on the tibia are used to construct a planar element grid on the medial and lateral sides. Contact surfaces on the tibia and femur represent the contacting articular geometry. Back surfaces only exist for natural knees and represent the subchondral bone surfaces on the tibia and femur needed to calculate local cartilage thickness. Projection of element boundaries from the element surfaces onto the tibial contact surfaces produces contact elements with unique normal directions, coordinate systems, thicknesses, and areas.

Fig. 3

Distance evaluation methods for contact calculations. (a) The tangent plane commonly used for nonconformal contact situations requires two sets of distance calculations—one from the tangent plane to each moving body (tibia and femur). (b) The implicit midsurface proposed in our approach works for conformal and nonconformal contact situations and requires half the number of distance calculations. For each contact element on the fixed body (tibia), the distance to the corresponding moving body (femur) contact surface is calculated directly and then corrected. (c) The corrected distance vector is calculated from trigonometry based on knowledge of the minimum distance and ray firing directions. Either minimum distance or ray firing can be used to calculate distances. The minimum distance vector (AB) is perpendicular to the point found on the moving body, while the ray firing vector (AC) is perpendicular to the normal on the fixed body. If ray firing is used, the minimum distance direction can be well approximated from the normal of the point found on the moving body. The distance along the midsurface normal, defined from the bisector between the minimum distance and ray firing directions, is larger than the minimum distance result and smaller than the ray firing result. If minimum distance is used, the distance is grown slightly by calculating the distance (AD) that produces the minimum distance (AB) when projected onto the minimum distance direction. If ray firing is used, the distance is shrunk slightly by projecting it onto the minimum distance direction (AE) and then growing the distance (AF) using the same approach as for minimum distance. In practice, corrected distances from the two approaches are nearly identical given typical element sizes.

Fig. 4

Static analysis results (via dynamic simulation) for the natural knee contact model using the large strain formulation. (a) Animation of dynamic simulation progression from initial pose to static pose. (b) Contact pressure contours on the tibial condyles at the final static pose. (c) Motions predicted by the dynamic simulation. Y translation is superior–inferior, Z translation is medial–lateral, and X rotation is varus–valgus.

Fig. 5

Dynamic simulation results for the artificial knee contact model using the linear material formulation. (a) Animation of dynamic simulation progression for one cycle of gait. (b) Contact pressure contours on the tibial insert at five locations in the gait cycle. (c) Motions predicted by the dynamic simulation (solid lines) compared to accuracy envelopes for the fluoroscopic measurements (gray bands) [10]. Y translation is superior–inferior, Z translation is medial–lateral, and X rotation is varus–valgus.

Fig. 6

Sensitivity of contact predictions to the number of element divisions in the x and z directions for the artificial knee model. (a) Contact force. (b) Contact torque. (c) Peak pressure. (d) Average pressure. The sensitivity can be different in the two directions. Contact force and torque are much less sensitive to the element grid than are peak and average contact pressure. The high average pressure sensitivity is due to contact area sensitivity (not shown).