FACTA UNIVERSITATIS Series: Mechanical Engineering Vol. 17, N o 1, 2019, pp. 29 - 38 https://doi.org/10.22190/FUME190122014E © 2019 by University of Niš, Serbia | Creative Commons License: CC BY-NC-ND Original scientific paper  MULTILEVEL NUMERICAL MODEL OF HIP JOINT ACCOUNTING FOR FRICTION IN THE HIP RESURFACING ENDOPROSTHESIS Galina M. Eremina, Alexey Yu. Smolin Institute of Strength Physics and Materials Science SB RAS, Tomsk, Russia, Tomsk State University, Tomsk, Russia Abstract. Friction between the moving parts of the endoprosthesis has a significant impact on the endoprosthesis operation time. Primarily, it concerns the endoprosthesis of hip and knee joints. To improve the tribological characteristics of the metal endoprosthesis, hardening nanostructured coatings are used. Usually, titanium and titanium alloys are used as metal, and titanium nitride is used as a coating. Herein, we propose an approach to multilevel modeling of the system “bone-endoprosthesis” which is based on the movable cellular automaton method and accounts for friction between the moving parts of the hip resurfacing endoprosthesis. We validated the models of the friction system materials using the instrumented scratch test simulation. Then, we simulated friction at the mesolevel, explicitly considering roughness of the coating. The results obtained at the mesolevel were used as tribological characteristics of the coating in the macroscopic model of the hip resurfacing endoprosthesis. Key Words: Hip Resurfacing Endoprosthesis, Coatings, Scratch Test, Friction, Simulation, Movable Cellular Automaton Method 1. INTRODUCTION Prosthetics is an important part of the modern traumatology and orthopedics, especially for large human joints such as hips and knees. Nowadays there are two options for hip joint endoprostheses: total hip replacement and hip resurfacing arthroplasty. In the first option the joint is replaced by endoprosthesis, in the second, just the surface of the femur head is replaced by the cap, which is specially shaped like a mushroom. In spite of the surgical complexity, the second approach is much better suited for young patients due Received January 22, 2019 / Accepted March 24, 2019 Corresponding author: Galina M. Eremina Institute of Strength Physics and Materials Science SB RAS, pr. Akademicheskii 2/4, Tomsk, 634055, Russia E-mail: anikeeva@ispms.ru 30 G.M. EREMINA, A.YU. SMOLIN to their activity. Mainly, both the resurfacing cup and the matching cup placed in the acetabulum are made of cobalt-chrome metal alloy. Sliding friction between the cups leads to wear of the metal. The wear debris generated and released from the metal parts of the endoprosthesis accumulates in the adjacent tissue. This may lead to necrosis of the tissue cells and small phagocytes; an expressed allergic vasculitis is also observed in the literature [1]. The surface layer structure of the contacting materials has a major influence on the wear process. Therefore, in the last decade, one can observe extensive research in the field of applying ultrathin hardening coatings and layers on titanium alloys used for resurfacing hip endoprosthesis. In this case, titanium alloys are used as metal and titanium nitride (TiN) is used as the coating. The structure of the coating is determined by the regimes of its deposition. Usually, the following technologies are used for this purpose: physical vapor deposition (PVD), chemical vapor deposition (CVD) and powder immersion reaction assisted coating (PIRAC) [2, 3]. The coating obtained by PVD and CVD techniques have low adhesion and may delaminate under dynamic loading, while adhesion of the coating obtained by PIRAC is much higher. Furthermore, PIRAC allows depositing a wide range of compositions such as carbides, borides, and nitrides [4, 5]. That is why PIRAC appears to be very promising for producing hardening coating of metal implants. Testing of prostheses has several stages including preclinical and clinical trials. Clinical trials are carried out through the prosthesis installation in a living human body and are performed at the final stage of the testing. It is important to understand that the installation of low quality or poorly designed prosthesis may lead to problems of the patient health. That is why the preclinical study is particularly important. Preclinical studies of the mechanical behavior of the prosthesis can be divided into experimental and theoretical ones. Experimental studies are tests using a special technological unit that simulates the dynamic loading experienced by prosthesis in a living body. To determine mechanical properties of the material surface such standard methods as depth-sensing indentation and scratch testing as well as three-point bending are used. It is worthy to note that usually experimental studies are performed in standard atmospheric conditions and their results may significantly differ from the material behavior in the conditions of the human body. That is why part of the experiments is conducted in the condition mimicking the living body (in vitro) or using special samples installed into an animal body (in vivo). However, it is obvious that the number of in vivo tests and their capabilities is considerably restricted. Theoretical studies of the mechanical behavior of the prosthesis using computer simulation make it possible to study the mechanical behavior of the prosthesis taking into account the influence of various factors. Therefore, in preclinical studies, computer simulation is used to predict the mechanical behavior of prostheses. For example, computer simulation based on numerical methods of continuum mechanics allows detailed investigations of the behavior of coating and surface layer under contact loading. The influence of the coating thickness on the load-displacement curve for the instrumented indentation of ceramic coatings on the metal substrate is studied in [6-8]. The role of the surface roughness of ceramic coating in the mechanical behavior of the system “coating- substrate” in instrumented indentation is numerically studied in [9-12]. Both instrumented indentation and scratch tests of the hardening coating are simulated in [13-17], where the authors analyze stress fields in the contact region as well as the peculiarities of the loading- displacement curves. Multilevel Numerical Model of Hip Joint Accounting for Friction in Hip Resurfacing Endoprosthesis 31 The majority of papers on simulation of the mechanical behavior of the system “bone- endoprosthesis” also use numerical methods of continuum mechanics. In [18], the authors study the stress fields in the hip bone for two types of the endoprostheses (total replacement and resurfacing) and show that resurfacing leads to minimal changes of the stress distribution in the bone compared with the original “healthy” bone. Dickinson et al. studied the mechanical behavior of the system “bone-endoprosthesis” for different structure and geometry of the implant [19] as well as its material [20]. In spite of a large number of the papers devoted to the modeling of the mechanical behavior of the hip endoprosthesis there is still no unified multiscale theoretical study of the system with explicit consideration of structural peculiarities of the materials. In such a model, the studies at each scale should solve their specific problems and their results should be used for modeling at the upper scale finally providing the capability to design personalized endoprosthesis. The purpose of this paper is the development of a multiscale numerical model of the mechanical behavior of the system “bone-endoprosthesis” for hip resurfacing arthroplasty that considers explicitly such geometrical parameters of the TiN coating as its thickness, roughness, and mechanical properties. To reach this purpose, we solved three tasks. Firstly, we validated the models of the materials used in the friction pair of the endoprosthesis. For this, we developed a numerical model for simulating instrumented indentation and scratch testing of the hardening coating on the titanium substrate. The simulation results were compared against available experimental data. Secondly, we developed a model for sliding friction at the mesoscale, which considers explicitly roughness and mechanical properties of the contacting surfaces. The simulation at this scale then provided an effective coefficient of friction that can be used at macroscale modeling. The final third task was the development of a numerical macromodel of the system “bone-endoprosthesis” and studying this model. 2. DESCRIPTION OF THE SIMULATION METHOD Our numerical models are based on the method of movable cellular automata (MCA), which is a representative of so-called discrete element methods (DEM) and differs substantially from numerical methods used in continuum mechanics [21]. The MCA considers a material as a set of particles of finite size (automata) interacting among each other and thereby simulating deformation and fracture of the real material. Translational motion and rotation of the movable automata are governed by the Newton-Euler equations. The main advantage of the MCA-method in comparison with the DEM are the generalized many-body formulas for the central interaction forces acting between the particle pair of s similar to the embedded atom force filed used in molecular dynamics. It is based on computing components of the average stress and strain tensors in the bulk of automaton according to the homogenization procedure described in [21]. The use of many-body interaction forces allows, within the discrete element approach, the correct simulation of important features of the mechanical behavior of solids such as Poisson effect and plastic flow. To describe the elastoplastic flow of the material, it was proposed to use the plastic flow theory (namely von Mises model) by adapting the algorithm of Wilkins in the MCA method [21, 22]. The radial return algorithm of Wilkins consists of the solution of the 32 G.M. EREMINA, A.YU. SMOLIN elastic problem in increments and the subsequent “drop” of components of stress deviator tensor to von Mises yield surface in the case that equivalent stress exceeds it. A pair of automata can be in one of two states: bound and unbound. Thus, in the MCA fracture and coupling of fragments (crack healing, microwelding etc) are simulated by the corresponding switching of the pair state. The switching criteria depend on physical mechanisms of material behavior [21, 22]. Note that the knowledge of the stress and strain tensor in the bulk of an automaton allows the direct application of conventional fracture criteria written in the tensor form. Herein, we used the von Mises criterion based on a threshold value of the equivalent stress. 3. VALIDATION OF THE MATERIALS MODEL From the literature [23], we chose the following values for the material properties of the titanium alloy Ti6Al4V: density ρ = 4420 kg/m 3 , shear modulus G = 41 GPa, bulk modulus K = 92 GPa, Young’s modulus E = 110 GPa, yield stress σy = 0.99 GPa, ultimate strength σb = 1.07 GPa and ultimate strain εb = 0.10. Geometric features of the TiN coating and its mechanical properties are determined by the deposition modes at PIRAC [24] forming. Thus, at a deposition temperature of 700º C and a treatment time of 48 hours, the coating on the titanium substrate has a thickness of 1.3 μm with an average roughness height of 0.15 μm (mode 1) and the elastic modulus, depending on the deposition regime, was equal to E1 = 258 GPa; at a deposition temperature of 800º C and a treatment time of 4 hours the coating on the titanium substrate has a thickness of 1.4 μm with an average roughness height of 0.132 μm (mode 2) and an elastic modulus of E2 = 258 GPa; at a deposition temperature of 900° C and a treatment time of 2 hours the coating on the titanium substrate has a thickness of 1.5 μm with an average roughness height of 0.265 μm (mode 3) and an elastic modulus of E3 = 321 GPa. According to this information, we chose the following values for the material properties of the TiN coating: ρ = 52200 kg/m 3 , G1 = 104 GPa, G2 = 104 GPa, G3 =129 GPa, K1 = 173 GPa, K2 = 173 GPa, K3 = 205 GPa. Data for yield stress and ultimate strength and strain were obtained using reverse analysis of the load-displacement curve for sensing indentation: σy = 4.50 GPa, σb = 5.50 GPa and εb = 0.075 [25]. A general view of the model for sensing indentation is shown in Fig. 1a. The model specimen was a parallelepiped consisting of titanium substrate, interface and TiN coating. Loading was simulated by moving all the automata of the Berkovich indenter with constant velocity Vz = −1 m/s for loading until the required penetration depth, and Vz = 1 m/s for unloading (Fig. 1, b). Applying the procedure of Oliver and Pharr [26] to the simulation results provided the dependence of the material hardness on the penetration depth (Fig. 2) for three kinds of coatings obtained in different regimes of deposition. It can be seen from Fig. 2 that the hardest is the coating obtained by mode 3, all the curves correspond to experimental data from [24]. Thus, we may conclude that the models of the materials behavior are validated and can be used in further steps in our work. The model specimen for scratch testing was also a parallelepiped consisting of titanium substrate, interface and TiN coating but elongated along axis Y (Fig. 3). To simulate the force acting on indenter along axis Z in the experiment, we set the velocity of indenter automata to VZ = −0.5 m/s (Fig. 3, b) until the indenter was immersed into a Multilevel Numerical Model of Hip Joint Accounting for Friction in Hip Resurfacing Endoprosthesis 33 predetermined penetration depth; after that the vertical velocity was set to zero. Here, we considered penetration of the indenter only up to the interface layer; the possibility of the coating to detach from the substrate was not allowed. To move the indenter along the sample surface, a constant velocity of the indenter automata along axis Y was set to VY = 1 m/s. a b Fig. 1 The model specimen for indentation (a) and its cross-section with loading parameters (b) represented by automata packing (the numbers indicate the model materials: 1 – titanium, 2 – interface, 3 – coating, 4 – diamond) Fig. 2 Dependence of the hardness on the depth of penetration (the numbers indicate the modes of coating deposition) Based on the results of our simulations, images of the deformed sample were created, and the values of the critical force characteristics at certain stages of fracture were obtained. It was found that the value of the critical force for the onset of fracture and cracking is largest for the coating thickness of 1.4 μm and roughness of 0.132 μm, and smallest for the thickness of 1.5 μm and roughness of 0.265 μm. At delamination of the coating, the maximum of the critical force is typical for the specimen with a coating thickness of 1.5 μm and roughness of 0.265 μm, and the minimum for the specimen with a coating thickness of 1.3 μm and roughness of 0.15 μm (Fig. 4). 34 G.M. EREMINA, A.YU. SMOLIN a b Fig. 3 The model specimen for scratch testing (a) and its cross-section with loading parameters (b) represented by automata packing (the numbers indicate the model materials: 1 – titanium, 2 – interface, 3 – coating, 4 – diamond) Fig. 4 The force acting on the indenter from coating P versus scratching path h of the indenter (the numbers indicate the modes of coating deposition) 4. MESOMODEL FOR FRICTION In the next step, we developed a mesomodel of sliding friction of the contacting surfaces of the endoprosthesis that explicitly accounts for surface roughness, the mean height of which is determined by the regime of the coating deposition (Fig. 5). At the macroscale we merely need the effective value of the friction coefficient. Therefore, the substrate and interface were not considered in this model. The scheme of loading is shown in Fig. 5b and is as follows. The automata of the bottom layer were fixed, the automata of the top layer were moved with constant horizontal velocity VY = 1 m/s and subjected to pressure P = 0.75σy. Multilevel Numerical Model of Hip Joint Accounting for Friction in Hip Resurfacing Endoprosthesis 35 a b Fig. 5 The model specimen for sliding friction at mesoscale (a) and its cross-section with loading parameters (b), the colors indicate different bodies Figure 6 depicts plots of the friction coefficients versus the calculation time obtained from the simulation at mesoscale for three different TiN coatings, whose parameters and roughness correspond to the different deposition regimes. Average values of these friction coefficients at the final steady stage were then used in the model at the macroscale. Fig. 6 Plots of the friction coefficient in the model at mesoscale versus the calculation time (the numbers indicate the modes of coating deposition) 5. MACROMODEL OF HIP RESURFACING ENDOPROSTHESIS Finally, we developed a macromodel of the hip resurfacing endoprosthesis with part of the femur bone shown in Fig. 7. The geometry of the model was based on the so-called 3rd generation composite femur [27], which provides geometries of cortical bone and cancellous bone as different solid bodies. For our purpose, we cut the top part of the bone 36 G.M. EREMINA, A.YU. SMOLIN geometry, added resurfacing endoprosthesis (colored in cyan, and its coating colored in blue) and special loading part (colored in red) (Fig. 7, b,c). a b c Fig. 7 The model of the hip resurfacing endoprosthesis (a) its cross-section (b) and scheme of loading (c) a b b d Fig. 8 Distributions of the mean stress in the system “bone-endoprosthesis” loaded by a force of 3 kN (a,c) and 10 kN (b,c) for titanium endoprosthesis (a,b) and titanium with TiN coating (c,d) Multilevel Numerical Model of Hip Joint Accounting for Friction in Hip Resurfacing Endoprosthesis 37 Based on the developed model, we simulated compression of the part of the femur with resurfacing endoprosthesis with hardening coating and without it. The loading was applied by setting a constant velocity in both horizontal and vertical directions as shown in Fig. 7c up to reaching the critical values of the resisting force of 3 kN and 10 kN [28]. The simulation results are presented in Fig. 8 as distributions of mean stress. It is reported that the strength limit for the cancellous bone of healthy people reaches about 10 MPa in compression and 5 MPa in tension; it may be lower than 5 MPa and 3 MPa for some diseases correspondingly [29]. Analysis of stress distribution in the model system showed that the maximum tensile stress was observed near the endoprosthesis and did not reach critical values. Therefore, we tried to analyze the compression stress in the scale up to 5 MPa. From Fig. 8, one can see that the maximum compression stress that can lead to fracture is observed in the femoral neck. 6. CONCLUSIONS A multiscale numerical model for analyzing mechanical behavior of the hip resurfacing endoprosthesis is proposed. The material models are validated against the available experimental data using the simulation of instrumented indentation and scratching. The macroscopic coefficient of friction in the friction pair of the endoprosthesis is proposed to be found from the model for sliding friction at mesoscale where roughness of the contacting surfaces is considered explicitly. Simulation at the macroscale includes a part of the femur bone with resurfacing endoprosthesis. An analysis of the stress distribution in the macromodel reveals that the critical compression stress is observed in the femoral neck. Further studies could include more realistic loading mimicking not only compression but also rotation of the femoral head in the acetabulum prosthesis cup. Acknowledgements: The reported study was funded by RFBR according to the research project № 18- 38-00323 (the models developed and the simulations) and framework of the Russian Fundamental Research Program of the State Academies of Sciences for 2013–2020 (Priority direction III.23, the modifications of the MCA software for importing CAD models and modeling friction). REFERENCES 1. Bitar, D., Parvizi, J., 2015, Biological response to prosthetic debris, World Journal of Orthopaedics, 6(2), pp. 172-89. 2. Sovak, G., Weiss, A., Gotman, I., 2000, Osseointegration of Ti6A14V alloy implants coated with titanium nitride by a new method, The Journal of Bone & Joint Surgery Series B, 82(2), pp. 290-296. 3. Gotman, I., Gutmanas, E.Y., Hunter, G., 2017, Wear-resistant ceramic films and coatings, Comprehensive Biomaterials II, pp.165-203. 4. Wu, S. J., Li, H., Wu, S. Y., Guo, Q. and Guo, B., 2014, Preparation of titanium carbide–titanium boride coatings on Ti6Al4V by PIRAC, Surface Engineering, 30(9), pp. 693-696. 5. Wu, S, Ma, S, Wu, S, Zhang, G, Dong, N., 2018, Composition, microstructure, and friction behavior of PIRAC chromium carbide coatings prepared on Q235 and T/P 24, International Journal of Applied Ceramic and Technology, 15, pp. 501–507. 6. Kot, M., Rakowski, W., Lackner, J.M., Major, T., 2013, Analysis of spherical indentations of coating- substrate systems: Experiments and finite element modeling, Materials and Design, 43, pp. 99-111. 7. Lofaj, F., Németh, D., 2017, The effects of tip sharpness and coating thickness on nanoindentation measurements in hard coatings on softer substrates by FEM, Thin Solid Films, 644, pp.173-181. 38 G.M. EREMINA, A.YU. SMOLIN 8. Marchiori, G., Lopomo, N., Boi, M., Berni, M., Bianchi, M., Gambardella, A., Visani, A., Russo, A., Marcacci, M., 2016, Optimizing thickness of ceramic coatings on plastic components for orthopedic applications: A finite element analysis, Materials Science and Engineering C, 58, pp. 381-388. 9. Bouzakis, K.-D., Michailidis, N., Hadjiyiannis, S., Skordaris, G., Erkens, G., 2002, The effect of specimen roughness and indenter tip geometry on the determination accuracy of thin hard coatings stress–strain laws by nanoindentation, Materials Characterization, 49(2), pp. 149-156. 10. Jiang, W.-G., Su, J.-J., Feng, X.-Q., 2008, Effect of surface roughness on nanoindentation test of thin films, Engineering Fracture Mechanics, 75(17), pp. 4965-4972. 11. Sliwa, A., Mikuła, J., Gołombek, K., Tanski, T., Kwasny, W., Bonek, M., Brytan, Z., 2016, Prediction of the properties of PVD/CVD coatings with the use of FEM analysis, Applied Surface Science, 388, Part A, pp. 281-287. 12. Skordaris, G., Bouzakis, K.D., Kotsanis, T., Charalampous, P., Bouzakis, E., Breidenstein, B., Bergmann, B., Denkena, B., 2017, Effect of PVD film's residual stresses on their mechanical properties, brittleness, adhesion and cutting performance of coated tools, CIRP Journal of Manufacturing Science and Technology, 18, pp. 145-150. 13. Zlotnikov, Ig., Dorogoy, A, Shilo, D., Gotman, I., Gutmanas, E., 2010, Nanoindentation, modeling, and toughening effects of zirconia/organic nanolaminates, Advanced engineering materials 12(9), pp. 935-941. 14. Li, J, Beres, W., 2006, Three-dimensional finite element modelling of the scratch test for a TiN coated titanium alloy substrate, Wear, 260, pp. 1232-1242. 15. Holmberg, K, Laukkanen, A, Ronkainen, H., Wallin, K., Varjus, S., 2003, A model for stresses, crack generation and fracture toughness calculation in scratched TiN-coated steel surfaces, Wear, 254, pp. 278-291. 16. Pandure, P., Jatti, V., Singh, T., 2014, Three dimensional FE modeling and simulation of nano-indentation and scratch test for TiN coated high speed steel substrate, International Journal of Applied Engineering Research, 9, pp. 2771-2777. 17. Toparlj, M., Sasaki, Sh., 2002, Evaluation of the adhesion of TiN films using nanoindentation and scratch testing, Philosophical Magazine A, 82(10), pp. 2191-2197. 18. Kuhl, E., Balle, F., 2005, Computational modeling of hip replacement surgery: total hip replacement vs. hip resurfacing, Technische Mechanik, 25(2), pp. 107-114. 19. Dickinson, A., Taylor, A., Browne, M., 2012, Implant-bone interface healing and adaptation in resurfacing hip replacement, Computer Methods in Biomechanics and Biomedical Engineering, 15(9), pp. 935-947. 20. Dickinson, A.S., Browne, M., Roques, A.C., Taylor, A.C., 2013, A fatigue assessment technique for modular and pre-stressed orthopaedic implants, Medical Engineering and Physics, 36(1), pp. 72-80. 21. Shilko, E.V., Psakhie, S.G., Schmauder, S., Popov, V.L., Astafurov, S.V., Smolin, A.Yu., 2015, Overcoming the limitations of distinct element method for multiscale modeling of materials with multimodal internal structure, Computational Materials Science, 102, pp. 267-285. 22. Smolin, A.Yu., Shilko, E.V., Astafurov, S.V., Kolubaev, E.A., Eremina, G.M., Psakhie, S.G., 2018, Understanding the mechanisms of friction stir welding based on computer simulation using particles, Defence Technology, 14, pp. 643-656. 23. Material Datasheet: TIMETAL 6-4s, Titanium Metals Corporation, 2000. 24. Bonello, T., Avelar-Batista Wilson, J.C., Housden, J., Gutmanas, E.Y., Gotman, I., Matthews, A., Leyland, A., Cassar, G., 2014, Evaluating the effects of PIRAC nitrogen-diffusion treatments on the mechanical performance of Ti–6Al–4V alloy, Materials Science & Engineering A, 619, pp. 300-311. 25. Giannakpoulos, A.E., Suresh, S., 1999, Determination of elastoplasic properties by instrumented sharp indentation, Scripta Materialia. 40(10), pp. 1191-1198. 26. Oliver, W., Pharr, G. M., 1992, An improved technique for determining hardness and elastic modulus using load and displacement sensing indentation experiments, Journal of Materials Research, 7(6), pp.1564-1583 27. Cheung, G., Zalzal, P., Bhandari, M., Spelt, J.K., Papini, M., 2004, Finite element analysis of a femoral retrograde intramedullary nail subject to gait loading, Medical Engineering & Physics, 26(2), pp. 93-108. 28. Todo, M., 2018, Biomechanical analysis of hip joint arthroplasties using CT-image based finite element method. Journal of Surgery and Research, 1, pp. 34-41. 29. Gerhardt, L.C., Boccaccini, A.R., 2010, Bioactive glass and glass-ceramic scaffolds for bone tissue engineering, Materials, 3(7), pp. 3867-3910.