- Email: [email protected]
Sensitivity of femoral strain pattern analyses to resultant and muscle forces at the hip joint
Med. Eng. Phys. Vol. 18, No. 1, Pp. 70-78, 1996 Copyright 0 1995 Elsevier Science Ltd for IPEMB
Printed in Great Britain. All rights reserved
Sensitivity of femoral strain pattern analyses to resultant and muscle forces at the hip joint M. Lengsfeld*,
J. Kaminsky”,
B. Merzt
and R. P. Franke’
*Department of Orthopaedic Surgery, Philipps-University of Marburg, Germany; tbrstitute of Biomedical Engineering and Medical Informatics, Swiss Federal Institute of Technology and University of Zurich, Switzerland; IInstitute of Biomedical Engineering, Department of Biomaterials, University of Ulm, Germany Received
14 December
1994, accepted
27 March
1995
AJWIRACT An automated geometrical preprocessor was dareloped with the aim of creating three-dimensional finite element models (IZM) of the human femur. On the basis of postprocessed computed tomography data, this preprocessor makes possible rapid, flexible and regular meshing with ‘brick’ elements. Three different material properties were modelled at the present stage of development. Sensitivity anabses demonstrated that the strain energy dasity (SELI) patterns of the dafferent femoral parts were most sensitive to the implementation of an iliotibial tract force. The variation of the resultant hip force and abductor force direction within the sagittal plane demonstrated a SED minimum at an anterim inclination of 17’; the variation of the resultant force direction within the frontal plane demonstrated a minimum SELI at a medial inclination of 21” relative to the mechanical axis of the lower limb. The orientation of the connecting line between the surface-SEL+eak.s in the horizontal uitw was found to be most sensitive to the variation of the resultant force within the sagittal plane.
Keywords: FEM, Med.
Eng.
Phys.,
preprocessing, 1996,
Vol.
human
femur,
force
The finite element method (FEM) is widely accepted as a powerful and also a rather sensitive numerical tool, which provides a way to simulate biomechanical phenomenala. FEM-software packages have powerful and convenient preprocessors, although not specifically adapted to the requirements of biomechanical modelling. An experienced investigator is required to develop models interactively as realistic as possible and to avoid errors concerning input model parameters and boundary conditions. In fact, irregular geometry, complex microstructure of biological tissues and loading situations are specific prob lems of the application of the FEM in biomechanits and are still difficult to model. Nevertheless, from our clinical point of view, it is desirable to provide specific biomechanical FEM-preprocessors, which should be easy to handle and which could allow a fast model development. This would also be a contribution to a faster investigation of individual variations and, in the future, could allow one to utilize the method in the clinical setting, preoperatively. to: Dr Markus PhilippsUniversity Germany.
sensitivity
analyses
18, 70-78, January
INTRODUCTION
Correspondence paedic Surgery, 35033 Marburg,
application,
Lengsfeld, Department of Ortbw of Marburg. Baldingerstrasse, D-
Therefore, our project has two objectives: (1) the introduction of a new automated FEM-preprocessor in order to create S-D-models of human femora, based upon postprocessed computed tomography data; and (2) a sensitivity study focusing on the influence of force transmission, meshing, model size and material properties on femoral strain energy density (SED) patterns. MATERIALS
AND METHOD
Description of the preprocessor In order to provide geometrical data, the inner and outer contours of the femoral cortical bone were extracted from computed tomograms (CT) with the method of thresholding. As cortical regions in cross-sections near the joints are very thin, modelling of the cortical bone by extracting the inner contour was neglected in these regions ‘,l”. Ten programs were developed, creating, step by step, a solvable 3D-FEM-model of the femur. The parametric design language of the ANSYS 5.0 general purpose FEM-package (Swanson Analysis Systems, Inc. Houston, PA, U.S.A) including ANSYS preprocessing commands was used to write these programs (macros) : 1. Definition
of the inner and outer femoral
sur-
Sasitiuity
2.
3.
4. 5.
6.
face by key-point-based interpolating spline functions. Description of the femoral curvature by a longitudinal mid-femoral line (vertical connecting line). This line represents the connection of the centres of the smallest circles fitting in the outer femoral contour (ZGg~re I). These centres including the centre of the femoral head were worked out slice after slice and were iteratively calculated (Eg~e 2). Input parameters were a tolerance value and a starting position vector within a transverse femoral slice. Definition of lines subdividing the transversal cross-sections into several sectors. The lines are orientated perpendicularly to the longitudinal mid-femoral line. If necessary, an optional program provides a smoothing of the outer femoral surface. Implementation of nodes along the sector lines and deletion of the lines afterwards. The centrally located intramedullary part is subdivided by nodes describing a rectangular geometry (Figure 1). Starting proximally, the predefined nodes
strain pattern
anal~srs:
M. Imgsfeld
et al.
were meshed with 3-D isoparametric 8-nodes brick elements (Figure 3). 7. Detection of distorted elements, which may cause numerical problems during the solution and tetrahedral-like reshaping of these elements without changing of the element 8. x!&ment of material properties to the different parts of the femur. 9. Definition of the boundary conditions. The nodes at the distal end of the model were constrained. The resultant hip force was applied by multiple force components at the nodes located within the geometrical extent of the load bearing surface l1 . The pelvi-trochanteric muscle force was applied at the tip of the greater trochanter in order to counterbalance the partial body weight in the single-leg stance phasel’-lt’~ The angle of the resultant force to the sagittal and frontal plane, the partial body weight and its lever arm are input parameters. 10. Start of the Solver. The standard model generated and tested here comprised 3308 elements, 3750 nodes and approximately 11,200 degrees of freedom. The thickness of the element layers was 12 mm. A finer mesh was modelled at the proximal femur (4 mm element layer). The calculations were run on an i486 DX 2-66 local bus personal computer and it took 3 h 48 min of computing time to generate the model and 1 h 22 min to run the strain analysis (solution). Input
Figure 1 The femoral curvature is described by a longitudinal midfemoral line. This line represents the connection of the centres of the smallest circles fitting in the outer femoral contour. Sector lines are orientated perpendicularly to the longitudinal mid-femoral line. The cross-section of the trochanteric region after deletion of the mid-femoral and sector lines demonstrates the architecture of the final mesh. A: anterior; L: lateral; M: medial; P: posteriox-
offpmoral
parameters
and description
of the model
Mesh density of the layers of brick elements is determined by the number of lines subdividing the layer into several sectors. The total number of element layers can be defined by changing their thicknesses. The thickness of the layers can be independendy selected from the table-feed of the computed tomograms. Finally, the user selects the Poisson number and three different Young’s moduli for the cortical and trabecular bone and the intramedullar tissue, which were previously separated by contour detection. Young’s moduli are averaged in case of elements located on the border between two differently stiff volumes by the eighth macro. If there is no further information given, the Young’s moduli assigned were 17,000 MPa (cortical bone), 1000 MPa (trabecular bone), 10 MPa (intramedullar tissue) and isotropic mechanical behaviour was assumed in the present study. A Poisson number of 0.33 was selected. About 85% of the femoral length was modelled, and the most distal layer of the model was located above the femoral condyles. The nodes within this plane were constrained in all directions. The model was based on CT scan slices taken from a right fresh-frozen femoral specimen. It was obtained from a male adult of 45 years with a weight of 77 kg”. An antero-posterior X-ray excluded pathological lesions. In case of an integration of the iliotibial tract
71
Sensitivity
offemoral
strain pattern
analyses: M. Lerpfeld
et al.
Starting position vector (C) within a horizontal femoral slice. Predefined tolerance value (a) of the centre to be determined.
1 Minimum distance vector (Dmin) between C and the outer femoral contour. The distance to the opposite outer femoral contour is expressed by S*Dmin.
j ! s/2
H s
Minimum distance (D) of the position vector : C+
s
* Dmin
Final centre : C+S*Dmin
H C
C
Figure 2 Flow chart of the iteration. A mid-femoral line is calculated, describing the 3D curvature of the femur. a: predefined tolerance value; convergence was based on Is*Dminl < a; C: starting position vector and position vector describing the centre of a transverse femoral slice, respectively; Dmin: minimum distance vector between C and the outer femoral contour; s: scalar; the product of the scalar s and the vector Dmin describes the distance vector between C and the opposite outer femoral contour; this vector has a direction opposite to that of Dmin
force in the model, it is applied at the middle of the greater trochanter and acts along the femoral axis according to Rohlmann et aL3 and Heimkes et al . I* . The standard load configuration was a resultant force of 2624 N, a gluteus medius muscle (abductor) force of 2036 N, a partial body weight
72
of 615 N and an iliotibial tract force of 611 N. According to Kummerlg, the latter represented 30% of the gluteus medius muscle force. The lever arm of the partial body weight was 108 mm. Forces, force components and force directions were described in this study relative to an orthog-
head
-__----
--
n&
______---
--
trochaater 1 -----trochanter 2 ------
shaft
2
-------
---
shaft 3 __-_---- --
IL
P-A
x
Jr
Y A
X
M-L
A-P
&aft 5 ____--______
L Y L-M
Figure 3 Element plots of the model. The standard load configuration of a resulrant force of 2624 N, a gluteus medius muscle (abductor) force of 2036 N, a partial body weight of 615 N and an iliotibial tract force of 611 N is displayed. A: anterior; IL lateral: M: medial; P: posterior
onal coordinate system, where its origin coincided with the centre of the femoral ‘head. The vertical axis was defined by the centre of the femoral head and the estimated knee joint centre within the femoral groove. The direction of the frontal axis was defined by the dorsal aspect of the femoral condyles. Relative to this coordinate system and according to the moment equilibrium of the single leg stance, the resultant force was 21” mediolaterally and 13” frontodorsally directed. The femoral shaft axis determined the iliotibial tract force orientation. A neck-shaft angle of 136” and an anteversion angle of 19” were measured at the femur modelled in this study. The results of the calculations testing the influence of force applications on the whole model were approximated (smoothed) by a non-interpolating third-order-polynomial using Excel (Microsoft Corporation, U.S.A.). The FEM-results were expressed by the SED (SED = strain energy/bone volume). It was hypothesized that the SED is suitable for comparative strain analysis of bone20-2’. The femur was divided into ten different parts in order to postprocess the results also relative to particular femoral parts: the head, neck, proximal trochanter and distal trochanter. The shaft was subdivided into six slices (F&WY 4).
shaft 6 -- - -- --
Figure 4 The femur was divided into ten dif’ferent parts in order to postprocess the results also relative to particular femoral parts: the head, neck, proximal trochanteland distal trochantrr. The shaft was \uhdividrd into six slices
RESULTS Mesh
density
In order to test the influence of the mesh density on the SED distribution, the slice thickness within the shaft was modelled with 4, 8, 12 and 16 mm. The slice thickness within the head, neck and trochanteric region was 4 mm. Consequently, the number of elements was observed to be between approximately 2500 and 7000. Almost identical results were found when testing the models with 4, 8, 12 and 16 mm slice thickness. Model
size
Using the standard boundary conditions described above, the model length, i.e. the height of the most distal layer, where all nodes were constrained, was varied. Small influence was found of the height of the constrained layer on the SED in the proximal femoral arts. The complete model is represented by 100 7o and the tests were performed assuming a model length between 100 and 50%. It was observed that only the two most distal layers were not constant. This distal part represented about 10% relative to the complete model.
73
Sensitivity
ojjaoral
Young’s
strain pattern
analyses: M. Lagsjeld
et al.
moduli
Variation of the resultant force and abductor force direction within the sag&al plane
Studying the SED of the bone when varying the Young’s modulus of the intramedullar tissue between 5 and 1000 MPa, there was no significant difference seen. Increasing the Young’s modulus of the cortical bone from 1000 to 20,000 MPa, a decrease of the SED was found and was particularly significant at the distal shaft. The SED ranged from 12 to 1 MPa/103. The variation of the cancellous Young’s modulus between 10 and 10,000 MPa showed a decreasing pattern of the SED at the head (2.5-0.4 MPa/lO ), neck (2.10.6 MPa/lO’) and the distal trochanter (3.22.2 MPa/103) and at the distal layers of the shaft (1.3-0.9 MPa/103). The other parts of the shaft were almost unaffected. Iliotibial
tract
The direction of the resultant force and the abductor force within the sagittal plane were systematically varied between -10” and 35” with an angle step of 5”, except between 0” and 20”, where an angle step of 1” was selected. The amount of the forces and their frontal components are constant. The negative sign refers to a dorso-frontal force direction (posterior inclination), and the positive sign to a fronto-dorsal direction (anterior inclination). In order to present the results (Figure I$‘, the SED of all elements that the model consisted of was averaged and the SED of the elements of various parts was averaged. The averaged SED values of the whole model were approximated by a third-order-polynoma (y:SED[MPa/103]; xangle[“]): y = 0.0002X7 + 0.01152 - 0.395x + 5.864 having a determination coefficient (72) of 0.93. The SED minimum of the polynomial was found at an anterior inclination of 13” and increased significantly at -10 and 35” inclination.
force
The influence of the iliotibial tract force on the distribution of the SED was tested by varying this force between 0 and 100% of the abductor force with a step of 5%. The SED values of the femoral head and neck were not altered by the introduction of the iliotibial tract force. The minimum SED of the distal femoral shaft (shaft 5 and 6) was observed between 0 and 15%. In contrast, the minimum SED of the greater trochanter and the proximal shaft were seen at 100% (trochanter 1 and 2, shaft 1). The minimum SED of the whole model was calculated at 45% (Figure 5).
Variation of the resultant the frontal plane
force
direction
within
The medio-lateral direction of the resultant force within the frontal plane was systematically varied between 7” and 42” with an angle step of 5”, except between 12” and 32’, where an angle step
12 +
head
-.x-neck -
nochanter I
-
trochsnter 2
+
shaft 2
--e-
shaft 3
tshsft4
+ -whole
shaft 6 model
0 0
5
10
15
20
25
30
35
40
45
50
55
60
65
70
75
80
85
90
95
100
Iliotibial tract force [ % of the abductor force ] Figure 5 The influence of the iliotibial tract force [N] on the distribution of the strain engery density (SED [MPa/10”]) was tested by varying this force between 0 and 100% of the abductor force with a step of 5%. The SED values of the femoral head and neck were not altered by the introduction of the iliotibial tract force. The minimum SED of the distal femoral shaft (shaft 5 and 6) was observed between 0 and 15%. In contrast, the minimum SED of the greater trochanter and the proximal shaft were seen at 100% (trochanter 1 and 2, shaft 1)
74
Sensitivity
oj,fpmmnl
@rain patk-rn
ana1y.sP.y: M. i.mgsfild
et al.
r
-+--head +-neck ---.
trochanter
I
--
rrochanter
2
-*-shaft2 -+-shaft
+
3
shaft 4
tshaff6 -whole -polynomial (whole 0
5
15
IO
20
Angle in YZ - plane of th; bone-coordinate
25
30
model
model)
35
system 1 O ]
Figure 6 The direction of the resultant force and the abductor force within the sagittal (yr) plane were systematically varied between -10 and 35”. The negative sign refers to a dorso-frontal force direction (posterior inclination), the positive sign to a fronto-dorsal direr&ion (anrerior inclination). The strain energy density (SED) of all elements which the model consisted of was averaged and the SED of the elements of various parts was averaged. The averaged SED values of the whole model were approximated by a third-order-polynomial (y:SED [MPailO’]; rangie[‘J)
of 1” was selected. An angle of 0” refers to a vertical direction of the resultant force and could not be tested because of a missing equilibrium. The direction of the resultant and the partial body weight were given. In order to match the equilibrium condition, the amount of the resultant force and the amount and direction of the abductor force had to be recalculated for every angle step, as described by Hamacher and Rbsler’ *I’. The resultant force direction within the sag&al plane was 13’. The averaged SED values of the whole model were approximated by a third-order-poly(y:SED[MPa/IO”]; xangle[“]): y = nomial 0.0004$ + 0.00552 - 0.0933x + 4.4487 having a determination coefficient (?) of 0.88. The minimum SED of this polynomial was found at an angle of 21°, and increased significantly at 36” (Fiflrf? 7). Orientation of the connecting lie between surfaceSED-peaks in the transverse view
the
When looking at transverse slices of the model, the medial and lateral elements of maximum SED were selected. At the proximal shaft the connecting line between the SED peaks had a slight antero-medial/postero-lateral orientation, while almost no deviation from the frontal direction was seen in the central &aft. A particular sensitivity of this parameter was found during the variation of the hip load within the sagittal plane. A large anterior inclination caused an internal rotation of the connect.ing line (&&Lw 8).
DISCUSSION Assuming a simplified experimental setup of femoral load application the finite element method has been validated by strain gauge measurements3,“.23,24 and frequently utilized. Moreover, the method has been applied to model the response of bone to a previously calculated stress field, which means a modification of the material property assignments during the remodelling simulation20,?1,2~,2~.
From our point of view, however, there are still two important aspects restricting the applicability and validity of this testing method: (a) generation of the models appears still rather time-consuming and requires engineering expertise making the introduction in the clinical setting and as a clinical screening tool in viuo almost impossible; and (b) the majority of researchers used a proximal load application according to the classical model of Pauwels’* of the single leg stance phase’.“,“. There is still not very much known about the sensitivity of the amount and direction of maximum stresses on the proximal load direction modelling. This and other input parameters have to be carefully tested before finite element analysis can also be safely and reliably utilized by the clinician. In addition to the use of FEM in the future in the clinic, it would still be desirable to learn more about the behaviour and sensitivity of femoral FEM-models before drawing clinical or design consequences from FEM-tests of implants. Based on a contour data file extracted from
Sensitivity
of femoral
strain pat&em analyses: M. Lengsjeld et al.
25 -+-head
-
ttochanter
1
-
trochanter
2
+
shaft 2
-+
shaft 3
+
shaft 4
-o-shaft -a-whole
Od....:....:....:. 5
, 10
15
.
20
.:.
, . .:. 2s
. . 30
Angle in XZ - plane of thd bone-coordinate
.:-, 35
5 shafl6 model
. .:.. 40
system [ O 1
Figure 7 The medio-lateral direction of the resultant force within the frontal (xz) plane was systematically varied between 7” and 42”. The direction of the resultant and the partial body weight were given. In order to match the equilibrium condition, the amount of the resultant force and the amount and direction of the abductor force had to be recalculated for every angle step. The resultant force direction within the sagittal plane was 13”. The averaged SED values of the whole model were approximated by a third-order-polynomial (y:SED [MPa/103]; xangler])
Angle in YZ - Diane of the bone-coordinate
-
system [ o ]
8 b Q 8 Q Q
head
1 t
neck trochanter
1
shaft 4 shaft 5
Q
@
0
Q
Qj
8
QJC)
shaft 6
Q
Q
Q
@
Q
d
a-
-A--l@ .---_-
whole model Figure 8 Orientation of the surfaceSED (strain energy density)-peaks in a schematic transverse 6, the direction of the resultant force and the abductor force within the sagittal (yz) plane were The left side of the circles refers to the medial, and the top of the circles refers to the ventral part the circumference of the dark half refer to the maximum SED (tensile strain) and of the white half
computed tomography, the finite element preprocessor presented here reduces the time needed to create a solvable 3-D-FEM model of the human femur to approximately 4 h. Mailer et al.‘*
76
view of different layers. According to F@zre systematically varied between -10” and 35”. of a right femur. The central points within to the maximum SED (compressive strain).
emphasized that the time needed to build up a detailed model of a specific specimen was reduced from 2-3 month (former manual method) to Z3 days using their subvolume technique. Our pre-
processing represents a further improvement in terms of modelling time in comparison to the method published by Miiller et ah” and is run Three different Young’s modnon-interactively. uli were assigned characterizing the cortical bone, the trabecular bone and the intramedullar tissue. Young’s moduli were homogenously modelled within these parts. However, using the empirical equations of Carter and Hayes” and Rice et aZ.‘“, it could be possible to assessthe density of bone at any particular location and to estimate its Young modulus. This specification was neglected at this stage of development, although the eighth macro of our preprocessor allows an element-specific assignment of material properties. This simplification does not significantly affect the primary objective of the study to identify the sensitivity of strain energy density (SED) patterns and the occurence of maxima or mimima. It is potentially significant, however, with respect to particular amounts of SED. Compared with the concept of automatic mesh generators (solid modelling), our approach has further advantages: The nodes are defined plane by plane in a regular and standardized way. Node numbers can be directly related to geometrical points. Usually, solid modelling preprocessing generates meshes with linear or parabolic tetrahedral elements. However, brick elements, which were selected can be expected to be of advantage in terms of computing time, accuracy and computer storage requirements. Keyak and coworkers-10.31 published an automated method of generating patient-specific models. A similar method was used by Hollister et aZ.““. It combines thresholding with an edge-following algorithm. The digital format of the data together with knowledge of the CT scan pixel size enables estimation of the bone dimensions. Their approach appears also to be a promising way for patient-specific modelling, although controversially discussedYg and will be compared with our geometry-based method in a future project in our group. Sawidis et ab”’ simulated load cases of hip joint flexion angles up to 60” and demonstrated FEMresults indicating a torsional deformation of the femur. Cheal et al.“” applied literature data of muscle forces in their femoral FEM-model and analyzed three phases of gait and four extreme loads. The load cases tested in the study presented here, especially the variations of the resultant force, represent not in all cases realistic values according to typical daily-life loadings. This kind of load application became necessary to assure consistency of the test& system and to restrict the analyses on the sensmvity of SED patterns. The tests on the influence of the model size provided information on how long a femur has to be modelled, in order to analyze a specific part of the bone accurately. It may be assumed that these results and the results of the various mesh densities will also be applicable when testing stems of hip endoprostheses of different lengths and help to minimize the computing costs. The variation of the Young’s moduli assigned showed that the intramedullar tissue is negligible
in a standard stress/strain analysis. In contrast, a considerable sensitivity of the results was seen when varying the cortical Young’s modulus. Hence, this parameter should be carefully measured whenever possible. The cancellous Young’s modulus influenced the results at the proximal and distal part of the femur. A final conclusion on the influence of the cancellous Young’s modulus, however, must not be made because inhomogeneity and anisotropy should be parameters of further tests. As regards the cortical bone, good agreement between experimental and theoretical stress analysis was reported when assuming a linear elastic, homogeneous and transversely isotropic behaviou?. The majority of investigators have frequently neglected the iliotibial tract force’.‘.-‘.i.“. However, the sensitivity results presented here emphasize a striking influence on the SED patterns. A decrease of the principal stresses when implementing an iliotibial tract force has been previously reported by Rohlmann et nl.“. The influence is not uniform because the minimum SED at the distal shaft was calculated when applying a small tract force. The minimum SED of the greater trochanter and the proximal shaft were computed after definition of The variation of the resultant a large tract force. force direction within the sagittal plane demonstrates a SED minimum at an anterior inclination of 13” and a striking increase in case of a smaller or larger angle. This value was approximately predicted already in 1935” as a physiological value in the single leg stance. Hence, it could be concluded that the morphology of the femur is reasonably adapted to this kind of sagittal force direction, and care should be taken to keep this direction when implanting a hip endoprosthesis. The variation of the resultant force direction within the frontal plane reveals a SED minimum at an angle of 21”. It has to be stressed that the SED was calculated at a relativelv low level between 5 and 25” and showed a marked increase at larger angles. Our interpretation of these results is that a small medial inclination of the action line does not very much affect the SED pattern and would not represent a primary risk in total hip arthroplasty. In contrast, large angles (>30”) should always be avoided. When experimentally or theoretically testing hip stems, load configurations of maximum SED values may be used. In case of a tolerable mechanical behaviour under these conditions, it may be concluded that the other load conditions would probably also be tolerated by the system. The orienlation of the connecting lines between the surface-SED-peaks in the horizontal view was found to be sensitive relative to the external load direction in the sagittal plane. Similar results were reported by Savvldis et nl.“’ and may be compared with morpholopcal findings of honzontal material distributions-“. This suggests that t.he relevance of a cross-sectional design adaptation of femoral stems to the pattern of surfaceSED-peaks connecting lines may be an important question of further tests.
77
Sensitivity
offemoral
strain pattern
analyses:M. Imgsfeld
et al.
ACKNOWLEDGEMENTS The study was supported by CADFEM GmbH, Grafing, Germany. The authors wish to thank Dr H. Sitter, Department of Theoretical Surgery, University of Marburg for his mathematical assistance.
19.
20. 21.
REFERENCES 1. Brekelmans WAM, Poort HW, Sloof TJH. A new method to analyse the mechanical behaviour of skeletal parts. Actu Orthop &and 1972; 43: 301-17. 2. Riihrle H, Scholten R, Sollbach W, Ritter G, Gtinert, A. Der Kraftfluss bei Hiiftendoprothesen. Arch orthop UnfallChir 1977; 89: 49-60. 3. Rohlmann A, Mossner U, Bergmann G, Kijlbel R. FiniteElement-Analysis and experimental investigation in a femur with hip endoprosthesis. J Biomech 1983; 16: 72742. 4. Huiskes R, Weinans H, Dalstra M. Adaptive bone remodelling and biomechanical design considerations. Orthopedics 1989; 12: 1255-67. 5. Walker PS, Poss R, Robertson DD, Reilly DT, Ewald FC, Thomas WH, Sledge CB. Design analysis of press-fit hip stems. Orthop Rel Sci 1990; 1: 75-85. und 6. Lengsfeld M, Weiss H, Kienapfel H. Morphologische numerische Untersuchungen zur Stresskompatibilimt von Patellaimplantaten. Biomed Z’echnik 1992; 37: 222-9. 7. Lengsfeld M. Stresses at the meniscofemoral joint. Elastostatic investigations on the applicability of interface elements. J Biomed Eng 1993; 15: 324-8 8. Rubin PJ, Rakotomanana RL, Leyvraz PF, Zysset PK, Curnier A, Heegaard JH. Frictional interface micromotions and anisotropic stress distribution in a femoral total hip component. J Biomech 1993; 26: 725-39. 9. Lengsfeld M, Kaminsky J, Merz B, Franke RP. Automatisierte Generierung von 3-D Finite Elemente Codes des menschlichen Femurs. Biomed Tecknik 1994; 39: 117-22 10. Mtiller R, Merz B, Rtiegsegger P. Interfacing quantitative computed tomography with finite element modeling. Spie Vol. 1395 Close-Range Photogrammetry Meets Machine Vision 1990; 1092-g. B. Die Beanspruchung des menschlichen 11. Tillmann Htiftgelenks. III. Die Form der Facies lunata. 2 Anat EntwickL-Gesch 1969; 128: 329-49. Ein mechanisches 12. Pauwels F. Der Schenkelhalsbruch. Problem. Beilageheft Z onhop Chir 1935; 63. 13. Pauwels F. Biomechanics of the Locomotor Apparatus. Berlin: Springer, 1980. 14. Hamacher P, Roesler H. Die Berechnung von Grosse und Richtung der Htiftgelenksresultierenden im Einzelfall. Arch arthop Unfall-Chir 1971; 70: 26-35. 15. Hamacher P, Roesler H. Ergebnisse der Berechnung von Grosse und Richtung der Hiiftgelenksresultierenden im Einzelfall. Arch orthop Unfall-Chir 1972; 72: 94-106. 16. Kummer B. Einfiihrung in die Biomechanik des Htiftgelenks. Heidelberg, New York: Springer 1985. 17. Merz B. Finite Element-Analyse von exzidierten Femora basierend auf der quantitativen Computer Tomographie. [Dissertation] Swiss Federal Institute of Technology Ziirich, Switzerland, 1993. 18. Heimkes B, Pose1 P, Plitz W, Jansson V. Forces acting on the juvenile hip joint in the one-legged stance. JPediatr Orthop 1993; 13: 431-6.
78
22. 23.
24.
25. 26.
27.
28. 29. 30. 31. 32.
33.
34.
35.
36.
Kummer B. Is the Pauwels’ theory of hip biomechanics still valid? A critical analysis, based on modern methods. Ann Anat 1993; 175: 203-10. Carter DR, Fyhrie DP, Whalen RT. Trabecular bone density and loading history: regulation of connective tissue biology by mechanical energy. J Biomech 1987; 20: 785-94. Huiskes R, Weinans H, Grootenboer HJ, Dalstra M, Fudala B, Sloof TJ. Adaptive bone-remodeling theory applied to prosthetic-design analysis. J Biomech 1987; 20: 1135-50. C&e, H. Werkstoffkundliche und mechanische Optimierung von zementfieien Htiftendo@othesen. Regensburg: Roderer 1992. Huiskes R, Janssen JD, Sloof TJ. A detailed comparison of experimental and theoretical stress analyses of a human femur. In: Cowin S.C. (Ed.) Mechanical Properties of Bone AMD 1981; 45: 211. The American Society of Mechanical Engineers, New York: American Sot. Medical Engineers. Keyak JH, Fourkas MG, Meagher JM, Skinner HB. Validation of an automated method of three-dimensional finite element modelling of bone. J Biomed Eng 1993; 15: 505-g. Beaupre GS, Orr TE, Carter DR. An approach for timedependent bone modeling and remodeling-theoretical development. J Orthop Res 1990; 8: 651-61. Van Rietbergen B, Huiskes R, Weinans H, Sumner DR, Turner TM, Galante JO. The mechanism of bone remodeling and resorption around press-fitted THA stems. J Biomech 1993; 26: 369-82. Merz B, Mtiller R, Rtiegsegger P, Niederer P. Stress analysis of excised human femora based on quantitative computed tomography. KVZ Meeting of the European Society of Biomechanics, Rome, 1992. Carter DR, Hayes WC. The compressive behaviour of bone as a two-phased porous structure. J. Bone Jt Surgery 1977; 59A: 954-62. Rice JC, Cowin SC, Burr DB. On the dependence of the elasticity and strength of cancellous bone on apparent density. J Biomech 1988; 16: 155-68. Keyak JH, Meagher JM, Skinner HB, Mote CD. Automated three-dimensional finite element modelling of bone: a new method. J Biomed Eng 1990; 12: 389-97. Keyak JH, Skinner HB. Three-dimensional finite element modelling of bone: effects of element size. f Biomed Eng 1992; 14: 483-9. Hollister SJ, Breean JM, Kikuchi N. A homogenization sampling procedure for calculating trabecular bone effective stiffness and tissue level stress. J Biomech 1994; 27: 43344. Marks LW, Gardner TN. The use of strain energy as a convergence criterion in the finite element modelling of bone and the effect of model geometry on stress convergence. J Biomed Eng 1993; 15: 474-6. Sawidis E, Lijer F, Grtiters H, Wiesener C. Analyse der Beanspruchung des proximalen Femur bei verschiedenen Arten der Belastung mit Hilfe der Finite-ElementMethode. Z Orthop 1991; 129: 268-77. Cheal EJ, Spector M, Hayes WC. Role of loads and prosthesis material properties on the mechanics of the proximal femur after total hip arthroplasty. J Orthop Res 1992; 10: 405-22. Breul R. Die Auswertung der Dichteverteilung im Riihrenknochenquerschnitt von menschlichen Femora fur die Bestimmung der Richtung der gunstigsten Biegebeanspruchung. Z Morph Anthrop 1985; 76: 63-76.
Printed in Great Britain. All rights reserved
Sensitivity of femoral strain pattern analyses to resultant and muscle forces at the hip joint M. Lengsfeld*,
J. Kaminsky”,
B. Merzt
and R. P. Franke’
*Department of Orthopaedic Surgery, Philipps-University of Marburg, Germany; tbrstitute of Biomedical Engineering and Medical Informatics, Swiss Federal Institute of Technology and University of Zurich, Switzerland; IInstitute of Biomedical Engineering, Department of Biomaterials, University of Ulm, Germany Received
14 December
1994, accepted
27 March
1995
AJWIRACT An automated geometrical preprocessor was dareloped with the aim of creating three-dimensional finite element models (IZM) of the human femur. On the basis of postprocessed computed tomography data, this preprocessor makes possible rapid, flexible and regular meshing with ‘brick’ elements. Three different material properties were modelled at the present stage of development. Sensitivity anabses demonstrated that the strain energy dasity (SELI) patterns of the dafferent femoral parts were most sensitive to the implementation of an iliotibial tract force. The variation of the resultant hip force and abductor force direction within the sagittal plane demonstrated a SED minimum at an anterim inclination of 17’; the variation of the resultant force direction within the frontal plane demonstrated a minimum SELI at a medial inclination of 21” relative to the mechanical axis of the lower limb. The orientation of the connecting line between the surface-SEL+eak.s in the horizontal uitw was found to be most sensitive to the variation of the resultant force within the sagittal plane.
Keywords: FEM, Med.
Eng.
Phys.,
preprocessing, 1996,
Vol.
human
femur,
force
The finite element method (FEM) is widely accepted as a powerful and also a rather sensitive numerical tool, which provides a way to simulate biomechanical phenomenala. FEM-software packages have powerful and convenient preprocessors, although not specifically adapted to the requirements of biomechanical modelling. An experienced investigator is required to develop models interactively as realistic as possible and to avoid errors concerning input model parameters and boundary conditions. In fact, irregular geometry, complex microstructure of biological tissues and loading situations are specific prob lems of the application of the FEM in biomechanits and are still difficult to model. Nevertheless, from our clinical point of view, it is desirable to provide specific biomechanical FEM-preprocessors, which should be easy to handle and which could allow a fast model development. This would also be a contribution to a faster investigation of individual variations and, in the future, could allow one to utilize the method in the clinical setting, preoperatively. to: Dr Markus PhilippsUniversity Germany.
sensitivity
analyses
18, 70-78, January
INTRODUCTION
Correspondence paedic Surgery, 35033 Marburg,
application,
Lengsfeld, Department of Ortbw of Marburg. Baldingerstrasse, D-
Therefore, our project has two objectives: (1) the introduction of a new automated FEM-preprocessor in order to create S-D-models of human femora, based upon postprocessed computed tomography data; and (2) a sensitivity study focusing on the influence of force transmission, meshing, model size and material properties on femoral strain energy density (SED) patterns. MATERIALS
AND METHOD
Description of the preprocessor In order to provide geometrical data, the inner and outer contours of the femoral cortical bone were extracted from computed tomograms (CT) with the method of thresholding. As cortical regions in cross-sections near the joints are very thin, modelling of the cortical bone by extracting the inner contour was neglected in these regions ‘,l”. Ten programs were developed, creating, step by step, a solvable 3D-FEM-model of the femur. The parametric design language of the ANSYS 5.0 general purpose FEM-package (Swanson Analysis Systems, Inc. Houston, PA, U.S.A) including ANSYS preprocessing commands was used to write these programs (macros) : 1. Definition
of the inner and outer femoral
sur-
Sasitiuity
2.
3.
4. 5.
6.
face by key-point-based interpolating spline functions. Description of the femoral curvature by a longitudinal mid-femoral line (vertical connecting line). This line represents the connection of the centres of the smallest circles fitting in the outer femoral contour (ZGg~re I). These centres including the centre of the femoral head were worked out slice after slice and were iteratively calculated (Eg~e 2). Input parameters were a tolerance value and a starting position vector within a transverse femoral slice. Definition of lines subdividing the transversal cross-sections into several sectors. The lines are orientated perpendicularly to the longitudinal mid-femoral line. If necessary, an optional program provides a smoothing of the outer femoral surface. Implementation of nodes along the sector lines and deletion of the lines afterwards. The centrally located intramedullary part is subdivided by nodes describing a rectangular geometry (Figure 1). Starting proximally, the predefined nodes
strain pattern
anal~srs:
M. Imgsfeld
et al.
were meshed with 3-D isoparametric 8-nodes brick elements (Figure 3). 7. Detection of distorted elements, which may cause numerical problems during the solution and tetrahedral-like reshaping of these elements without changing of the element 8. x!&ment of material properties to the different parts of the femur. 9. Definition of the boundary conditions. The nodes at the distal end of the model were constrained. The resultant hip force was applied by multiple force components at the nodes located within the geometrical extent of the load bearing surface l1 . The pelvi-trochanteric muscle force was applied at the tip of the greater trochanter in order to counterbalance the partial body weight in the single-leg stance phasel’-lt’~ The angle of the resultant force to the sagittal and frontal plane, the partial body weight and its lever arm are input parameters. 10. Start of the Solver. The standard model generated and tested here comprised 3308 elements, 3750 nodes and approximately 11,200 degrees of freedom. The thickness of the element layers was 12 mm. A finer mesh was modelled at the proximal femur (4 mm element layer). The calculations were run on an i486 DX 2-66 local bus personal computer and it took 3 h 48 min of computing time to generate the model and 1 h 22 min to run the strain analysis (solution). Input
Figure 1 The femoral curvature is described by a longitudinal midfemoral line. This line represents the connection of the centres of the smallest circles fitting in the outer femoral contour. Sector lines are orientated perpendicularly to the longitudinal mid-femoral line. The cross-section of the trochanteric region after deletion of the mid-femoral and sector lines demonstrates the architecture of the final mesh. A: anterior; L: lateral; M: medial; P: posteriox-
offpmoral
parameters
and description
of the model
Mesh density of the layers of brick elements is determined by the number of lines subdividing the layer into several sectors. The total number of element layers can be defined by changing their thicknesses. The thickness of the layers can be independendy selected from the table-feed of the computed tomograms. Finally, the user selects the Poisson number and three different Young’s moduli for the cortical and trabecular bone and the intramedullar tissue, which were previously separated by contour detection. Young’s moduli are averaged in case of elements located on the border between two differently stiff volumes by the eighth macro. If there is no further information given, the Young’s moduli assigned were 17,000 MPa (cortical bone), 1000 MPa (trabecular bone), 10 MPa (intramedullar tissue) and isotropic mechanical behaviour was assumed in the present study. A Poisson number of 0.33 was selected. About 85% of the femoral length was modelled, and the most distal layer of the model was located above the femoral condyles. The nodes within this plane were constrained in all directions. The model was based on CT scan slices taken from a right fresh-frozen femoral specimen. It was obtained from a male adult of 45 years with a weight of 77 kg”. An antero-posterior X-ray excluded pathological lesions. In case of an integration of the iliotibial tract
71
Sensitivity
offemoral
strain pattern
analyses: M. Lerpfeld
et al.
Starting position vector (C) within a horizontal femoral slice. Predefined tolerance value (a) of the centre to be determined.
1 Minimum distance vector (Dmin) between C and the outer femoral contour. The distance to the opposite outer femoral contour is expressed by S*Dmin.
j ! s/2
H s
Minimum distance (D) of the position vector : C+
s
* Dmin
Final centre : C+S*Dmin
H C
C
Figure 2 Flow chart of the iteration. A mid-femoral line is calculated, describing the 3D curvature of the femur. a: predefined tolerance value; convergence was based on Is*Dminl < a; C: starting position vector and position vector describing the centre of a transverse femoral slice, respectively; Dmin: minimum distance vector between C and the outer femoral contour; s: scalar; the product of the scalar s and the vector Dmin describes the distance vector between C and the opposite outer femoral contour; this vector has a direction opposite to that of Dmin
force in the model, it is applied at the middle of the greater trochanter and acts along the femoral axis according to Rohlmann et aL3 and Heimkes et al . I* . The standard load configuration was a resultant force of 2624 N, a gluteus medius muscle (abductor) force of 2036 N, a partial body weight
72
of 615 N and an iliotibial tract force of 611 N. According to Kummerlg, the latter represented 30% of the gluteus medius muscle force. The lever arm of the partial body weight was 108 mm. Forces, force components and force directions were described in this study relative to an orthog-
head
-__----
--
n&
______---
--
trochaater 1 -----trochanter 2 ------
shaft
2
-------
---
shaft 3 __-_---- --
IL
P-A
x
Jr
Y A
X
M-L
A-P
&aft 5 ____--______
L Y L-M
Figure 3 Element plots of the model. The standard load configuration of a resulrant force of 2624 N, a gluteus medius muscle (abductor) force of 2036 N, a partial body weight of 615 N and an iliotibial tract force of 611 N is displayed. A: anterior; IL lateral: M: medial; P: posterior
onal coordinate system, where its origin coincided with the centre of the femoral ‘head. The vertical axis was defined by the centre of the femoral head and the estimated knee joint centre within the femoral groove. The direction of the frontal axis was defined by the dorsal aspect of the femoral condyles. Relative to this coordinate system and according to the moment equilibrium of the single leg stance, the resultant force was 21” mediolaterally and 13” frontodorsally directed. The femoral shaft axis determined the iliotibial tract force orientation. A neck-shaft angle of 136” and an anteversion angle of 19” were measured at the femur modelled in this study. The results of the calculations testing the influence of force applications on the whole model were approximated (smoothed) by a non-interpolating third-order-polynomial using Excel (Microsoft Corporation, U.S.A.). The FEM-results were expressed by the SED (SED = strain energy/bone volume). It was hypothesized that the SED is suitable for comparative strain analysis of bone20-2’. The femur was divided into ten different parts in order to postprocess the results also relative to particular femoral parts: the head, neck, proximal trochanter and distal trochanter. The shaft was subdivided into six slices (F&WY 4).
shaft 6 -- - -- --
Figure 4 The femur was divided into ten dif’ferent parts in order to postprocess the results also relative to particular femoral parts: the head, neck, proximal trochanteland distal trochantrr. The shaft was \uhdividrd into six slices
RESULTS Mesh
density
In order to test the influence of the mesh density on the SED distribution, the slice thickness within the shaft was modelled with 4, 8, 12 and 16 mm. The slice thickness within the head, neck and trochanteric region was 4 mm. Consequently, the number of elements was observed to be between approximately 2500 and 7000. Almost identical results were found when testing the models with 4, 8, 12 and 16 mm slice thickness. Model
size
Using the standard boundary conditions described above, the model length, i.e. the height of the most distal layer, where all nodes were constrained, was varied. Small influence was found of the height of the constrained layer on the SED in the proximal femoral arts. The complete model is represented by 100 7o and the tests were performed assuming a model length between 100 and 50%. It was observed that only the two most distal layers were not constant. This distal part represented about 10% relative to the complete model.
73
Sensitivity
ojjaoral
Young’s
strain pattern
analyses: M. Lagsjeld
et al.
moduli
Variation of the resultant force and abductor force direction within the sag&al plane
Studying the SED of the bone when varying the Young’s modulus of the intramedullar tissue between 5 and 1000 MPa, there was no significant difference seen. Increasing the Young’s modulus of the cortical bone from 1000 to 20,000 MPa, a decrease of the SED was found and was particularly significant at the distal shaft. The SED ranged from 12 to 1 MPa/103. The variation of the cancellous Young’s modulus between 10 and 10,000 MPa showed a decreasing pattern of the SED at the head (2.5-0.4 MPa/lO ), neck (2.10.6 MPa/lO’) and the distal trochanter (3.22.2 MPa/103) and at the distal layers of the shaft (1.3-0.9 MPa/103). The other parts of the shaft were almost unaffected. Iliotibial
tract
The direction of the resultant force and the abductor force within the sagittal plane were systematically varied between -10” and 35” with an angle step of 5”, except between 0” and 20”, where an angle step of 1” was selected. The amount of the forces and their frontal components are constant. The negative sign refers to a dorso-frontal force direction (posterior inclination), and the positive sign to a fronto-dorsal direction (anterior inclination). In order to present the results (Figure I$‘, the SED of all elements that the model consisted of was averaged and the SED of the elements of various parts was averaged. The averaged SED values of the whole model were approximated by a third-order-polynoma (y:SED[MPa/103]; xangle[“]): y = 0.0002X7 + 0.01152 - 0.395x + 5.864 having a determination coefficient (72) of 0.93. The SED minimum of the polynomial was found at an anterior inclination of 13” and increased significantly at -10 and 35” inclination.
force
The influence of the iliotibial tract force on the distribution of the SED was tested by varying this force between 0 and 100% of the abductor force with a step of 5%. The SED values of the femoral head and neck were not altered by the introduction of the iliotibial tract force. The minimum SED of the distal femoral shaft (shaft 5 and 6) was observed between 0 and 15%. In contrast, the minimum SED of the greater trochanter and the proximal shaft were seen at 100% (trochanter 1 and 2, shaft 1). The minimum SED of the whole model was calculated at 45% (Figure 5).
Variation of the resultant the frontal plane
force
direction
within
The medio-lateral direction of the resultant force within the frontal plane was systematically varied between 7” and 42” with an angle step of 5”, except between 12” and 32’, where an angle step
12 +
head
-.x-neck -
nochanter I
-
trochsnter 2
+
shaft 2
--e-
shaft 3
tshsft4
+ -whole
shaft 6 model
0 0
5
10
15
20
25
30
35
40
45
50
55
60
65
70
75
80
85
90
95
100
Iliotibial tract force [ % of the abductor force ] Figure 5 The influence of the iliotibial tract force [N] on the distribution of the strain engery density (SED [MPa/10”]) was tested by varying this force between 0 and 100% of the abductor force with a step of 5%. The SED values of the femoral head and neck were not altered by the introduction of the iliotibial tract force. The minimum SED of the distal femoral shaft (shaft 5 and 6) was observed between 0 and 15%. In contrast, the minimum SED of the greater trochanter and the proximal shaft were seen at 100% (trochanter 1 and 2, shaft 1)
74
Sensitivity
oj,fpmmnl
@rain patk-rn
ana1y.sP.y: M. i.mgsfild
et al.
r
-+--head +-neck ---.
trochanter
I
--
rrochanter
2
-*-shaft2 -+-shaft
+
3
shaft 4
tshaff6 -whole -polynomial (whole 0
5
15
IO
20
Angle in YZ - plane of th; bone-coordinate
25
30
model
model)
35
system 1 O ]
Figure 6 The direction of the resultant force and the abductor force within the sagittal (yr) plane were systematically varied between -10 and 35”. The negative sign refers to a dorso-frontal force direction (posterior inclination), the positive sign to a fronto-dorsal direr&ion (anrerior inclination). The strain energy density (SED) of all elements which the model consisted of was averaged and the SED of the elements of various parts was averaged. The averaged SED values of the whole model were approximated by a third-order-polynomial (y:SED [MPailO’]; rangie[‘J)
of 1” was selected. An angle of 0” refers to a vertical direction of the resultant force and could not be tested because of a missing equilibrium. The direction of the resultant and the partial body weight were given. In order to match the equilibrium condition, the amount of the resultant force and the amount and direction of the abductor force had to be recalculated for every angle step, as described by Hamacher and Rbsler’ *I’. The resultant force direction within the sag&al plane was 13’. The averaged SED values of the whole model were approximated by a third-order-poly(y:SED[MPa/IO”]; xangle[“]): y = nomial 0.0004$ + 0.00552 - 0.0933x + 4.4487 having a determination coefficient (?) of 0.88. The minimum SED of this polynomial was found at an angle of 21°, and increased significantly at 36” (Fiflrf? 7). Orientation of the connecting lie between surfaceSED-peaks in the transverse view
the
When looking at transverse slices of the model, the medial and lateral elements of maximum SED were selected. At the proximal shaft the connecting line between the SED peaks had a slight antero-medial/postero-lateral orientation, while almost no deviation from the frontal direction was seen in the central &aft. A particular sensitivity of this parameter was found during the variation of the hip load within the sagittal plane. A large anterior inclination caused an internal rotation of the connect.ing line (&&Lw 8).
DISCUSSION Assuming a simplified experimental setup of femoral load application the finite element method has been validated by strain gauge measurements3,“.23,24 and frequently utilized. Moreover, the method has been applied to model the response of bone to a previously calculated stress field, which means a modification of the material property assignments during the remodelling simulation20,?1,2~,2~.
From our point of view, however, there are still two important aspects restricting the applicability and validity of this testing method: (a) generation of the models appears still rather time-consuming and requires engineering expertise making the introduction in the clinical setting and as a clinical screening tool in viuo almost impossible; and (b) the majority of researchers used a proximal load application according to the classical model of Pauwels’* of the single leg stance phase’.“,“. There is still not very much known about the sensitivity of the amount and direction of maximum stresses on the proximal load direction modelling. This and other input parameters have to be carefully tested before finite element analysis can also be safely and reliably utilized by the clinician. In addition to the use of FEM in the future in the clinic, it would still be desirable to learn more about the behaviour and sensitivity of femoral FEM-models before drawing clinical or design consequences from FEM-tests of implants. Based on a contour data file extracted from
Sensitivity
of femoral
strain pat&em analyses: M. Lengsjeld et al.
25 -+-head
-
ttochanter
1
-
trochanter
2
+
shaft 2
-+
shaft 3
+
shaft 4
-o-shaft -a-whole
Od....:....:....:. 5
, 10
15
.
20
.:.
, . .:. 2s
. . 30
Angle in XZ - plane of thd bone-coordinate
.:-, 35
5 shafl6 model
. .:.. 40
system [ O 1
Figure 7 The medio-lateral direction of the resultant force within the frontal (xz) plane was systematically varied between 7” and 42”. The direction of the resultant and the partial body weight were given. In order to match the equilibrium condition, the amount of the resultant force and the amount and direction of the abductor force had to be recalculated for every angle step. The resultant force direction within the sagittal plane was 13”. The averaged SED values of the whole model were approximated by a third-order-polynomial (y:SED [MPa/103]; xangler])
Angle in YZ - Diane of the bone-coordinate
-
system [ o ]
8 b Q 8 Q Q
head
1 t
neck trochanter
1
shaft 4 shaft 5
Q
@
0
Q
Qj
8
QJC)
shaft 6
Q
Q
Q
@
Q
d
a-
-A--l@ .---_-
whole model Figure 8 Orientation of the surfaceSED (strain energy density)-peaks in a schematic transverse 6, the direction of the resultant force and the abductor force within the sagittal (yz) plane were The left side of the circles refers to the medial, and the top of the circles refers to the ventral part the circumference of the dark half refer to the maximum SED (tensile strain) and of the white half
computed tomography, the finite element preprocessor presented here reduces the time needed to create a solvable 3-D-FEM model of the human femur to approximately 4 h. Mailer et al.‘*
76
view of different layers. According to F@zre systematically varied between -10” and 35”. of a right femur. The central points within to the maximum SED (compressive strain).
emphasized that the time needed to build up a detailed model of a specific specimen was reduced from 2-3 month (former manual method) to Z3 days using their subvolume technique. Our pre-
processing represents a further improvement in terms of modelling time in comparison to the method published by Miiller et ah” and is run Three different Young’s modnon-interactively. uli were assigned characterizing the cortical bone, the trabecular bone and the intramedullar tissue. Young’s moduli were homogenously modelled within these parts. However, using the empirical equations of Carter and Hayes” and Rice et aZ.‘“, it could be possible to assessthe density of bone at any particular location and to estimate its Young modulus. This specification was neglected at this stage of development, although the eighth macro of our preprocessor allows an element-specific assignment of material properties. This simplification does not significantly affect the primary objective of the study to identify the sensitivity of strain energy density (SED) patterns and the occurence of maxima or mimima. It is potentially significant, however, with respect to particular amounts of SED. Compared with the concept of automatic mesh generators (solid modelling), our approach has further advantages: The nodes are defined plane by plane in a regular and standardized way. Node numbers can be directly related to geometrical points. Usually, solid modelling preprocessing generates meshes with linear or parabolic tetrahedral elements. However, brick elements, which were selected can be expected to be of advantage in terms of computing time, accuracy and computer storage requirements. Keyak and coworkers-10.31 published an automated method of generating patient-specific models. A similar method was used by Hollister et aZ.““. It combines thresholding with an edge-following algorithm. The digital format of the data together with knowledge of the CT scan pixel size enables estimation of the bone dimensions. Their approach appears also to be a promising way for patient-specific modelling, although controversially discussedYg and will be compared with our geometry-based method in a future project in our group. Sawidis et ab”’ simulated load cases of hip joint flexion angles up to 60” and demonstrated FEMresults indicating a torsional deformation of the femur. Cheal et al.“” applied literature data of muscle forces in their femoral FEM-model and analyzed three phases of gait and four extreme loads. The load cases tested in the study presented here, especially the variations of the resultant force, represent not in all cases realistic values according to typical daily-life loadings. This kind of load application became necessary to assure consistency of the test& system and to restrict the analyses on the sensmvity of SED patterns. The tests on the influence of the model size provided information on how long a femur has to be modelled, in order to analyze a specific part of the bone accurately. It may be assumed that these results and the results of the various mesh densities will also be applicable when testing stems of hip endoprostheses of different lengths and help to minimize the computing costs. The variation of the Young’s moduli assigned showed that the intramedullar tissue is negligible
in a standard stress/strain analysis. In contrast, a considerable sensitivity of the results was seen when varying the cortical Young’s modulus. Hence, this parameter should be carefully measured whenever possible. The cancellous Young’s modulus influenced the results at the proximal and distal part of the femur. A final conclusion on the influence of the cancellous Young’s modulus, however, must not be made because inhomogeneity and anisotropy should be parameters of further tests. As regards the cortical bone, good agreement between experimental and theoretical stress analysis was reported when assuming a linear elastic, homogeneous and transversely isotropic behaviou?. The majority of investigators have frequently neglected the iliotibial tract force’.‘.-‘.i.“. However, the sensitivity results presented here emphasize a striking influence on the SED patterns. A decrease of the principal stresses when implementing an iliotibial tract force has been previously reported by Rohlmann et nl.“. The influence is not uniform because the minimum SED at the distal shaft was calculated when applying a small tract force. The minimum SED of the greater trochanter and the proximal shaft were computed after definition of The variation of the resultant a large tract force. force direction within the sagittal plane demonstrates a SED minimum at an anterior inclination of 13” and a striking increase in case of a smaller or larger angle. This value was approximately predicted already in 1935” as a physiological value in the single leg stance. Hence, it could be concluded that the morphology of the femur is reasonably adapted to this kind of sagittal force direction, and care should be taken to keep this direction when implanting a hip endoprosthesis. The variation of the resultant force direction within the frontal plane reveals a SED minimum at an angle of 21”. It has to be stressed that the SED was calculated at a relativelv low level between 5 and 25” and showed a marked increase at larger angles. Our interpretation of these results is that a small medial inclination of the action line does not very much affect the SED pattern and would not represent a primary risk in total hip arthroplasty. In contrast, large angles (>30”) should always be avoided. When experimentally or theoretically testing hip stems, load configurations of maximum SED values may be used. In case of a tolerable mechanical behaviour under these conditions, it may be concluded that the other load conditions would probably also be tolerated by the system. The orienlation of the connecting lines between the surface-SED-peaks in the horizontal view was found to be sensitive relative to the external load direction in the sagittal plane. Similar results were reported by Savvldis et nl.“’ and may be compared with morpholopcal findings of honzontal material distributions-“. This suggests that t.he relevance of a cross-sectional design adaptation of femoral stems to the pattern of surfaceSED-peaks connecting lines may be an important question of further tests.
77
Sensitivity
offemoral
strain pattern
analyses:M. Imgsfeld
et al.
ACKNOWLEDGEMENTS The study was supported by CADFEM GmbH, Grafing, Germany. The authors wish to thank Dr H. Sitter, Department of Theoretical Surgery, University of Marburg for his mathematical assistance.
19.
20. 21.
REFERENCES 1. Brekelmans WAM, Poort HW, Sloof TJH. A new method to analyse the mechanical behaviour of skeletal parts. Actu Orthop &and 1972; 43: 301-17. 2. Riihrle H, Scholten R, Sollbach W, Ritter G, Gtinert, A. Der Kraftfluss bei Hiiftendoprothesen. Arch orthop UnfallChir 1977; 89: 49-60. 3. Rohlmann A, Mossner U, Bergmann G, Kijlbel R. FiniteElement-Analysis and experimental investigation in a femur with hip endoprosthesis. J Biomech 1983; 16: 72742. 4. Huiskes R, Weinans H, Dalstra M. Adaptive bone remodelling and biomechanical design considerations. Orthopedics 1989; 12: 1255-67. 5. Walker PS, Poss R, Robertson DD, Reilly DT, Ewald FC, Thomas WH, Sledge CB. Design analysis of press-fit hip stems. Orthop Rel Sci 1990; 1: 75-85. und 6. Lengsfeld M, Weiss H, Kienapfel H. Morphologische numerische Untersuchungen zur Stresskompatibilimt von Patellaimplantaten. Biomed Z’echnik 1992; 37: 222-9. 7. Lengsfeld M. Stresses at the meniscofemoral joint. Elastostatic investigations on the applicability of interface elements. J Biomed Eng 1993; 15: 324-8 8. Rubin PJ, Rakotomanana RL, Leyvraz PF, Zysset PK, Curnier A, Heegaard JH. Frictional interface micromotions and anisotropic stress distribution in a femoral total hip component. J Biomech 1993; 26: 725-39. 9. Lengsfeld M, Kaminsky J, Merz B, Franke RP. Automatisierte Generierung von 3-D Finite Elemente Codes des menschlichen Femurs. Biomed Tecknik 1994; 39: 117-22 10. Mtiller R, Merz B, Rtiegsegger P. Interfacing quantitative computed tomography with finite element modeling. Spie Vol. 1395 Close-Range Photogrammetry Meets Machine Vision 1990; 1092-g. B. Die Beanspruchung des menschlichen 11. Tillmann Htiftgelenks. III. Die Form der Facies lunata. 2 Anat EntwickL-Gesch 1969; 128: 329-49. Ein mechanisches 12. Pauwels F. Der Schenkelhalsbruch. Problem. Beilageheft Z onhop Chir 1935; 63. 13. Pauwels F. Biomechanics of the Locomotor Apparatus. Berlin: Springer, 1980. 14. Hamacher P, Roesler H. Die Berechnung von Grosse und Richtung der Htiftgelenksresultierenden im Einzelfall. Arch arthop Unfall-Chir 1971; 70: 26-35. 15. Hamacher P, Roesler H. Ergebnisse der Berechnung von Grosse und Richtung der Hiiftgelenksresultierenden im Einzelfall. Arch orthop Unfall-Chir 1972; 72: 94-106. 16. Kummer B. Einfiihrung in die Biomechanik des Htiftgelenks. Heidelberg, New York: Springer 1985. 17. Merz B. Finite Element-Analyse von exzidierten Femora basierend auf der quantitativen Computer Tomographie. [Dissertation] Swiss Federal Institute of Technology Ziirich, Switzerland, 1993. 18. Heimkes B, Pose1 P, Plitz W, Jansson V. Forces acting on the juvenile hip joint in the one-legged stance. JPediatr Orthop 1993; 13: 431-6.
78
22. 23.
24.
25. 26.
27.
28. 29. 30. 31. 32.
33.
34.
35.
36.
Kummer B. Is the Pauwels’ theory of hip biomechanics still valid? A critical analysis, based on modern methods. Ann Anat 1993; 175: 203-10. Carter DR, Fyhrie DP, Whalen RT. Trabecular bone density and loading history: regulation of connective tissue biology by mechanical energy. J Biomech 1987; 20: 785-94. Huiskes R, Weinans H, Grootenboer HJ, Dalstra M, Fudala B, Sloof TJ. Adaptive bone-remodeling theory applied to prosthetic-design analysis. J Biomech 1987; 20: 1135-50. C&e, H. Werkstoffkundliche und mechanische Optimierung von zementfieien Htiftendo@othesen. Regensburg: Roderer 1992. Huiskes R, Janssen JD, Sloof TJ. A detailed comparison of experimental and theoretical stress analyses of a human femur. In: Cowin S.C. (Ed.) Mechanical Properties of Bone AMD 1981; 45: 211. The American Society of Mechanical Engineers, New York: American Sot. Medical Engineers. Keyak JH, Fourkas MG, Meagher JM, Skinner HB. Validation of an automated method of three-dimensional finite element modelling of bone. J Biomed Eng 1993; 15: 505-g. Beaupre GS, Orr TE, Carter DR. An approach for timedependent bone modeling and remodeling-theoretical development. J Orthop Res 1990; 8: 651-61. Van Rietbergen B, Huiskes R, Weinans H, Sumner DR, Turner TM, Galante JO. The mechanism of bone remodeling and resorption around press-fitted THA stems. J Biomech 1993; 26: 369-82. Merz B, Mtiller R, Rtiegsegger P, Niederer P. Stress analysis of excised human femora based on quantitative computed tomography. KVZ Meeting of the European Society of Biomechanics, Rome, 1992. Carter DR, Hayes WC. The compressive behaviour of bone as a two-phased porous structure. J. Bone Jt Surgery 1977; 59A: 954-62. Rice JC, Cowin SC, Burr DB. On the dependence of the elasticity and strength of cancellous bone on apparent density. J Biomech 1988; 16: 155-68. Keyak JH, Meagher JM, Skinner HB, Mote CD. Automated three-dimensional finite element modelling of bone: a new method. J Biomed Eng 1990; 12: 389-97. Keyak JH, Skinner HB. Three-dimensional finite element modelling of bone: effects of element size. f Biomed Eng 1992; 14: 483-9. Hollister SJ, Breean JM, Kikuchi N. A homogenization sampling procedure for calculating trabecular bone effective stiffness and tissue level stress. J Biomech 1994; 27: 43344. Marks LW, Gardner TN. The use of strain energy as a convergence criterion in the finite element modelling of bone and the effect of model geometry on stress convergence. J Biomed Eng 1993; 15: 474-6. Sawidis E, Lijer F, Grtiters H, Wiesener C. Analyse der Beanspruchung des proximalen Femur bei verschiedenen Arten der Belastung mit Hilfe der Finite-ElementMethode. Z Orthop 1991; 129: 268-77. Cheal EJ, Spector M, Hayes WC. Role of loads and prosthesis material properties on the mechanics of the proximal femur after total hip arthroplasty. J Orthop Res 1992; 10: 405-22. Breul R. Die Auswertung der Dichteverteilung im Riihrenknochenquerschnitt von menschlichen Femora fur die Bestimmung der Richtung der gunstigsten Biegebeanspruchung. Z Morph Anthrop 1985; 76: 63-76.