0021-9290,W S3.C0+.CKI Pcrgamon PIma plc
J. Biowchanics Vol. 23, No. 10, PP. 1013-1020, 1990. Printed in Great Britain
PRESSURE DISTRIBUTION ON ARTICULAR SURFACES: APPLICATION TO JOINT STABILITY EVALUATION K. N. AN, S. HIMENO, H. TSUMURA, T. KAWAI and E. Y. S.
CHAO
Biomechanics I ‘aboratory, Department of Orthopedics, Mayo Clinic/Mayo Foundation, Rochester, MN 55905, U.S.A. Abutract-A generalized method, based on the rigid body spring model, is developed to determine the twodimensional joint contact pressure of any shape of articular surface and loading condition. The formulation of the model and procedures for the determination of parameters used in the model are presented. The model was used first to analyze the pressure distribution of an ideal hinge joint for comparison with reported results in the literature, and with those using the conventional finite element method. Implications of the results based on the calculated ‘virtual displacement’ to the joint stability are hypothesized and discussed. The problem of joint pressure distribution of the elbow joint is then analyzed, along with the discussion of agonistic and antagonistic muscle activities.
as a result of the analysis itself, and this distribution was treated as external load on each of the joint Precise knowledge of pressure distribution in joints is components. not only a prerequisite for the determination of stresses Brown and DiGioia (1984) discussed a contactin articular cartilage, but might also be the determincoupled finite element analysis of the hip joint which ing factor of the joint stability and thus influences the could calculate the contact force distribution and distribution of muscle force under synergistic and internal stress simultaneously, given the incremental antagonistic conditions. displacements of the acetabular component as the The concept of joint pressure distribution and its boundary condition and iterative calculation to deterapplication was discussed extensively by Pauwels mine the surface members in contact. This was, un(1980) by considering a simple, hinged elbow joint. doubtedly, a great improvement of the situation; Mathematical analyses of force transmission were however, there were still several problems. A long carried out by Greenwald and O’Connor (1971) computing time was required to accomplish this anathrough the hip joint and by Wynarsky and Greenlysis, even though a rather coarse zoning was adopted wald (1983) through the ankle joint. In these two to lessen the amount of calculation. Displacements mathematical analyses, it was assumed that deformawere given rather than loads as boundary conditions. tion took place only in the cartilage layer, treated as a This feature may face some difficulties when the simple elastic material, and that displacement of a method is applied to multibody contact problems, point was independent of the normal stress applied at such as in the wrist or foot joints, because the a priori any other point. They concluded that the distribution displacements of intermediate components cannot be of pressure depended critically on variations in cartil- given. age thickness, and that it also depended on the initial In certain cases where the actual cartilage stress incongruity of the joint at low loads. In their analysis, distribution is of little direct interest, as opposed to the ideal joint geometry was assumed. However, for joints contact stress on the articulating surface, a simpler of irregular configuration, these analytic methods are approach could be performed. In this study, a gendifficult to apply and additional assumptions are eralized method to solve such problems of joint presrequired. sure distribution was developed. The derivation of this Finite element analysis has thus been adopted. It is model was based on the minimum strain energy a well established method. Geometrical complexities principle as used in most of the formulation of finite are easily solved by using fine mesh division. Bound- element codes. An analytic and numeric model was ary conditions are easy to modify, and most of the formulated by using the rigid body spring method of nonlinearities can be handled. Flexibility in modeling Kawai (1977) and Kawai and Takeuchi (1981). The is excellent in this method. As for articulating contact, results were first verified by comparison with those however, most studies assumed the contact force analytic results in the literature for the special case of a distribution on the joint surface a priori rather than hinge joint (Pauwels, 1980). The effects of size of the articular surface and the location of the resultant joint force on the pressure distribution were then examined. Receivedin final form 21 March 1990. The implication of the results on the joint stability as Correspondence to: Kai-Nan An, Ph.D., Orthopedic Bio- well as their influence on the activities of ago&tic and mechanics Lab., Mayo Clinic/Mayo Foundation, 200 First Street, Southwest, Rochester, MN 55905, U.S.A. Tel.: (507) antagonistic muscles are also discussed by considering the elbow joint. 284-2589. 1013
K. N.
1014
AN et al.
METHODS
To establish the analytic and numeric model for joint pressure determination, the shapes of articular surfaces were either mathematically or numerically described. The rigid body spring model was then formulated by considering the rigid body to be in equilibrium with external loading. Reaction forces between adjacent bodies were produced by the spring system distributed over the possible contact surface between the two adjacent bodies. Displacements of the simulated compressive spring, normal and tangential to the articular surface, were described as functions of the displacements of the centroid and rotation of the associated rigid bodies. The capsule and ligaments Fig. 1. Definition of parameters used for joint pressure diswere also modeled with tensile springs, whose defor- tribution study of an ideal circular arc. The available articumation could also be described based on the displace- lating surfaces are expressed by the subtended angle of the ment of the centroids, and rotation of the rigid bodies arc $. The externally applied force, R, is located at angle 0 from the center of the circular arc. as well. Strain energy due to relative displacement of the spring system was formulated using a quadratic function of the displacement vectors of the centroid of Y adjacent rigid bodies. By applying the principle of minimum potential energy, the stiffness equation was obtained by differentiating strain energy with respect to displacement of centroids of adjacent rigid bodies. A system of equations was obtained for solving Y.&T* centroid displacement of rigid bodies. Displacement YP and force in the spring system which describe the pressure distribution on articular surfaces were then Ygr determined. Furthermore, in the calculation, the ‘virtual displacement’ or the ‘attempted displacement’, as named by Pauwels (1980), of the centroid of rigid X bodies was also obtained. The direction of ‘virtual xg1 xp xgt displacement’ was found to correspond with the direction of maximum pressure and could be a useful Fig. 2. Global and local coordinates on the two rigid bodies. parameter in describing joint stability. Reference points g, and g2 for each rigid body can be defined Formulation
of rigid-body spring models
The formulation and procedure of the rigid-body spring model for joint contact pressure distribution was illustrated by using a simple hemicircular joint (Fig. 1). Both the convex and concave components of the joint articulation were assumed to be. rigid elements. On the articular surfaces, linear springs of a given stiffness were distributed. Two kinds of springs were considered, normal and shear. The contour of the surface contact area was approximated with any number of polygonal sides. An arbitrary point, usually the center of gravity or location of external force, was taken as a reference point for each rigid-element. Member stiffness matrices for each contact side were calculated and then summed up into the global stiffness matrix for the whole model. Initially, all the possible contact areas were assumed to be in contact. There are six degrees of freedom in this particular twodimensional, two-body contact problem. All of the member stiffness matrices were overlaid on the same six-by-six matrix. On every contact surface (Fig. 2), in addition to the global coordinate system (X, Y), a local coordinate
at the centroids. Point p is located on the contact surface before loading.
system (x, y) was considered. Because of the rigid body assumption, the displacement, UP, of point p on the contact surface of bodies 1 and 2 to p1 and pz was expressed as a function of the displacement of the selected reference point U, (Fig. 3).
U,=CQlU,
(1)
where
and the transformation
r
matrix
1 0
(-Y,-y,,)
0 1
(X,-X&
0 0
0
0
0
0
1
1 0 -cyp-r,,, 0 1 (x,-x,,) where (X,, Y,) and (X,, , Y,,) and (X,, and Y,,) are
1015
Pressure distribution on articular surfaces Y
f/
X
I-
Fig. 3. After loading, the amount of displacement of adjacent rigid bodies can be expressed by the amount of translation (U,I, U,,, V,i, V,*) of the reference points and rotation (C3, and e,). The associated contact point, p, moved top, and pz, on rigid bodies 1 and 2 by amounts of U,,, V,,, and U,,, VP, respectively.
Fig. 4. The relative displacement between the matching points p, and p2 can be expressed on the local coordinate system (x, y). The components normal and tangent to the contact surface are b,, and 6, respectively.
coordinates of point p and reference points, gl and 62, respectively, based on the global coordinate system before deformation. For each rigid body, the transformation between the global coordinate, U,, to the local coordinate, u,, was achieved through a transformation matrix, k.
where
“p =
CR1 up
where I1
CR1 =
and I,, (x, Xx The points
ml
0
0
0
1, ml
1
[M-j=
[
-l
0
0
-10
1
0 1
1
and 6, and 6, represent the components normal and tangent to the contact surface respectively. Combining equations (l), (2), and (3), we have
~=C~lCRICQIU,=C~lU,. Assuming that the normal and shear springs have the linear spring constant, k,, and k,, respectively, the strain energy store, W, is ” W=; !s&T[D]6dS J =;Ugr
CBITIDICBldSU, I
S
k
[B]‘[D][B]dSU,= s
(3)
where
[ 0
and the integration takes place along the active possible contact surfaces, S. Based upon the principle of minimum potential energy, the external force, F, related to the member stiffness matrix, K, was formu-
F=g=
0 0 1, m2 [ ml, I,, and m2 are directional cosines between 6, Y), (Y, X), and (Y, Y), respectively. relative displacement between the matched pi and pl (Fig. 4) was
s=[s,s,]r=[M]UP
kd 01
lated:
0 0
I2 m, 0
(2)
D=
(4)
Is
[K]U,.
(5)
This member stiffness matrix represents two-dimensionally distributed springs along a contact surface, and hence not only normal force balance but also shear and moment balances were maintained. Boundary conditions were then specified to complete the stiffness equation. In this hemicireular joint model (Fig. 1) the concave element was rigidly fixed. A load was applied on the reference point of the convex element. A rotatiotlal degree of freedom of the latter element was also fixed to stabilize the model. The stiffness equation was solved to determine displacements of reference points of each rigid body. For a small number of articulating bodies, and thus a low degree of freedom, the equation can be solved very quickly. This step is usually most time consuming in ordinary finite element analysis. The relative displacement between origin and insertion of a spring was calculated from the location and displacement of involved reference points. Stress on a spring was then calculated by multiplying its elongation with the associated stiffness. Inadmissible solutions on the spring, such as tension on cartilage springs and compression on ligament springs, were removed from the model by subtracting their member stiffness matrices from the global stiffness matrix and repeating the process until no
K. N. AN et al.
1016
inadmissible stress was detected. In most cases, less than five iterations wereenough to obtain the final solution. Determination of spring constant
In the above-described model, the material properties of both cartilage and capsule-ligamentous structures were represented by two spring constants, k, and k,. The values of these constants were determined systematically. For cartilage, on the contact surface S, the normal and shear stresses u and T, are related to the normal and shear strain components, E and y respectively. CJ=
E(l -v)
E and
(l+v)(l-2v)
T=G~
(6)
where E, G, v are the compression modulus of elasticity, shear modulus of elasticity, and Poisson’s ratio of the cartilage, respectively. Assume the undeformed cartilage has thickness h, the strain components Eand y can be approximated by the following finite difference expression: 8=6,/h
and
y=&/h.
Parameters affecting pressure distribution
To examine the effects of various parameters defining the joint contact surfaces as well as the loading configuration, the articulating pressure distribution between congruent circular surfaces of 10 mm radius and 1 mm thickness was considered. As shown in Fig. 1, the arc length of the available articulating surface can be expressed by the subtended angle of the arc, $. The location of the externally applied force with reference to the joint contact surface can be expressed by an angle, 0, measured from the center of the circular arc. In the calculation, a spring constant k, of 1 N mm- 1mmW2 and applied load of 10 N were assumed. The output of the calculation from the model consists of the magnitude of the pressure at any given point on the articulating surface as well as the direction and amount of the relative displacement between the two rigid bodies (Fig. 5). This displacement or motion, U, is referred to as the ‘virtual displacement’ or the ‘attempted displacement’ of the articulating body which also corresponds to the location where
(7)
On the other hand, the following relations are obtained from the definition of spring constants: c= k,6,
and
t=k,6,.
(8)
Therefore, the spring constants could be derived as kd=
E(l -v) (1 +v)(l-2v)h
and
k,=E.
(9)
Normally, in normal human articular joints, friction is negligible. In such conditions, the effect of shear stress on the contact surface could be neglected, and the value for k, is assigned as zero. In other words, the tangential component of spring deformation will not contribute to the spring force calculation. The spring constants for the capsule-ligamentous structure were determined in the same fashion. In this case, the spring constant was determined by dividing the tensile modules of the material by the initial length of the soft tissue under unloaded conditions.
ol~‘l”““““““~J 10
*
180
RESULTS
In this study, application of the formulated model to two problems of joint articulating force distribution was considered. Firstly, the model was applied to a special case of hinge joint. The pressure distribution between two ideal congruent circular arc surfaces were studied so that the results could be compared with those available in the literature and with those using the conventional finite element method. In addition, the effects of characteristic parameters on the magnitude and pattern of pressure distribution could be examined.
Fig. 5. Peak pressure, P,,,, on articular surface and direction of ‘virtual displacement’ or ‘attempted displacement’, U, of articulating bodies are related to available contact surface represented by enclosing angle, +, and location, 0, of the resultant joint force on the hemispheric articular surface (bottom). The angle I#Irepresents the direction of the ‘virtual displacement’ as measured from the center line of the circular arc. The effect of $ with O=o” is shown on the top and the effect of 0 with $ = 180” is shown on the bottom.
Pressure distribution on articular surfaces
1017
dramatically as the subtended angle decreased below 30”. The effects of off-ceoter loading of the pressure distribution were also studied by considering a hemicircular supporting arc. The magnitude of the peak pressure and the direction of the ‘virtual displacement’ as functions of the eccentricity of the load, expressed by angle 0, were obtained (Fig. 5, bottom). Again, the magnitude of the possible peak pressure increased dramatically as the point of load application approached the edge of the supporting arc. The direction P=PDI*cosa of the intended motion or virtual displacement of the where P, is the magnitude of the possible peak rigid body coincided with the direction of external load application when the load was centrally located. pressure. The results were also examined and compared with However, as the load was applied eccentrically, the direction of intended motion deviated from the directhose of analytical elasticity solution by using convention of load application and came closer to the edge of tional finite element methods. A two-dimensional finite element model of circular surface contact (Fig. 1) the arc. With increased eccentric loading, the deviation between these two directions increased, leaving was formulated by using ANSYS package (Swanson Analysis Systems, Inc., Houston, PA, U.S.A.). A layer the possibility for the direction of the intended motion of elastic cartilage surface on each of the rigid bones to lie beyond the articulating arc. was modeled by using two-dimensional isoparametric and triangular solid elements. At the contact surface, Implication in muscle force distribution Application of the above analysis to human joints the two-dimensional interface element was used. This element allows the two surfaces to either maintain or leads to interesting observations. The pattern of break physical contact, and they may slide relative to muscle activity determines the direction and magnieach other. In this simulation this element was capable tude of the resultant joint force on the articulating of supporting only compression in the direction nor- surface and therefore affects the pattern of the joint pressure distribution. In this study, the analysis was mal to the surface, and no shear force in the tangential direction. The normal force in this element is related to applied to the elbow joint pressure distribution in a the joint contact pressure. The values of peak pressure simplified sag&al plane, assuming that a simple force, P, is perpendicularly applied at the distal end of the as well as the pressure distribution on the contact surface using the finite element model were identical to forearm (Fig. 6). An extension moment was therefore generated at the elbow joint. In general, a group of those of using the rigid-body spring model. The effect of the amount of available supporting arc flexor muscles is thought to provide forces for balanlength on the pressure distribution was studied for the cing this extension moment. However, the purpose of case of applied load at the center of the arc in which this part of the study was to examine the pressure case the peak or maximum pressure was also located distribution if the antagonistic extensor were to be recruited. The coronoid and olecranon facets and the at the center of the arc. However, the magnitude depends very much upon the available arc length for troclilea were assumed to be articulated on two segsupport (Fig. 5, top). This peak pressure increased ments of the circular arc as shown.
potential peak pressure may take place. The direction of this ‘virtual displacement’ may have significant implication on the joint stability. For this ideal circular arc contact, the pressure distribution on the surface calculated by the model was identical to those analytic results reported in the literature (Pauwels, 1980) in which the pressure at a gven location on the contact surface applied at an angle a from the point of peak pressure was expressed as a function of the angle
Fig. 6. Pressure distribution on elbow joint surface as external load, P. is applied at the distal end of the ulna. Distribution of muscle force, Fe and F,, influences the magnitude, R, and the direction, 0, of mukant force on the elbow joint. U represent the ‘virtual displacement’ of the humerus relative to the ulna.
1018
K. N.
AN et al.
1
1
I
I
I
.2
.3
.4
.5
.6
Fe
/Ff
Pm (N/Cm2)
10 -
- -40 I
1
I
1
I
.2
.3
.4
.5
.6
Fe/
Ff
Fig. 7. For the given loading condition in Fig. 6, resultant joint force increased with increasing involvement of extensor muscle as represented by the ratio of extensor force F, to flexor force FI (top). Peak articular pressure and the direction of ‘virtual displacement’ of the humerus are also affected by the level of involvement of extensor muscle (bottom).
The distribution of the simplified resultant muscle forces between the flexor, F,, and extensor F,, and thus the resultant joint forces, to counterbalance the externally applied extension moment, was calculated through simple free-body analysis. There are infinite possibilities of muscle force distributions among these two muscles in resisting loads, P, applied to the forearm. For each of the muscle force distributions, a specific magnitude and location of the resultant joint force on the articulating surfaces will result. The joint pressure distribution was thereafter calculated by using the established model. When one unit of force (P= 1 N) was applied on the distal forearm, the magnitude of the resultant joint force, R, and its location on the articular surface, 0, (shown in Fig. 6) were calculated (Fig. 7, top) for different distributions of muscle forces as represented by the ratio between the forces of extensor to flexor. If the externally applied extension moment was resisted by merely the agonistic muscle, the magnitude of resultant joint force, R, was found to be minimal. On
the other hand, with additional involvement of the antagonistic extensor muscle, F,, the magnitude of resultant joint force increased proportionally. However, due to the directions of the lines of action of these two muscle groups as shown in Fig. 6, the associated direction of-the resultant joint force was oriented closer to the center of the trochlear notch with involvement of the extensor. Based on the magnitude and location of this resultant joint force, the direction of the associated ‘virtual displacement’ (U in Fig. 6), and the magnitude of peak joint pressure, Pm, as a function of the relative involvement of the extensor and flexor were also obtained (Fig. 7, bottom). The magnitude of possible peak pressure depends not only on the magnitude of the resultant joint force, but also on how much articular surface was available to support this force. For example, when the resultant joint force is located closer to the edge of the trochlear notch, in this case only the flexor being active, the available contact area was relatively small
Pressure distribution on articular surfaces compared with that when the resultant force was centrally located. Therefore, the magnitude of peak pressure was relatively high even when the magnitude of the overall resultant joint force was low. In this particular example, the ideal articular pressure distribution, with minimum peak pressure, could probably be achieved with activation of both extensor and flexor muscles at a ratio in the range of 0.4-0.5.
DISCUSSlON
Modeling
Estimation of the force distribution on the articular surfaces in the joint is important for better understanding of the etiology of articular cartilage wear and possible osteoarthritis. In addition, we believe it also has some important implications to joint stability. Numerous experimental and analytic methods have been employed in an attempt to obtain such information. In this study, the method based on the rigid-body spring model of Kawai and Takeuchi (1981) was adopted for such analysis. This method was originally developed in civil engineering for analysis of foundation strengthening methods where the magnitude of inter-element displacement is much larger than the intraelement deformation. This method has been widely used and validated in that field. The rationale of adopting this method to determine the force distribution in human articulating joints is that the deformation of the cartilage on the joint surface and the movement or displacement of articulating bone are relatively larger than the actual deformation of the bone. The force-distribution problem, therefore, could be considered by assuming the rigid bony elements to be connected by both compressive and shear springs for simulating cartilage, and tensile springs for simulating capsule and ligamentous structures. When using the model, it is necessary to know the values of the spring constants in order to simulate both the articular cartilage and the capsuloligamentous elements if they are all considered in the model. The precise values of the spring constants can be estimated from the experimental data. However, since the model deve@ed so far is a linear one, only the relative scales of these spring constants are important. Using the model of a congruent circular supporting arc, the effects of the material properties were also examined. The material properties affect the magnitude of the displacement but not the pressure distribution. In addition, a special case of ‘sticky contact surface’, where the tangential deformation is important for the spring force determination, was also studied. With a large variation of Poisson’s ratio, the resultant pressure distributions are almost unchanged. A similar observation was also noticed using the conventional finite element method. The results of the articulating force distribution analysis match quite well with those of Pauwels’ using
1019
the analytic solution of the ideal circular arc as well as those predicted based on the conventional finite element method. While the model of this study in its current state of development compared to Pauwols’ model is more general, it could be used for dotermining articulating pressure distribution of any complicated geometry, such as an elliptical arc or any anatomical geometry. This model could also easily be applied for conditions of discrete segmental articulating arcs such as those of the olecranon and coronoid facets in the elbow joint. Joint pressure
and stability
The model has been utilized to study the transverse stability carpal arch in the wrist. The predicted displacement (Garcia-Elias et u1., 1989a) of the carpal structure under normal and various simulated pathological conditions match quite well with those of experimentally measured (Garcia-Elias et al., 1989b). In addition, the method developed has also been applied for the prediction of intercarpal force transmission (Horii et al., 1990). Again, the general results agree quite favorably with those experimental results reported in the literature. It is therefore felt that the rigid-body spring model formulated could be a reliable tool for the applications discussed. It is clearly illustrated that, if the resultant joint force is centrally located on the articular surface, the pressure distribution about a small arc will be almost uniform. The ‘virtual displacement’ will be close to the direction of the resultant force. The peak pressure decreases with the increase of available contact area. However, when the location of the resultant joint force moves away from the center of possible contact area, the magnitudes of peak pressure increase. In addition, the directiou of ‘virtual displacement’ moves away from the direction of the resultant joint force and toward the rim of the joint contact surface. It is our hypothesis that, for a given loading condition, various types of muscle force distributions could load the joint into three possible conditions. Fig. 8(a) shows a stable condition, when both the resultant joint force and direction of ‘virtual displacement’ are within the arc of contact area. In this case, the joint surfaces will be articulated in their arc of support with
,i3 c3 R a
‘+
R
c
Fig. 8. Based on the relationship of the directions of resultant joint force, R, and ‘virtual displacement’, U, of the aritculating body, loading on the joint could lead to (a) stable, (b) subluxation, (c) dislocation conditions.
1020
K. N. AN et al.
minimum displacement; (b) shows a subluxation condition, when resultant joint force vector is within, but the direction of ‘virtual displacement’ is outside, the articular surface. In this case, the joint surface will have a tendency to offset and a loss of part of the previously available contact surface; and (c) illustrates a dislocated condition, when the resultant joint force vector is pointing outside the articular surface. In this case, the resultant joint force will create a moment to dislocate the joint by pivoting around the rim of the joint surface. Usually, resultant joint force is minimal when there is no involvement of antagonistic muscles, and it progressively increases with the addition of antagonistic muscles. However, for the given configuration in the example, involvement of antagonistic muscles may favor the peak pressure and direction of ‘virtual displacement’. This is because the line of action of resultant muscle force is modified by adding antagonistic muscle, thus bringing the resultant joint force more toward the center of the available contact surface. Acknowledgement-This study is funded in part by NIH grant AR 26287 and CA 40583.
REFERENCES
Brown, T. D. and DiGioia, A. M. (1984)A contact-coupled finite element analysis of the natural adult hip. J. Biomechanics 6,437448.
Garcia-E& M., An, K. N., Cooney, W. P., Linscheid, R. L. and Chao, E. Y. (1989a)Transverse stability of the carpus. An analytic study. J. orthop. Res. 7, 738-743. Garcia-Elias, M., An, K. N., Cooney, W. P., Linscheid, R. L. and Chao, E. Y. (1989b) Stability of the transverse carpal arch. An experimental study. J. Hand Surg. (Am.) 14, 277-282. Greenwald, A. S. and O’Connor, J. J. (1971) The transmission of load through the human hip joint. J. Biomechanics 4, 507-528. Horii, E., Garcia-Elias, M., An, K. N., Bishop, A. T., Cooney, W. P., Linscheid, R. L. and Chao, E. Y. (1990) Effect on force transmission across the carpus of procedures designed to treat Keinb6ck’s disease. An analytic study. J. Hand Surg. (Am.) (in press). Kawai, T. (1977) A new element in discrete analysis of plane strain problems. Seisan Kenkyu 29, 2W207. Kawai, T. and Takeuchi, N. (1981) A discrete method of limit analysis with simplified elements. AXE Int. Con& Cornput. Ciu. Engng. New York. Pauwels, F. (1980) Biomechanics of the Locomotor Apparatus. Springer, Berlin. Wynarsky, T. T. and Greenwald, A. S. (1983) Mathematical model of the human ankle joint. J. Biomechanics 16, 241-251.
J. Biowchanics Vol. 23, No. 10, PP. 1013-1020, 1990. Printed in Great Britain
PRESSURE DISTRIBUTION ON ARTICULAR SURFACES: APPLICATION TO JOINT STABILITY EVALUATION K. N. AN, S. HIMENO, H. TSUMURA, T. KAWAI and E. Y. S.
CHAO
Biomechanics I ‘aboratory, Department of Orthopedics, Mayo Clinic/Mayo Foundation, Rochester, MN 55905, U.S.A. Abutract-A generalized method, based on the rigid body spring model, is developed to determine the twodimensional joint contact pressure of any shape of articular surface and loading condition. The formulation of the model and procedures for the determination of parameters used in the model are presented. The model was used first to analyze the pressure distribution of an ideal hinge joint for comparison with reported results in the literature, and with those using the conventional finite element method. Implications of the results based on the calculated ‘virtual displacement’ to the joint stability are hypothesized and discussed. The problem of joint pressure distribution of the elbow joint is then analyzed, along with the discussion of agonistic and antagonistic muscle activities.
as a result of the analysis itself, and this distribution was treated as external load on each of the joint Precise knowledge of pressure distribution in joints is components. not only a prerequisite for the determination of stresses Brown and DiGioia (1984) discussed a contactin articular cartilage, but might also be the determincoupled finite element analysis of the hip joint which ing factor of the joint stability and thus influences the could calculate the contact force distribution and distribution of muscle force under synergistic and internal stress simultaneously, given the incremental antagonistic conditions. displacements of the acetabular component as the The concept of joint pressure distribution and its boundary condition and iterative calculation to deterapplication was discussed extensively by Pauwels mine the surface members in contact. This was, un(1980) by considering a simple, hinged elbow joint. doubtedly, a great improvement of the situation; Mathematical analyses of force transmission were however, there were still several problems. A long carried out by Greenwald and O’Connor (1971) computing time was required to accomplish this anathrough the hip joint and by Wynarsky and Greenlysis, even though a rather coarse zoning was adopted wald (1983) through the ankle joint. In these two to lessen the amount of calculation. Displacements mathematical analyses, it was assumed that deformawere given rather than loads as boundary conditions. tion took place only in the cartilage layer, treated as a This feature may face some difficulties when the simple elastic material, and that displacement of a method is applied to multibody contact problems, point was independent of the normal stress applied at such as in the wrist or foot joints, because the a priori any other point. They concluded that the distribution displacements of intermediate components cannot be of pressure depended critically on variations in cartil- given. age thickness, and that it also depended on the initial In certain cases where the actual cartilage stress incongruity of the joint at low loads. In their analysis, distribution is of little direct interest, as opposed to the ideal joint geometry was assumed. However, for joints contact stress on the articulating surface, a simpler of irregular configuration, these analytic methods are approach could be performed. In this study, a gendifficult to apply and additional assumptions are eralized method to solve such problems of joint presrequired. sure distribution was developed. The derivation of this Finite element analysis has thus been adopted. It is model was based on the minimum strain energy a well established method. Geometrical complexities principle as used in most of the formulation of finite are easily solved by using fine mesh division. Bound- element codes. An analytic and numeric model was ary conditions are easy to modify, and most of the formulated by using the rigid body spring method of nonlinearities can be handled. Flexibility in modeling Kawai (1977) and Kawai and Takeuchi (1981). The is excellent in this method. As for articulating contact, results were first verified by comparison with those however, most studies assumed the contact force analytic results in the literature for the special case of a distribution on the joint surface a priori rather than hinge joint (Pauwels, 1980). The effects of size of the articular surface and the location of the resultant joint force on the pressure distribution were then examined. Receivedin final form 21 March 1990. The implication of the results on the joint stability as Correspondence to: Kai-Nan An, Ph.D., Orthopedic Bio- well as their influence on the activities of ago&tic and mechanics Lab., Mayo Clinic/Mayo Foundation, 200 First Street, Southwest, Rochester, MN 55905, U.S.A. Tel.: (507) antagonistic muscles are also discussed by considering the elbow joint. 284-2589. 1013
K. N.
1014
AN et al.
METHODS
To establish the analytic and numeric model for joint pressure determination, the shapes of articular surfaces were either mathematically or numerically described. The rigid body spring model was then formulated by considering the rigid body to be in equilibrium with external loading. Reaction forces between adjacent bodies were produced by the spring system distributed over the possible contact surface between the two adjacent bodies. Displacements of the simulated compressive spring, normal and tangential to the articular surface, were described as functions of the displacements of the centroid and rotation of the associated rigid bodies. The capsule and ligaments Fig. 1. Definition of parameters used for joint pressure diswere also modeled with tensile springs, whose defor- tribution study of an ideal circular arc. The available articumation could also be described based on the displace- lating surfaces are expressed by the subtended angle of the ment of the centroids, and rotation of the rigid bodies arc $. The externally applied force, R, is located at angle 0 from the center of the circular arc. as well. Strain energy due to relative displacement of the spring system was formulated using a quadratic function of the displacement vectors of the centroid of Y adjacent rigid bodies. By applying the principle of minimum potential energy, the stiffness equation was obtained by differentiating strain energy with respect to displacement of centroids of adjacent rigid bodies. A system of equations was obtained for solving Y.&T* centroid displacement of rigid bodies. Displacement YP and force in the spring system which describe the pressure distribution on articular surfaces were then Ygr determined. Furthermore, in the calculation, the ‘virtual displacement’ or the ‘attempted displacement’, as named by Pauwels (1980), of the centroid of rigid X bodies was also obtained. The direction of ‘virtual xg1 xp xgt displacement’ was found to correspond with the direction of maximum pressure and could be a useful Fig. 2. Global and local coordinates on the two rigid bodies. parameter in describing joint stability. Reference points g, and g2 for each rigid body can be defined Formulation
of rigid-body spring models
The formulation and procedure of the rigid-body spring model for joint contact pressure distribution was illustrated by using a simple hemicircular joint (Fig. 1). Both the convex and concave components of the joint articulation were assumed to be. rigid elements. On the articular surfaces, linear springs of a given stiffness were distributed. Two kinds of springs were considered, normal and shear. The contour of the surface contact area was approximated with any number of polygonal sides. An arbitrary point, usually the center of gravity or location of external force, was taken as a reference point for each rigid-element. Member stiffness matrices for each contact side were calculated and then summed up into the global stiffness matrix for the whole model. Initially, all the possible contact areas were assumed to be in contact. There are six degrees of freedom in this particular twodimensional, two-body contact problem. All of the member stiffness matrices were overlaid on the same six-by-six matrix. On every contact surface (Fig. 2), in addition to the global coordinate system (X, Y), a local coordinate
at the centroids. Point p is located on the contact surface before loading.
system (x, y) was considered. Because of the rigid body assumption, the displacement, UP, of point p on the contact surface of bodies 1 and 2 to p1 and pz was expressed as a function of the displacement of the selected reference point U, (Fig. 3).
U,=CQlU,
(1)
where
and the transformation
r
matrix
1 0
(-Y,-y,,)
0 1
(X,-X&
0 0
0
0
0
0
1
1 0 -cyp-r,,, 0 1 (x,-x,,) where (X,, Y,) and (X,, , Y,,) and (X,, and Y,,) are
1015
Pressure distribution on articular surfaces Y
f/
X
I-
Fig. 3. After loading, the amount of displacement of adjacent rigid bodies can be expressed by the amount of translation (U,I, U,,, V,i, V,*) of the reference points and rotation (C3, and e,). The associated contact point, p, moved top, and pz, on rigid bodies 1 and 2 by amounts of U,,, V,,, and U,,, VP, respectively.
Fig. 4. The relative displacement between the matching points p, and p2 can be expressed on the local coordinate system (x, y). The components normal and tangent to the contact surface are b,, and 6, respectively.
coordinates of point p and reference points, gl and 62, respectively, based on the global coordinate system before deformation. For each rigid body, the transformation between the global coordinate, U,, to the local coordinate, u,, was achieved through a transformation matrix, k.
where
“p =
CR1 up
where I1
CR1 =
and I,, (x, Xx The points
ml
0
0
0
1, ml
1
[M-j=
[
-l
0
0
-10
1
0 1
1
and 6, and 6, represent the components normal and tangent to the contact surface respectively. Combining equations (l), (2), and (3), we have
~=C~lCRICQIU,=C~lU,. Assuming that the normal and shear springs have the linear spring constant, k,, and k,, respectively, the strain energy store, W, is ” W=; !s&T[D]6dS J =;Ugr
CBITIDICBldSU, I
S
k
[B]‘[D][B]dSU,= s
(3)
where
[ 0
and the integration takes place along the active possible contact surfaces, S. Based upon the principle of minimum potential energy, the external force, F, related to the member stiffness matrix, K, was formu-
F=g=
0 0 1, m2 [ ml, I,, and m2 are directional cosines between 6, Y), (Y, X), and (Y, Y), respectively. relative displacement between the matched pi and pl (Fig. 4) was
s=[s,s,]r=[M]UP
kd 01
lated:
0 0
I2 m, 0
(2)
D=
(4)
Is
[K]U,.
(5)
This member stiffness matrix represents two-dimensionally distributed springs along a contact surface, and hence not only normal force balance but also shear and moment balances were maintained. Boundary conditions were then specified to complete the stiffness equation. In this hemicireular joint model (Fig. 1) the concave element was rigidly fixed. A load was applied on the reference point of the convex element. A rotatiotlal degree of freedom of the latter element was also fixed to stabilize the model. The stiffness equation was solved to determine displacements of reference points of each rigid body. For a small number of articulating bodies, and thus a low degree of freedom, the equation can be solved very quickly. This step is usually most time consuming in ordinary finite element analysis. The relative displacement between origin and insertion of a spring was calculated from the location and displacement of involved reference points. Stress on a spring was then calculated by multiplying its elongation with the associated stiffness. Inadmissible solutions on the spring, such as tension on cartilage springs and compression on ligament springs, were removed from the model by subtracting their member stiffness matrices from the global stiffness matrix and repeating the process until no
K. N. AN et al.
1016
inadmissible stress was detected. In most cases, less than five iterations wereenough to obtain the final solution. Determination of spring constant
In the above-described model, the material properties of both cartilage and capsule-ligamentous structures were represented by two spring constants, k, and k,. The values of these constants were determined systematically. For cartilage, on the contact surface S, the normal and shear stresses u and T, are related to the normal and shear strain components, E and y respectively. CJ=
E(l -v)
E and
(l+v)(l-2v)
T=G~
(6)
where E, G, v are the compression modulus of elasticity, shear modulus of elasticity, and Poisson’s ratio of the cartilage, respectively. Assume the undeformed cartilage has thickness h, the strain components Eand y can be approximated by the following finite difference expression: 8=6,/h
and
y=&/h.
Parameters affecting pressure distribution
To examine the effects of various parameters defining the joint contact surfaces as well as the loading configuration, the articulating pressure distribution between congruent circular surfaces of 10 mm radius and 1 mm thickness was considered. As shown in Fig. 1, the arc length of the available articulating surface can be expressed by the subtended angle of the arc, $. The location of the externally applied force with reference to the joint contact surface can be expressed by an angle, 0, measured from the center of the circular arc. In the calculation, a spring constant k, of 1 N mm- 1mmW2 and applied load of 10 N were assumed. The output of the calculation from the model consists of the magnitude of the pressure at any given point on the articulating surface as well as the direction and amount of the relative displacement between the two rigid bodies (Fig. 5). This displacement or motion, U, is referred to as the ‘virtual displacement’ or the ‘attempted displacement’ of the articulating body which also corresponds to the location where
(7)
On the other hand, the following relations are obtained from the definition of spring constants: c= k,6,
and
t=k,6,.
(8)
Therefore, the spring constants could be derived as kd=
E(l -v) (1 +v)(l-2v)h
and
k,=E.
(9)
Normally, in normal human articular joints, friction is negligible. In such conditions, the effect of shear stress on the contact surface could be neglected, and the value for k, is assigned as zero. In other words, the tangential component of spring deformation will not contribute to the spring force calculation. The spring constants for the capsule-ligamentous structure were determined in the same fashion. In this case, the spring constant was determined by dividing the tensile modules of the material by the initial length of the soft tissue under unloaded conditions.
ol~‘l”““““““~J 10
*
180
RESULTS
In this study, application of the formulated model to two problems of joint articulating force distribution was considered. Firstly, the model was applied to a special case of hinge joint. The pressure distribution between two ideal congruent circular arc surfaces were studied so that the results could be compared with those available in the literature and with those using the conventional finite element method. In addition, the effects of characteristic parameters on the magnitude and pattern of pressure distribution could be examined.
Fig. 5. Peak pressure, P,,,, on articular surface and direction of ‘virtual displacement’ or ‘attempted displacement’, U, of articulating bodies are related to available contact surface represented by enclosing angle, +, and location, 0, of the resultant joint force on the hemispheric articular surface (bottom). The angle I#Irepresents the direction of the ‘virtual displacement’ as measured from the center line of the circular arc. The effect of $ with O=o” is shown on the top and the effect of 0 with $ = 180” is shown on the bottom.
Pressure distribution on articular surfaces
1017
dramatically as the subtended angle decreased below 30”. The effects of off-ceoter loading of the pressure distribution were also studied by considering a hemicircular supporting arc. The magnitude of the peak pressure and the direction of the ‘virtual displacement’ as functions of the eccentricity of the load, expressed by angle 0, were obtained (Fig. 5, bottom). Again, the magnitude of the possible peak pressure increased dramatically as the point of load application approached the edge of the supporting arc. The direction P=PDI*cosa of the intended motion or virtual displacement of the where P, is the magnitude of the possible peak rigid body coincided with the direction of external load application when the load was centrally located. pressure. The results were also examined and compared with However, as the load was applied eccentrically, the direction of intended motion deviated from the directhose of analytical elasticity solution by using convention of load application and came closer to the edge of tional finite element methods. A two-dimensional finite element model of circular surface contact (Fig. 1) the arc. With increased eccentric loading, the deviation between these two directions increased, leaving was formulated by using ANSYS package (Swanson Analysis Systems, Inc., Houston, PA, U.S.A.). A layer the possibility for the direction of the intended motion of elastic cartilage surface on each of the rigid bones to lie beyond the articulating arc. was modeled by using two-dimensional isoparametric and triangular solid elements. At the contact surface, Implication in muscle force distribution Application of the above analysis to human joints the two-dimensional interface element was used. This element allows the two surfaces to either maintain or leads to interesting observations. The pattern of break physical contact, and they may slide relative to muscle activity determines the direction and magnieach other. In this simulation this element was capable tude of the resultant joint force on the articulating of supporting only compression in the direction nor- surface and therefore affects the pattern of the joint pressure distribution. In this study, the analysis was mal to the surface, and no shear force in the tangential direction. The normal force in this element is related to applied to the elbow joint pressure distribution in a the joint contact pressure. The values of peak pressure simplified sag&al plane, assuming that a simple force, P, is perpendicularly applied at the distal end of the as well as the pressure distribution on the contact surface using the finite element model were identical to forearm (Fig. 6). An extension moment was therefore generated at the elbow joint. In general, a group of those of using the rigid-body spring model. The effect of the amount of available supporting arc flexor muscles is thought to provide forces for balanlength on the pressure distribution was studied for the cing this extension moment. However, the purpose of case of applied load at the center of the arc in which this part of the study was to examine the pressure case the peak or maximum pressure was also located distribution if the antagonistic extensor were to be recruited. The coronoid and olecranon facets and the at the center of the arc. However, the magnitude depends very much upon the available arc length for troclilea were assumed to be articulated on two segsupport (Fig. 5, top). This peak pressure increased ments of the circular arc as shown.
potential peak pressure may take place. The direction of this ‘virtual displacement’ may have significant implication on the joint stability. For this ideal circular arc contact, the pressure distribution on the surface calculated by the model was identical to those analytic results reported in the literature (Pauwels, 1980) in which the pressure at a gven location on the contact surface applied at an angle a from the point of peak pressure was expressed as a function of the angle
Fig. 6. Pressure distribution on elbow joint surface as external load, P. is applied at the distal end of the ulna. Distribution of muscle force, Fe and F,, influences the magnitude, R, and the direction, 0, of mukant force on the elbow joint. U represent the ‘virtual displacement’ of the humerus relative to the ulna.
1018
K. N.
AN et al.
1
1
I
I
I
.2
.3
.4
.5
.6
Fe
/Ff
Pm (N/Cm2)
10 -
- -40 I
1
I
1
I
.2
.3
.4
.5
.6
Fe/
Ff
Fig. 7. For the given loading condition in Fig. 6, resultant joint force increased with increasing involvement of extensor muscle as represented by the ratio of extensor force F, to flexor force FI (top). Peak articular pressure and the direction of ‘virtual displacement’ of the humerus are also affected by the level of involvement of extensor muscle (bottom).
The distribution of the simplified resultant muscle forces between the flexor, F,, and extensor F,, and thus the resultant joint forces, to counterbalance the externally applied extension moment, was calculated through simple free-body analysis. There are infinite possibilities of muscle force distributions among these two muscles in resisting loads, P, applied to the forearm. For each of the muscle force distributions, a specific magnitude and location of the resultant joint force on the articulating surfaces will result. The joint pressure distribution was thereafter calculated by using the established model. When one unit of force (P= 1 N) was applied on the distal forearm, the magnitude of the resultant joint force, R, and its location on the articular surface, 0, (shown in Fig. 6) were calculated (Fig. 7, top) for different distributions of muscle forces as represented by the ratio between the forces of extensor to flexor. If the externally applied extension moment was resisted by merely the agonistic muscle, the magnitude of resultant joint force, R, was found to be minimal. On
the other hand, with additional involvement of the antagonistic extensor muscle, F,, the magnitude of resultant joint force increased proportionally. However, due to the directions of the lines of action of these two muscle groups as shown in Fig. 6, the associated direction of-the resultant joint force was oriented closer to the center of the trochlear notch with involvement of the extensor. Based on the magnitude and location of this resultant joint force, the direction of the associated ‘virtual displacement’ (U in Fig. 6), and the magnitude of peak joint pressure, Pm, as a function of the relative involvement of the extensor and flexor were also obtained (Fig. 7, bottom). The magnitude of possible peak pressure depends not only on the magnitude of the resultant joint force, but also on how much articular surface was available to support this force. For example, when the resultant joint force is located closer to the edge of the trochlear notch, in this case only the flexor being active, the available contact area was relatively small
Pressure distribution on articular surfaces compared with that when the resultant force was centrally located. Therefore, the magnitude of peak pressure was relatively high even when the magnitude of the overall resultant joint force was low. In this particular example, the ideal articular pressure distribution, with minimum peak pressure, could probably be achieved with activation of both extensor and flexor muscles at a ratio in the range of 0.4-0.5.
DISCUSSlON
Modeling
Estimation of the force distribution on the articular surfaces in the joint is important for better understanding of the etiology of articular cartilage wear and possible osteoarthritis. In addition, we believe it also has some important implications to joint stability. Numerous experimental and analytic methods have been employed in an attempt to obtain such information. In this study, the method based on the rigid-body spring model of Kawai and Takeuchi (1981) was adopted for such analysis. This method was originally developed in civil engineering for analysis of foundation strengthening methods where the magnitude of inter-element displacement is much larger than the intraelement deformation. This method has been widely used and validated in that field. The rationale of adopting this method to determine the force distribution in human articulating joints is that the deformation of the cartilage on the joint surface and the movement or displacement of articulating bone are relatively larger than the actual deformation of the bone. The force-distribution problem, therefore, could be considered by assuming the rigid bony elements to be connected by both compressive and shear springs for simulating cartilage, and tensile springs for simulating capsule and ligamentous structures. When using the model, it is necessary to know the values of the spring constants in order to simulate both the articular cartilage and the capsuloligamentous elements if they are all considered in the model. The precise values of the spring constants can be estimated from the experimental data. However, since the model deve@ed so far is a linear one, only the relative scales of these spring constants are important. Using the model of a congruent circular supporting arc, the effects of the material properties were also examined. The material properties affect the magnitude of the displacement but not the pressure distribution. In addition, a special case of ‘sticky contact surface’, where the tangential deformation is important for the spring force determination, was also studied. With a large variation of Poisson’s ratio, the resultant pressure distributions are almost unchanged. A similar observation was also noticed using the conventional finite element method. The results of the articulating force distribution analysis match quite well with those of Pauwels’ using
1019
the analytic solution of the ideal circular arc as well as those predicted based on the conventional finite element method. While the model of this study in its current state of development compared to Pauwols’ model is more general, it could be used for dotermining articulating pressure distribution of any complicated geometry, such as an elliptical arc or any anatomical geometry. This model could also easily be applied for conditions of discrete segmental articulating arcs such as those of the olecranon and coronoid facets in the elbow joint. Joint pressure
and stability
The model has been utilized to study the transverse stability carpal arch in the wrist. The predicted displacement (Garcia-Elias et u1., 1989a) of the carpal structure under normal and various simulated pathological conditions match quite well with those of experimentally measured (Garcia-Elias et al., 1989b). In addition, the method developed has also been applied for the prediction of intercarpal force transmission (Horii et al., 1990). Again, the general results agree quite favorably with those experimental results reported in the literature. It is therefore felt that the rigid-body spring model formulated could be a reliable tool for the applications discussed. It is clearly illustrated that, if the resultant joint force is centrally located on the articular surface, the pressure distribution about a small arc will be almost uniform. The ‘virtual displacement’ will be close to the direction of the resultant force. The peak pressure decreases with the increase of available contact area. However, when the location of the resultant joint force moves away from the center of possible contact area, the magnitudes of peak pressure increase. In addition, the directiou of ‘virtual displacement’ moves away from the direction of the resultant joint force and toward the rim of the joint contact surface. It is our hypothesis that, for a given loading condition, various types of muscle force distributions could load the joint into three possible conditions. Fig. 8(a) shows a stable condition, when both the resultant joint force and direction of ‘virtual displacement’ are within the arc of contact area. In this case, the joint surfaces will be articulated in their arc of support with
,i3 c3 R a
‘+
R
c
Fig. 8. Based on the relationship of the directions of resultant joint force, R, and ‘virtual displacement’, U, of the aritculating body, loading on the joint could lead to (a) stable, (b) subluxation, (c) dislocation conditions.
1020
K. N. AN et al.
minimum displacement; (b) shows a subluxation condition, when resultant joint force vector is within, but the direction of ‘virtual displacement’ is outside, the articular surface. In this case, the joint surface will have a tendency to offset and a loss of part of the previously available contact surface; and (c) illustrates a dislocated condition, when the resultant joint force vector is pointing outside the articular surface. In this case, the resultant joint force will create a moment to dislocate the joint by pivoting around the rim of the joint surface. Usually, resultant joint force is minimal when there is no involvement of antagonistic muscles, and it progressively increases with the addition of antagonistic muscles. However, for the given configuration in the example, involvement of antagonistic muscles may favor the peak pressure and direction of ‘virtual displacement’. This is because the line of action of resultant muscle force is modified by adding antagonistic muscle, thus bringing the resultant joint force more toward the center of the available contact surface. Acknowledgement-This study is funded in part by NIH grant AR 26287 and CA 40583.
REFERENCES
Brown, T. D. and DiGioia, A. M. (1984)A contact-coupled finite element analysis of the natural adult hip. J. Biomechanics 6,437448.
Garcia-E& M., An, K. N., Cooney, W. P., Linscheid, R. L. and Chao, E. Y. (1989a)Transverse stability of the carpus. An analytic study. J. orthop. Res. 7, 738-743. Garcia-Elias, M., An, K. N., Cooney, W. P., Linscheid, R. L. and Chao, E. Y. (1989b) Stability of the transverse carpal arch. An experimental study. J. Hand Surg. (Am.) 14, 277-282. Greenwald, A. S. and O’Connor, J. J. (1971) The transmission of load through the human hip joint. J. Biomechanics 4, 507-528. Horii, E., Garcia-Elias, M., An, K. N., Bishop, A. T., Cooney, W. P., Linscheid, R. L. and Chao, E. Y. (1990) Effect on force transmission across the carpus of procedures designed to treat Keinb6ck’s disease. An analytic study. J. Hand Surg. (Am.) (in press). Kawai, T. (1977) A new element in discrete analysis of plane strain problems. Seisan Kenkyu 29, 2W207. Kawai, T. and Takeuchi, N. (1981) A discrete method of limit analysis with simplified elements. AXE Int. Con& Cornput. Ciu. Engng. New York. Pauwels, F. (1980) Biomechanics of the Locomotor Apparatus. Springer, Berlin. Wynarsky, T. T. and Greenwald, A. S. (1983) Mathematical model of the human ankle joint. J. Biomechanics 16, 241-251.
