Artículo Científico / Scientific Paper https://doi.org/10.17163/ings.n20.2018.07 pISSN: 1390-650X / eISSN: 1390-860X EQUIVALENT ELASTIC MODULUS FOR PREDICTION OF DEFORMATIONS IN JOINTS MÓDULOS ELÁSTICOS EQUIVALENTES PARA PREDICCIÓN DE DEFORMACIONES EN ARTICULACIONES Franco Marinelli1, Brenda A. Weiss1,*, Marcelo E. Berli1, José Di Paolo1

Table 1. Parameters used for the simulation of knee joints

The hertzian load acting on dry contact in Figure 1b [16] follows a semi-ellipse distribution for the specific load per unit area, as shown in Equation (1):

 (1)

Where q0 is the maximum specific load per unit area, acting in the center of the contact, that is, in the center of the segment of length 2b, and is given by:

 (2)

P is the load per unit length in the direction of the cylinder axis. The value of P is given by the relation between the total load F on the contact and the length of the cylinder H:

 (3)

The term  it is taken as  for tending R1 to infinity. The constants ki of Equation 2 contain the properties of the linear elastic material: Young’s modulus and Poisson’s ratio. Given the parameters in Figure 1b, k2 is null, therefore:

 (4)

Finally, the load q0 is defined by:

 (5)

The half-width b of the contact is given by:

 (6)

and taking into account the previous considerations

 (7)

where the radius R’ [16] is given by:

 (7)

2.2. SoP model

The hip or coxofemoral joint can be physically abstracted as a full spherical body in point contact with a hollow spherical body as shown in Figure 2a. This geometric model is usually called a ball-in-socket type and is characterized by the radius of the femoral head (convex spherical surface), the radius of the acetabulum (concave spherical surface), and the space between both surfaces (clearence c*). The equations that govern this model are expressed in spherical coordinates, therefore, they require a 3-dimensional analysis [11].

An equivalent model, which allows for the expression of the governing equations in Cartesian coordinates— and for the simplification of its resolution by converting it into a problem in 2 dimensions—, is the so-called “ball-on-plane” as shown in Figure 2b [19].

Both a healthy hip joint and a prosthetic joint can be modeled geometrically as a sphere on a deformable plane. This approach is valid when it comes to small contact regions compared to the radius of the spherical surfaces. The equivalent sphere has a R' radius which is obtained from the radius of the femoral head (R2) and the inner radius of the acetabulum (R1), so that the sphere in contact with the plane possesses the combined curvature of the spheres of Figure 2a.

(a)

(b)

Figure 2. The geometry has polar symmetry. (a) Total hip prosthesis [18]. The femoral head and acetabulum are highlighted in pink. (b) Equivalent representation in a SoP contact.

Table 2. Parameters used for the simulation of hip joints

For the determination of the charge on the dry contact surface, the hertzian model was used for a sphere on a plane.

 (9)

Where qo is the maximum load per unit area, acting in the center of the circular contact surface of radius b.

 (10)

P, in this case, is the total charge exerted on the sphere. The radius b, is given by:

 (11)

K being

 (12)

Since E2 is 2 orders of magnitude greater than E1, it could be assumed that E2 tends to infinity, so K becomes:

 (13)

2.3. Equivalent elastic modulus

2.3.1. Column model for layer material (PE)

This model is based on the assumption that the deformable material is composed of columns of cross section to the load (see Figure 3) with differential size. Each column is in lateral contact with columns of the same characteristics, the tangential forces between them being negligible [6].

Figure 3. Layer solid subject to a compression in the direction of the layer.

Assuming the solid is a linear elastic, the relationship between tensions (T) and deformations (E) is given by Hooke’s law which, in terms of the constants of Lamé (μ and ) and in indexical notation, is expressed:

 (14)

where i,k is the Kronecker delta and E and T are the tensors of deformations and tensions, respectively. If it is assumed that the material considered is linear elastic, of low stiffness and layer, loaded to compression in the direction of its thickness (direction y or 2 in indexes), then, the normal tension in the direction y is much greater than the others, which can therefore be neglected. Likewise, deformations other than the relative shortening in the y direction are negligible. Hooke’s law is:

 (15)

which is equivalent to Hooke’s law for a single column of elastic modulus (2μ + ), then (see Figure 1b):

 (16)

 (17)

Where p is the normal tension of the surfaces in contact, d is the surface displacement and Alt is the thickness of the low rigidity material.

 (18)

Where   is the equivalent elastic modulus, which in terms of Young’s modulus (E1) and Poisson’s ratio (v1) of the tibial component is expressed:

 (19)

2.3.2. Winkler equivalent elastic modulus (W) and corrected semi-infinite solid (SIC)

The Winkler equivalent module arises from the case of a solid

limited by a plane and with infinite depth in the direction normal to that plane. In this normal direction, a distributed contact load acts on the plane, which produces displacements at every point of the solid. The elastic problem is described by a potential function from which the tensions and, by virtue of Hooke’s law, the displacements at each point are obtained. The potential function or function of Airy, has one expression for linear contact and another for point contact.

For a load distribution, both for the case of linear contact [16] and point contact [17], the combined displacements of the surfaces that delimit two different solids in contact, are equivalent to twice the displacement of a solid with a «reduced» or «equivalent» elastic modulus, which links the elastic properties of both solids:

 (20)

In the particular case in which both solids in contact are equal, the equivalent modulus E' takes the expression (21) and it is called Winkler:

 (21)

The Winkler modulus is therefore applicable only in the case of natural joints, where both surfaces correspond to the articular cartilage and has been used in scientific publications [7].

As previously discussed, in the case of prosthetic joints, the modulus of elasticity of the femoral component (E2) is 2 orders of magnitude greater than E1, therefore the modulus of elasticity obtained from Equation (20) is.

 (22)

E'SIC is identified as elastic modulus of corrected semi-infinite solid (SIC) since it is another particular case of Equation (20).

Note that  y  are less than for all values of E and v that are considered. Hence, one of the objectives of this work is to determine the validity of using ,  and  in a specific application.

2.4. Displacement models

To determine the displacement suffered by the surfaces in each case to be analyzed, expression (18) was used, for which the contact pressure is the specific load per unit area q of the hertzian contact. In the case of the knee prosthesis, the load varies longitudinally with the coordinate transverse to the axis of the cylinder

(Equation 1), while for the hip joint, the load varies radially (Equation 9). Thus, the displacements in each case are:

 (23)

 (24)

Where d is the displacement of the surface and  is the equivalent elastic modulus, which is varied for the analysis between those obtained for: material of layer (), semi-infinite or Winkler solid () and corrected semi-infinite solid ().

2.5. Parameters for the simulation of knee joints

In the case of knee joints, the parameters shown in Table 1 were used.

It is usually considered that the average maximum strength that a lower limb joint supports is approximately three times the person’s body weight in most daily activities such as walking, climbing stairs and standing on one leg [20, 21]. Considering a person with a body mass (MC) of 80 kg, the total load F on the contact area will be 2352 [N] (F = MC[kg] × 3 × 9, 8[m/s2]). The load per unit length P is calculated according to Equation (3), assuming that the cylinder length H is 0.03 [m].

For the healthy joint, the modulus of elasticity of the articular cartilage is 10 MPa [22], and for the case of the prosthetic joint, it is 1 GPa for the UHMWPE [19]. Both cartilage and UHMWPE have a Poisson’s ratio of 0.4.

For this simulation, three thicknesses are used for the deformable material: 1,2, 2 and 2,5 mm for the joint thickness of both layers of cartilage [23] and 8, 10 and 12 mm for the case of the prosthetic joint made of UHMWPE [24].

2.6. Hip joint

In the case of hip joints, the parameters shown in Table 2 were used.

As in the case of the knee joint, a person with a body mass (MC) of 80 kg is considered and the total load borne by the joint is approximately 3 times the body weight of the person [25].

It should be clarified that the constant c* is the difference between the radius of the acetabulum and the radius of the femoral head (R2R1).

For this simulation, three thicknesses are used for the deformable material: 1,2, 2 and 2,5 mm for the joint thickness of both layers of cartilage and 8, 10 and 12 mm for the case of the prosthetic joint made of UHMWPE.

2.7. Obtaining results through FEM

The calculation by finite elements is done using the Structural Mechanics module of the COMSOL Multiphysics 5.3 licensed software. The analyzed geometries were discretized with Lagrangequadratic elements, tetrahedral for the volumes and triangular for the surfaces. In any case and in the corresponding geometries, the linear elasticity differential equation and the compatibility equations, preloaded in the Structural Mechanics module were solved using the MUMPS direct solver, available in COMSOL Multiphysics. The boundary conditions used are described below.

For the knee joint, the simulation is carried out in a 2D domain (indicated in gray in Figure 4) since the load does not vary in the direction of the axis of the cylinder. The boundary conditions are

 established as follows: the edges L1, L3, L4 and L5 are defined as tension-free, L6 is fixed, that is, with zero displacements, and at the edge L2 the acting tension is the specific hertzian charge q(x) defined in Equation (1) (Figure 4). It should be noted that the used domain is larger than the sector of interest (the one limited by L2), because the high Poisson’s ratio of the simulated materials makes them act almost as incompressible, deforming in the opposite direction at some points in sections L1 and L3. Therefore, a larger domain allows the simulation to comply with the imposed boundary conditions. Figure 4. Scheme of the domain for knee joint where the elastic problem is solved through the FEM. The edges are presented where different edge conditions are applied. In the case of the hip joint, the simulation is performed in a 3D domain (indicated in gray in Figure 5) because the load varies radially. For this case the base of the prism representing the deformable material was fixed and the rest of the faces are free of tensions, except the central circle, of diameter 2b, where the load is exerted. There, the acting tension is the specific hertzian charge q (r) defined in Equation (9). Figure 5. Scheme of the domain for hip joint where the elastic problem is solved through the FEM. 3. Results and discussion 3.1. Knee joint For the case of the healthy joint, where two layers of thin material are assumed to cover the femoral condyles and the tibial plate, the results of displacement of the surfaces in contact are shown in Figure 6. It should be noted that the total displacements of both coatings are equal to one of only two times the thickness (see Equation 23). Therefore, the results shown in Figure 6 correspond to a solid with the thickness of two cartilages. It is worth highlighting the positive displacements predicted by the FEM solution taken as accurate shown in Figure 6, due to the quasi-incompressibility of the assumed material due to the high value of Poisson’s ratio, as discussed in item 2.6. While the simplified model with the Winkler equivalent elastic modulus estimates displacements excessively greater than that those given by the FEM solution,the calculations using the corrected elastic modulus of semi-infinite solid underestimates them for the analyzed cases. The results for the simplified model with the equivalent elastic modulus corresponding to the layer solid are the ones that have better approximation with respect to the values thrown by the FEM solution, in Figure 6 the curves for PE and FEM are superimposed. In Table 3, the maximum displacements obtained for each equivalent elastic modulus and the value obtained by the FEM solution are shown. Table 4 shows the percentage differences corresponding to Table 3, always with respect to the FEM solution. Figure 7 shows the results obtained for the case of the knee prosthesis. Table 5 shows the maximum displacements for each equivalent elastic modulus and those obtained through FEM. Table 6 shows the corresponding percentage difference. The overestimation of the displacements when the Winkler equivalent Young’s modulus is used in the case of the prosthesis is notewhorthy. This is because this modulus of elasticity corresponds to a particular case of the equivalent modulus of elasticity for a semiinfinite solid in which both surfaces in contact have the same mechanical properties, which does not happen in prosthetic joints. However, the validity of the Winkler modulus applied to healthy joints (Figure 6) can also be questioned since it assumes a semi-infinite solid and not a thin layer as it actually the case with natural cartilages. Table 3. Maximum displacement obtained by the FEM and for the simplified model with each equivalent elastic modulus considered for the natural knee joint Table 4. Percentage difference of PE, W and SIC results compared to that obtained by FEM for natural knee joints

 Figure 6. Comparison of the displacement curves of the surfaces in contact obtained for each equivalent elastic module, in contrast to the solution obtained by the FEM in the case of a healthy joint. Table 5. Maximum displacement obtained by FEM and for the simplified model with each equivalent elastic modulus considered in knee prosthesis Figure 7. Comparison of the displacement curves of the surfaces in contact obtained for each equivalent elastic module, in contrast to the solution obtained by the FEM in the case of a prosthetic joint. The comparison of the results of Figures 6 and 7 also reveal the impoverishment of the approximation of layer solid when the elastic modulus of the material and the thickness of the material increase. It is worth noting that UHMWPE’s Young’s modulus is 100 times greater than that corresponding to natural cartilages.

 Table 6. Percentage difference of PE, W and SIC results compared to those obtained by FEM in knee prosthesis 3.2. Hip joint For a healthy hip joint, where two coatings of equal thickness are assumed for which the displacements are added, the surfaces of Figures 8 to 10 were obtained. Table 7 shows the maximum displacements obtained by the simplified model and the different equivalent Young’s modules considered, together with the results obtained by the FEM. Table 8 shows the percentage difference of the results of the simplified model with respect to those obtained by the FEM. Figures 8, 9 and 10 show that for the simplified model, for the cases of layer material and corrected semi-infinite solid, the results have a good approximation to those obtained by the FEM. In the case of the equivalent elastic Winkler modulus, it is shown that the displacement is highly overestimated. This is seen more clearly in Table 8 where, for the layer material, the approximate error is 5% by default with respect to that estimated by the FEM, for the corrected semifinite material the error is close to 15%, also by default, while for the Winkler elastic modulus the displacement values are overestimated with a very high error, which exceeds 70% (Table 9). Table 7. Maximum displacement obtained by FEM and for the simplified natural hip joint model with each equivalent elastic modulus considered Figure 8. Displacement surfaces obtained for each equivalent Young’s modulus and that correspond to the FEM, for healthy natural joint. Alt = 1.5 mm. E = 10 MPa.

 Figure 9. Displacement surfaces obtained for each equivalent Young’s modulus and that correspond to the FEM, for healthy natural joint. Alt = 2 mm. E = 10 MPa. Figure 10. Displacement surfaces obtained for each equivalent Young’s modulus and that correspond to the FEM, for healthy natural joint. Alt = 2.5 mm. E = 10 MPa.

 Figure 11. Displacement surfaces obtained for each equivalent Young’s modulus and that correspond to the FEM, for prosthetic articulation. Alt = 8 mm. E = 1 GPa. Figure 12. Displacement surfaces obtained for each equivalent Young’s modulus and that correspond to the FEM, for prosthetic articulation. Alt = 10 mm. E = 1 GPa.

 Figure 13. Displacement surfaces obtained for each equivalent Young’s modulus and that correspond to the FEM, for prosthetic articulation. Alt = 12 mm. E = 1 GPa. Table 8. Percentage difference of PE, W and SIC results compared to those obtained by FEM for natural hip joint Table 9. Maximum displacement obtained by FEM and for each equivalent elastic modulus considered for the case of prosthetic joint Figures 11 to 13 show the results corresponding to hip prostheses. In this case, discussions for healthy hip articulation are also valid, except that the difference or percentage error relationships are different. Comparing the maximum displacements that occur in the center of the contact, it is observed that for the layer and corrected semi-infinite materials the results of the simplified model differ around 20% and 30%, respectively, by default with respect to that obtained through the FEM. For the Winkler’s equivalent elastic modulus, on the other hand, an error of 38% and 48% is presented (Table 10). This is due to the fact that with the increase of the thickness of the solid, the approximation of semi-infinite solid is improved, while the SIC and the layer approach worsen. In addition, the PE approach, as in the case of the knee, is worsened by the increase in the elastic modulus between the natural cartilage and the UHMWPE of the prosthetic joint. Table 10. Percentage difference of the PE, W and SIC results with respect to that obtained by the FEM for the case of prosthetic joint 4. Conclusions An analysis of deformations has been presented from the calculation of displacements of the surfaces of two solids in charged contact, to evaluate the relevance of using simplified models when constructing contact models. The study was motivated by human joints, where the displacement of surfaces in contact is fundamental to enable a lubrication flow that acts to minimize wear. For the prediction of the displacements, three equivalent Young’s modulus were evaluated in a simplified column model: for a layer solid, for a semi-infinite or Winkler solid, and for a corrected semi-infinite solid. The results obtained for the Young equivalent module corresponding to a material of layer, are those that better approximate those obtained by the FEM in both geometries considered: linear contact (equivalent to the knee) and point contact (equivalent to the hip). However, the percentage differences with respect to the FEM solution vary in each case and for each material analyzed: natural cartilage (less than 5% in the hip and less than 1% in the knee) or UHMWPE in the prosthetic joint (less than 23, 3% in the

 is worsened by the increase in the elastic modulus between the natural cartilage and the UHMWPE of the prosthetic joint. It was demonstrated that the equivalent Young modules derived from the semi-infinite solid approximation, such as the Winkler modulus and the SIC modulus, are inappropriate due to the excessive overestimation of the calculated displacements obtained when considering the first (greater than 38% and reaching the 70% in several cases) and the underestimation of the displacements when considering the second (higher than 10% and reaching 30% in some cases). This, by virtue of what was said in the previous paragraph, puts the analysis at risk for the Winkler’s case since, in complete lubrication models, there is a risk of predicting separations between the surfaces in contact that are greater than what would actually occur. Finally, the work was very useful for the training in the arts of Computational Mechanics of an advanced student of the Bioengineering major, in their role as a research fellow at Universidad Nacional de Entre Ríos, Argentina. Acknowledgements To Universidad Nacional de Entre Ríos for financing through PID 6162 and to CONICET for funding through a doctoral grant. References [1] K. J. Bozic and M. D. Ries, “Wear and osteolysis in total hip arthroplasty,” Seminars in Arthroplasty, vol. 16, no. 2, pp. 142–152, 2005, cOMPLICATIONS AFTER TOTAL HIP ARTHROPLASTY. [2] R. Pivec, A. J. Johnson, S. C. Mears, and M. A. Mont, “Hip arthroplasty,” The Lancet, vol. 380, no. 9855, pp. 1768–1777, 2012. [3] S. Kurtz, K. Ong, E. Lau, F. Mowat, and M. Halpern, “Projections of primary and revision hip and knee arthroplasty in the united states from 2005 to 2030,” The Journal of bone and joint surgery, vol. 89, pp. 780–785, 04 2007. [Online]. Available: https://goo.gl/1qXADn [4] Y. Su, P. Yang, Z. Fu, Z. Jin, and C. Wang, “Time-dependent elastohydrodynamic lubrication analysis of total knee replacement under walking conditions,” Computer Methods in Biomechanics and Biomedical Engineering, vol. 14, no. 6, pp. 539–548, 2011. [5] M. Mongkolwongrojn, K. Wongseedakaew, and F. E. Kennedy, “Transient elastohydrodynamic lubrication in artificial knee joint with nonnewtonian fluids,” Tribology International, vol. 43, no. 5, pp. 1017–1026, 2010, special Issue on Second International Conference on Advanced Tribology (iCAT2008). [6] J. Di Paolo, D. M. Campana, S. Ubal, and M. E. Berli, “Flujos de lubricación en canales elásticos. una experiencia didáctica en clases de mecánica del continuo para bioingeniería,” in VI Congreso Argentino de Enseñanza en Ingeniería, 2008. [Online]. Available: https://goo.gl/LLLpNC [7] M. Yousfi, “Etude biomécanique de l’articulation du genou humain. caractérisation mécanique et modélisation de l’ecoulement du fluide synovial en ecrasement lors d’un cycle de marche,” Master’s thesis, Université du 20 Août 1955 Skikda, 2014. [Online]. Available: https://goo.gl/GQJmsr [8] T. J. R. Hugues, The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, C. Corporation, Ed., 2012. [Online]. Available: https://goo.gl/ktCM6S [9] K. V. Nemade and V. K. Tripathi, “A mathematical model to calculate contact stresses in artificial human hip joint,” International Journal of Engineering Research and Development, vol. 6, no. 12, pp. 119–123, 2013. [Online]. Available: https://goo.gl/FyQQpX [10] F. Li, H. Chen, and K. Mao, “Computational simulation analysis for torus radius of edge contact in hip prostheses,” Acta of bioengineering and biomechanics, vol. 17, no. 3, pp. 67–73, 2015. [Online]. Available: https://goo.gl/8AJ7nz [11] Q. Meng, F. Liu, J. Fisher, and Z. Jin, “Contact mechanics and lubrication analyses of ceramic-onmetal total hip replacements,” Tribology International, vol. 63, pp. 51–60, 2013, the International Conference on BioTribology 2011. [12] J. Armas Sánchez, “Análisis de la fractura por contacto mecánico y su solución por análisis numérico,” Master’s thesis, Instituto Politécnico Nacional. México, 2015. [Online]. Available: https://goo.gl/LnEdPf [13] J. Di Paolo and M. E. Berli, “Numerical analysis of the effects of material parameters on the lubrication mechanism for knee prosthesis,” Computer Methods in Biomechanics and Biomedical Engineering, vol. 9, no. 2, pp. 79–89, 2006, pMID: 16880159. [14] M. E. Berli, D. M. Campana, S. Ubal, and J. Di Paolo, “Lubrication model of a knee prosthesis, with non newtonian fluid and porous rough material,” Latin American applied research, vol. 39, no. 2, pp. 105–111, 2009. [Online]. Available: https://goo.gl/aKL9Ky [15] B. A. Weiss, M. E. Berli, S. Ubal, and J. Di Paolo, “umerical solution of a 2d lubrication model with sommerfeld boundary conditions for hip prostheses,” in VI Latin American Congress on Biomedical Engineering CLAIB, Paraná, Argentina 29, 30 & 31 October, vol. 49, 2014. [16] D. Dowson and G. R. Higginson, Elasto-Hydrodinamic Lubrication. Pergamon Press, 1977. [Online]. Available: https://goo.gl/NV522s [17] H. P. Evans and W. Snidle, “The elastohydrodynamic lubrication of point contacts at heavy loads,” Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, vol. 382, no. 1782, pp. 183–199, 1982.

 [18] MBA. (2017) More power in your hands. MBA Surgical Empowerment. [Online]. Available: https://goo.gl/aPngur [19] F. Mattei, Land Di Puccio, B. Piccigallo, and E. Ciulli, “Lubrication and wear modelling of artificial hip joints: A review,” Tribology International, vol. 44, no. 5, pp. 532–549, 2011, special Issue: ECOTRIB 2009. [20] I. Kutzner, B. Heinlein, F. Graichen, A. Bender, A. Rohlmann, A. Halder, A. Beier, and G. Bergmann, “Loading of the knee joint during activities of daily living measured in vivo in five subjects,” Journal of Biomechanics, vol. 43, no. 11, pp. 2164–2173, 2010. [21] A. Sharma, “A method to calculate the femoropolyethylene contact pressures in total knee arthroplasty in vivo,” Master’s thesis, University of Tennessee, United States, 2005. [Online]. Available: https://goo.gl/QucK2c [22] J. S. Jurvelin, M. D. Buschmann, and E. B. Hunziker, “Mechanical anisotropy of the human knee articular cartilage in compression,” Proceedings of the Institution of Mechanical Engineers, Part H: Journal of Engineering in Medicine, vol. 217, no. 3, pp. 215–219, 2003, pMID: 12807162. [23] D. E. T. Shepherd and B. B. Seedhom, “Thickness of human articular cartilage in joints of the lower limb,” Annals of the Rheumatic Diseases, vol. 58, no. 1, pp. 27–34, 1999. [24] H. Oonishi, H. Iwaki, N. Kin, S. Kushitani, N. Murata, S. Wakitani, and K. Imoto, “The effects of polyethylene cup thickness on wear of total hip prostheses,” Journal of Materials Science: Materials in Medicine, vol. 9, no. 8, pp. 475–478, 1998. [25] G. Bergmann, G. Deuretzbacher, M. Heller, F. Graichen, A. Rohlmann, J. Strauss, and G. N. Duda, “Hip contact forces and gait patterns from routine activities,” Journal of Biomechanics, vol. 34, no. 7, pp. 859–871, 2001.