A novel modelling and simulation method of hip joint surface contact stress

ABSTRACT

Understanding the hip joint surface contact stress distribution characteristics is helpful to determine hip joint biomechanical features and abnormal pathological behavior. Firstly, a 3-dimensional static hip joint biomechanical model is built using analytical method of model in order to study biomechanical properties including bearing area, stress distribution and the peak value of the contact stress of the femoral head, which reveals the relationship between the biomechanical properties and its geometric parameters. Secondly, based on the finite element analysis of the hip joint model, the contact stress distribution on the surface of femoral head is acquired under the condition of the different joint force and the acetabulum coverage rate. Finally, according to the evaluation of the femoral head surface stress and contact stress peak under different load distribution, accuracy and universality of the biomechanical model is verified.

KEYWORDS:

• biomechanical model
• finite element analysis
• hip joint
• surface contact stress

Introduction

The hip joint is the body's largest weight-bearing joint. The ball-and-socket structure plays a role of conducting gravity. Among the multiple active joints in the body, the hip joint is the most stable and allows a wide range of human action to be completed. In addition, it is considered the most common multiple lesions joint in clinical studies. For instance, in arthritis caused by hip dysplasia, if the pathological change cannot be corrected in a timely manner, the acetabulum will divorce from the original position, which may cause dislocation. As a result, the acetabulum cannot cover the femoral head in the original ball-and-socket manner. Thus, the weight-bearing area on the articular surface is reduced, and the stress is increased. In addition, the stress is focused on the outer edge of the acetabulum. Eventually, hip osteoarthritis will occur. Therefore, establishment of a simulation model for contact stress on the surface of the hip joint will aid in the timely detection of the condition of the diseased joint. This will provide important reference data for the determination of treatment programs.

Research methods for analyzing hip biomechanical properties mainly include the experimental methods and various numerical methods.^1-4 For a long time, the biomechanical properties of human bones were analyzed by in vitro research. For example, Taylor et al. pointed out that, under the influence of the hip joint or gluteal muscles, there were significant difference between the femoral interior compressive stress and exterior tensile stress;⁵ Zhang et al. tested fresh femur specimens of healthy adults using an MTS testing machine, and found that, under longitudinal compression, 3-point bending, and axial torsion, the mechanical distribution borne by the femoral shaft was different in different sections under the same load.⁶ This method of determining bone biomechanical properties through mechanical experimentation was limited to obtaining the distribution model for contact stress in individual living tissue. Most researchers used numerical analysis simulation methods to simulate the biomechanical behavior of living tissues, of which the finite element method was the most widely used. It was initially applied in studies of mechanical problems in the cardiovascular system, and began to be used in biomechanical studies on orthopedics in the 1970s;⁷ in the 1980s, its application spread to various parts of the skeleton.⁸ For example, researchers at Pennsylvania State University established a simulation model of the legs using the finite element method.⁹ A fast finite element model was established by the Human-Computer Interaction Technology Laboratory in Washington University.¹⁰ These finite element analyses were developed to simulate the biomechanical properties of human tissues, with the advantages of complex modeling and higher accuracy in the results. However, it offered no advantage in terms of computing speed and the complex analysis process was simultaneously more demanding for the operator.

In summary, this study intended to solve 2 problems. Firstly, we were dedicated to finding a simple and special calculation method for contact stress in the articular surface of the hip joint. Secondly, we were dedicated to validating the accuracy of this calculation method in characterizing the contact stress in the articular surface of the hip joint. Based on the above objectives, this study established a precise 3-dimensional static mathematical model for contact stress in the articular surface of the hip joint. The model was used to calculate characteristics, such as the bearing area, stress distribution, and the peak value of contact stress, in the articular surface of the hip joint. The model yielded an innovative idea of the relationship between the mechanical properties and the geometric parameters of the hip surface. We further proposed a way of using the finite element method to establish a virtual model of the hip joint for analysis of contact stress distribution in the articular surface of the hip joint, as well as the peak value of the contact stress under different loads. The results of the analysis were compared to the constructed results of the mathematical model for verification of the accuracy and universality of the mathematical model for contact stress established in this study.

Materials and methods

Establishment of the mathematical model for contact stress in the articular surface of the hip joint

The hip spherical coordinate system was established based on the right hip as the research object. The center of the spherical joint was placed at the origin of the coordinate system. The X and Z-axes were in the coronal plane of the body, and Y-axis was in the sagittal plane of the body. The direction of the X-axis was defined from the right pointing to the left of the body; the Y-axis was defined as pointing to the back; and the Z-axis was defined as pointing vertically. Figure 1 shows the coordinate system of the established model of the hip joint.

Figure 1. Coordinate system of hip joint model.

The tissue fluid between the bone and cartilage could play a role of lubrication, thus in the model the influence of friction was ignored, and only the stress was considered. At this moment, the tangential stress acting on the spherical joint sphere was zero. The resultant force acting between the head and acetabulum was the radial pressure R, which was located on the contact surface of the hip joint, pointing toward the center of the sphere. Suppose the spheroidal coordinates of any point on the articular surface were (${\displaystyle \phi }$, ${\displaystyle \theta }$, r), and the spheroidal coordinates of the acting point of the resultant force R were (${\displaystyle \beta }$, ${\displaystyle \alpha }$, r), wherein, ${\displaystyle \theta }$ and ${\displaystyle \alpha }$ were the azimuthal angles; ${\displaystyle \phi }$ and ${\displaystyle \beta }$ were the polar angles; and r was the radius of the sphere. The resultant force R acting on the contact surface and the contact pressure p of any point on the spheroidal joint should meet the following integral relationship:(1) ${\displaystyle {\displaystyle {\int }_S^{}p\text{d}}S=R}$(1) Where in, $d\text{S}=r^2(\text{s}\text{i}\text{n}\theta \text{c}\text{o}\text{s}\phi ,\text{s}\text{i}\text{n}\theta \text{s}\text{i}\text{n}\phi ,\text{c}\text{o}\text{s}\theta )\text{s}\text{i}\text{n}\theta \text{d}\theta \text{d}\phi$, $R=R\left(\text{sin}\alpha \text{cos}\beta ,\text{sin}\alpha \text{sin}\beta ,\text{cos}\alpha \right).$

The elastic modulus of the articular cartilage was significantly smaller than that of the femur and acetabulum. The femoral head was considered as a rigid sphere, and the acetabulum was as a rigid hemispherical shell. These two components were concentric to each other, and the articular cartilage was located centrally between them. When a load was applied to the hip joint, the articular cartilage was extruded, becoming thinner. However, the deformation was within the elastic range. Under the physiological load, the radial strain of the articular cartilage was significantly smaller than the radius of the femoral head. The articular cartilage was considered as a linear elastic body. Suppose the line of direction through the origin of the coordinate system, and parallel to the load direction, was intersected by the acetabular hemispherical surface. According to the principles of theoretical mechanics, this point of intersection coincided with the acting point of the resultant force. At this moment, the maximum stress ${\displaystyle p_{\text{max}}}$ on the spherical surface of the joint was the contact stress ${\displaystyle p_R}$ on the acting point of the resultant force. Therefore, the stress of contact pressure p of any point on the articular surface and the contact pressure ${\displaystyle p_R}$ on the acting point of the resultant force should meet the following linear relationship:(2) ${\displaystyle p=p_R\text{cos}\gamma }$(2) Where in, ${\displaystyle \gamma }$ is the central angle between any point on the contact surface and the acting point of the resultant force.

Before discussing the bearing area of the femoral head surface, the concept of central-edge angle (CE angle) should be introduced. It is the included angle between the vertical of the center of the femoral head and the outer edge of the acetabulum. The CE angle was an important indicator to measure the contact area of the hip joint. The contact area decreased with the reduction of the CE angle. Suppose the acetabulum was a regular hemispherical concave, then the contact area and the CE angle have a linear relationship. The CE angle value can be gotten by using angle measuring function in Mimics software. The result of angle measurement is 51.04°. For building acetabulum and cartilage model of different CE angle, other angle should be marked, as shown in Fig. 2.

Figure 2. CE angle measurement.

While standing motionless on a single foot, the component of contact force of the hip joint perpendicular to the coronal plane was only 3% of the resultant force in the coronal plane. Thus, this could be ignored.¹¹ Therefore, the contact force of the hip joint only included 2 components, namely ${\displaystyle R=(R\text{sin}\alpha ,0,R\text{cos}\alpha )}$. The angles ${\displaystyle \alpha }$, ${\displaystyle \beta }$, ${\displaystyle \theta }$, ${\displaystyle \phi }$, and θ_CE measured counterclockwise around the axis were stipulated as positive when they orientated forward along the axis from the origin of the coordinate system. The coordinate axis was rotated so that the acting point of the resultant force was located on the vertex of the spherical joint. The integral relations of the mathematical model in the newly established coordinate system were as follows:(3) ${\displaystyle \{\begin{array}{l}\text{d}\text{S}=r^2(\text{cos}\theta \text{sin}\phi ,\text{sin}\theta ,\text{cos}\theta \text{cos}\phi )\text{cos}\theta \text{d}\theta \text{d}\phi \\ R_x=p_R{r}^2{\displaystyle \int {\text{cos}}^3\theta \text{d}\theta {\displaystyle \int \text{cos}\phi \text{sin}\phi \text{d}\phi =R\text{sin}\alpha }}\\ R_y=p_R{r}^2{\displaystyle \int {\text{cos}}^2\theta \text{sin}\theta \text{d}\theta }{\displaystyle \int \text{cos}\phi \text{sin}\phi \text{d}\phi }=0\\ R_z=p_R{r}^2{\displaystyle \int {\text{cos}}^3\theta \text{d}\theta }{\displaystyle \int {\text{cos}}^2\phi \text{d}\phi =R\text{cos}\beta }\end{array}}$(3) Where $\theta \in \left[-\frac{\pi }{2},{\theta }_{\text{CE}}-\alpha \right]$, $\phi \in \left[-\frac{\pi }{2},\frac{\pi }{2}\right]$.

These integrations were solved, and the stress of the contact pressure of the acting point on the spherical surface was found to be:(4) ${\displaystyle p_R=\frac{3R}{2\left[\frac{\pi }{2}+{\theta }_{\text{C}\text{E}}-\alpha +\text{sin}({\theta }_{\text{CE}}-\alpha )\text{cos}({\theta }_{\text{CE}}-\alpha )\right]r^2}}$(4)

This analysis showed that when the line of direction through the coordinate system origin, and parallel to the load direction, was intersected with the acetabular hemispherical surface, this point of intersection coincided with the acting point of the resultant force. In other words, the point of maximum stress of contact pressure on the articular surface acted on the point of intersection. At this moment, ${\displaystyle p_{\text{max}}}$ = ${\displaystyle p_R}$. When the same line of direction was intersected with the articular spherical surface and the intersection point is out of contact surface, the resultant force R did not cross the spherical center. The peak value ${\displaystyle p_{\text{max}}}$ of contact pressure in the articular spherical surface was the stress value of the nearest point to the point of intersection on the edge of the bearing area. At this moment, the maximum contact stress satisfied Eq. (2), i.e., ${\displaystyle p_{\text{max}}=p_R\text{cos}({\theta }_{\text{CE}}-\alpha )}$.

Finite element model and contact stress in the articular surface of the hip joint

Figure 3 shows finite element modeling process. With the total hip being the object of study, the spiral CT was carried out in normal adult males. The scan length was 346 mm and the space between scan layers was 1.5mm. CT images were output in DICOM format and stored in the computer. CT data in DICOM format was imported into the medical image processing software Mimics, and thus the tomographic images could be obtained.

Figure 3. Finite element modeling process.

After CT images were imported into Mimics, the range of image editing was determined by setting the threshold. The data of bone tissues were then extracted based on the difference in gray value between bone tissues and cartilage. The most basic segmentation could be achieved. The segmented 2 parts underwent 3-dimensional reconstruction, and formed the 3-dimensional model of the femur and acetabulum. The cartilage could not be clearly displayed in the CT image for its smaller gray value, so that MRI images were applied for 3-dimensional reconstruction to the articular cartilage.

Firstly, the generated 3-dimensional model was imported into the reverse engineering software Geomagic Studio 2012. The model was modified and optimized. Secondly, the model was matched using a number of curve patches with different sizes and curvatures. The final 3-dimensional geometric model was stored in IGS format. Finally, the femur, articular and acetabular IGS model was imported into Unigraphics NX 8.0 for assembly. The assembly model of the hip joint was then exported in IGS format.

The geometrical model that had been assembled was imported into ANSYS 14.0 s to generate the solid model automatically. Before meshing, appropriate element attributes should be set according to the material properties of bone. Herein, both the cortical bone and articular cartilage were defined as continuous and uniform isotropic elastic materials. Surface-surface contact was assigned between cartilage and cortical bone. The elastic modulus of the cortical bone was defined as 17300MPa, with a Poisson's ratio of 0.29; the elastic modulus of the cartilage was 10MPa, with a Poisson's ratio of 0.4.¹² The element type was defined as a 10-node tetrahedral solid element (10solid187). After the material properties and element type were set, the model underwent automatic meshing. Then the entire sphere of the femoral head underwent mesh refinement by the Remesh command to improve the accuracy of analysis. The final 3-dimensional finite element model was generated, with a total of 6397 elements and 12378 nodes, as shown in Fig. 4.

Figure 4. Finite element model of the hip joint.

The three-dimensional finite element model was used to simulate the state of standing motionless on a single foot. The distal femur was completely fixed, while the acetabulum and the articular cartilage were partially fixed, with movement limited to vertical movement in the coronal plane. The contact force of the hip joint under these conditions was approximately 3 times that of the mass of the body i.e., for a body mass of 70 kg, the size of load forced on the acetabulum was 2100 N.¹³ The load direction was located in the OXZ plane through the center of the femoral head, and the CE angle was set according to the CE angle of the hip joint model.

The finite element model of the total hip was established in this study in order to verify the mathematical model for contact stress in the articular surface of the hip joint. The finite element analysis method was used for stimulation of contact stress in the articular surface of the hip joint. This model was used to obtain the surface stress distribution and peak. This result was compared to the surface stress distribution and peak calculated by the mathematical model in order to evaluate the accuracy of the mathematical model.

Results and discussion

Analysis results of the peak value of contact stress in the articular surface of the hip joint

Table 1 shows the peak values of contact stress in the articular surface of the hip joint with different load directions and CE angles, by the finite element model simulation analysis.

Table 1. Peak stress of different CE angle femoral head surface under different load direction.

θCE=50°		θCE=40°		θCE=30°		θCE=20°		θCE=10°
α/°	pmax/MPa	α/°	pmax/MPa	α/°	pmax/MPa	α/°	pmax/MPa	α/°	pmax/MPa
50	6.6267	—	—	—	—	—	—	—	—
40	5.3350	40	6.2362	—	—	—	—	—	—
30	4.0552	30	5.5851	30	5.4777	—	—	—	—
20	3.7235	20	4.1964	20	4.5185	20	6.2275	—	—
10	3.5930	10	3.5512	10	3.9340	10	5.2382	10	6.8989
0	3.3823	0	3.6331	0	3.6511	0	4.9581	0	4.6097
−10	3.1417	−10	3.7278	−10	3.7802	−10	4.1990	−10	4.0628
−20	3.0415	−20	3.3861	−20	3.4105	−20	4.2246	−20	4.2321
−30	3.2665	−30	3.2475	−30	3.1417	−30	3.4431	−30	3.4669
−40	2.8177	−40	3.4230	−40	2.9128	−40	2.9789	−40	2.9675

Figure 5 shows the changes in the peak values of contact stress in the articular surface of the hip joint with different θ_CE, with changes load direction. Wherein, the data obtained by the finite element analysis was sequentially connected by a straight line. The smooth curve represented the theoretical curve of the mathematical model.

Figure 5. Contrast of simulation results (actual value) and calculated results (theoretical value) in the mathematical model.

The mathematical model showed that for different θ_CE, when the difference of θ_CE –${\displaystyle \alpha }$ was equal, the results of the peak value of contact stress was similar, which was consistent with the calculated results of the simulation model. Meanwhile, the contrast of the mathematical model and simulation results showed that the theoretical curve of the mathematical model was relatively similar to that of the finite element analysis curve, regardless of the trend or value. This verified the accuracy of the mathematical model established in this study.

Distribution of contact stress in the articular surface of the hip joint

From the results in Fig. 6 (a), (b), (c) and (e), when loads were applied in different directions to the finite element model with different θ_CE, all of the peak values of contact stress in the articular surface were near to the intersection of the femoral head surface and the load direction. In addition, they are distributed radially to the periphery, decreasing from the center to the periphery. When the θ_CE was 40°, ${\displaystyle \alpha }$ = −30°, the peak value of the stress was close to the peak value of the stress when the θ_CE was 30°, ${\displaystyle \alpha }$ = −40°. This conclusion was consistent with the results of the mathematical model. However, from the comparison of (c) and (d), it was shown that at a θ_CE of 30°, there was almost no difference in the stress distribution between ${\displaystyle \alpha }$ = − 40° and ${\displaystyle \alpha }$ = −50°. When ${\displaystyle \alpha }$ = −50°, there was still no significant change in the position of the peak value of the stress in the articular surface of the femoral head. When θ_CE was 10°, the deflection angle ${\displaystyle \alpha }$ of the load increased to −60°. As shown in (f), the surface stress was concentrated on the contact edge of the hip acetabulum and the femoral head, which was in contrast to the calculated results of the mathematical model established in this study. This reason for this was that the geometrical model of the hip joint established in this paper was thick at the top and thin in the surrounding area, making the top of the femoral head and femoral articular cartilage in full contact with each other. Simultaneously, the gaps were formed in the peripheral parts. When the deflection angle of the load was gradually increased to a certain value, interspaces began to form in the joint capsule in the acting position of the load. Thus, there was no contact pressure stress generated at the intersection location of the femoral head surface and the load direction. The model was invalid for the contact of the joint capsule, thus the stress in the articular surface of the femoral head was concentrated at the edge of the acetabulum, in contact with the femoral head.

Figure 6. von Mises stress diagram in the articular surface of the hip joint. (A) von Mises stress diagram in the articular surface of the hip joint with the θ_CE = 50° and ${\displaystyle \alpha }$ = 0°; (B) von Mises stress diagram in the articular surface of the hip joint with the θ_CE = 40° and ${\displaystyle \alpha }$ = −30°; (C) von Mises stress diagram in the articular surface of the hip joint with the θ_CE = 30° and ${\displaystyle \alpha }$ = −40°; (D) von Mises stress diagram in the articular surface of the hip joint with the θ_CE =30° and ${\displaystyle \alpha }$ = −50°; (E) von Mises stress diagram in the articular surface of the hip joint with the θ_CE = 20° and ${\displaystyle \alpha }$ =10°; (F) von Mises stress diagram in the articular surface of the hip joint with the θ_CE = 10° and ${\displaystyle \alpha }$ = −60°.

Conclusions

The present study aims to establish a simple model for contact stress on the articular surface of the hip joint, which meets the accuracy requirements. We have mainly solved 2 problems to achieve the objective. Firstly, we determined the expression of simple contact stress in the articular surface of the hip joint. The mathematical model for contact stress in the articular surface of the hip joint was established on the basis of the relationship between contact stress and the geometrical parameters. Secondly, this model was accurately verified by finite element analysis. The following conclusions were achieved through analysis and verification:

The mathematical model could correctly reflect the distribution of stress in the articular surface of the hip joint. When the load deflection angle ${\displaystyle \alpha }$ was greater than −40°, the same characteristics of stress distribution were shown for the mathematical and finite element model. When the load deflection angle ${\displaystyle \alpha }$ was less than −40°, there was a larger deviation between the 2 models due to the contact failure in the finite element model. Typically, the body's normal activities are conducted within in a load deflection range of −20°to 10°, thus the mathematical model could correct reflect the distribution of the stress in the articular surface of the hip joint.

The mathematical model could correctly calculate the peak values of contact stress in the articular surface of the hip joint at different force directions and CE angles. The simulation and mathematical model results of the peak values of contact stress in the articular surface of the hip joint were compared at different directions, including θ_CE = 10°, θ_CE = 20°, θ_CE = 30°, θ_CE = 40°, and θ_CE = 50°. The results were similar in terms of the trends and values.

The mathematical model was simple and satisfied the accuracy requirement. The simulation of the contact stress in the articular surface of the hip joint using the mathematical model could avoid the complex process of the finite element modeling and analysis. In addition, it could also meet the requirement of clinical accuracy in 2 aspects of the values and trends. It could realize the target of using a simple method to achieve the contact stress data without in vitro biomechanical experiments.

Due to the limit of the finite element model, the conclusion “When the direction of the resultant force R is perpendicular to the plane of the outer edge of the bearing area, there was the minimum peak value of contact stress in the articular surface of the femoral head” derived by the mathematical model had not been validated.

Disclosure of potential conflicts of interest

No potential conflicts of interest were disclosed.

Funding

The authors would like to give thanks for the financial support of NSFC (61272387), NCET (NCET-13-0756), State Key Laboratory of Robotics and System (HIT) and Chang Jiang Scholar Candidates Program for Provincial Universities in Heilongjiang (2013CJHB007).