Next Article in Journal
A Numerical Model for the Analysis of the Bearings of a Diesel Engine Subjected to Conditions of Wear and Misalignment
Previous Article in Journal
Grease Performance Requirements and Future Perspectives for Electric and Hybrid Vehicle Applications
Previous Article in Special Issue
Analysis of Friction in Total Knee Prosthesis during a Standard Gait Cycle
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Mathematical Development of a Novel Discrete Hip Deformation Algorithm for the In Silico Elasto-Hydrodynamic Lubrication Modelling of Total Hip Replacements

by
Alessandro Ruggiero
* and
Alessandro Sicilia
Department of Industrial Engineering, University of Salerno, 84084 Fisciano, Salerno, Italy
*
Author to whom correspondence should be addressed.
Lubricants 2021, 9(4), 41; https://doi.org/10.3390/lubricants9040041
Submission received: 9 February 2021 / Revised: 23 March 2021 / Accepted: 7 April 2021 / Published: 9 April 2021

Abstract

:
In this paper, the procedure to achieve an accurate deformation model of a total hip replacement (THR) was proposed with the aim to obtain a numerical tool to be simply merged into THR elasto-hydrodynamic computational synovial lubrication algorithms. The approach was based on the Finite Element Method (FEM) and was developed in a Matlab code, allowing the definition of the influence matrix and of a boundary conditions vector. It works with linear tetrahedra and performs the displacement calculation for both the acetabular cup and the femoral head, taking into account the anatomical hip relative motion, by coupling them with a cubic interpolation matrix. Two simulations were conducted in order to validate the algorithm and the results were compared with the ones obtained by the commercial software Ansys. The comparison provides a satisfactory agreement in terms of surface deformation, Von Mises stress and strain energy, proving the reliability of the model and the possibility to use the model in the in silico prostheses tribological simulations, avoiding the complexity and the high computational resource requirement coming from the coupling between complex lubrication algorithms and FEM commercial software, and with the possibility to directly act on many key parameter characteristics of the investigated problem.

1. Introduction

In biomechanics, the human joints can be physically modelled as coupled surfaces separated by the lubricant synovial fluid, so the lubrication theory is fundamental to develop numerical algorithms able to predict interesting tribological phenomena, like the joint eccentricity, synovial pressure and thickness, surface wear, etc [1,2,3,4].
The operating conditions in which the lubricated joints operate are represented by the load acting on the joint (dynamics) and the relative motion between the linked bodies (kinematics). Based on the operating conditions, three lubrication modes can be distinguished [5,6,7]:
  • When the joint is subjected to intense relative motion and light loading, the hydrodynamic lubrication occurs, so the coupled surfaces are totally separated by the lubricating fluid along the whole contact area;
  • The boundary lubrication is established when the high load and the slow relative motion do not allow to obtain the surfaces’ separation, so the lubrication effect is governed by chemical reactions at the contact interface; and the slow relative motion don’t allow to obtain the surfaces’ separation, so the lubrication effect is governed by chemical reactions at the contact interface;
  • In the intermediate condition, the load is supported both by asperities in contact and by the lubricating fluid pressure, so it is defined as mixed lubrication.
A particular case of the hydrodynamic mode, the Elasto-Hydrodynamic Lubrication (EHL), occurs when the lubricating fluid pressure is so high that it causes the coupled surfaces’ elastic deformation, guaranteeing the survival of a small gap between the surfaces filled by the synovial fluid even in hard loading conditions. In this case, the calculations of the lubricant’s gap within a numerical lubrication algorithm needs a deformation model able to evaluate the surfaces’ deflection due to fluid pressure [5].
In the framework of human joint replacements, the mechanical behaviour of the implant’s surfaces is determined by the coupled biomaterials and the geometry, so several solutions are available in terms of hard materials (ceramics or metal alloys) and soft materials (polymers). In the particular case of the hip arthroplasty, the implant can be geometrically modelled as a spherical cup for the acetabular component and a sphere for the femoral one.
Nowadays, an accurate description of the global tribological behaviour of the prosthesis surfaces’ interactions is an important research topic, in order to predict with an in silico approach the implant wear due to particular kinematics [8,9]: the loss of material within the prosthesis is strictly related to its duration and, as a consequence, its estimation is necessary to minimize the revision surgical procedures after the first installation.
In a total hip replacement (THR), elasto-hydrodynamic lubrication modelling several deformation models is usually adopted in the solution algorithms. When the hip implant is a type of hard-on-soft prosthesis, the constrained column model can be used in order to neglect the deformation of the harder part (generally the metallic or ceramic femoral head) and to evaluate the dominant deformation of the softer one (the acetabular cup made of polyethylene). This approach consists of the modelling of the local cup deformation as proportional to the local fluid pressure, through constants which depend on the mechanical and geometrical properties of the acetabular cup (Young modulus, Poisson ratio, inner radius and thickness) [10,11,12,13]. An approach similar to the constrained column model is to consider the analysed surface as an Elastic Foundation, in order to keep the independence of the local deformation on the surrounding pressures and obtaining a relationship useful for the contact modelling [14,15,16,17,18,19]. The deformation of other types of implants, composed of materials with comparable Young modulus, cannot be approximated by a spring model (like the constrained column one) because there is not a dominant deflection between the two surfaces and generally the local deformation could not depend only on the local pressure but also on the surroundings of the analysed points. In order to overcome these issues and to model the deformation of linear elastic materials with comparable stiffness, some authors use an equivalent spherical discrete convolution method based on the Fast Fourier Transform of the influence coefficients coming from the theory of the semi-infinite bodies in point contact, while others take advantage of several finite element simulations to build the entire matrix of influence coefficients entry by entry [20,21,22,23].
The aim of this paper is to propose, in the framework of the tribological modelling of a THR, a deformation model based on the finite element approach, able to provide the influence coefficient matrix, which evaluates the surfaces’ deformation, due to a particular pressure field acting on them and which, above all, can be written and run in the same computational environment of more complex lubrication algorithms for the accurate in silico wear assessment of the prostheses. The developed finite element model considers the two whole parts of the hip joint (acetabular cup and femoral head) and the relative motion between them, through a cubic interpolation matrix. In order to validate the model, the results are compared with the ones calculated by the software Ansys, of course for the same problem inputs. The proposal represents a succeeding improvement of previous authors’ works [2,3,4].

2. Materials and Methods

2.1. Finite Element Model

According to the classical finite elements approach, the model is based on the discretization of the analysed volumes in several linear tetrahedra with a regular mesh algorithm. In order to obtain the influence matrix, nodal forces due to the pressure are taken into account through nodal pressures acting on a face of the tetrahedron. With reference to Figure 1, the generic linear tetrahedron, composed of the nodes I , J , K and H , is subjected to the pressure action represented by the nodal pressures p I , p J and p K on the respective face with outward-pointing normal n ^ located on the face’s centre C .
Starting from the nodal pressures, the relative nodal force vector Φ p was calculated with the virtual work principle—the virtual work δ W p done by the pressure p acting on the tetrahedron’s face which produces a virtual displacement u is given in (1), considering the analysed face area A and the nodal displacement vector q .
δ W p =   p   n ^ T u   d A = Φ p T q
The pressure acting on the face is considered to be equal to the arithmetic average of the nodal pressures p m in (2).
p m = p I + p J + p K 3
By knowing the coordinates x of the tetrahedron’s nodes, the calculation of a vector n normal to the I J K face, together with the coordinates of the face centre x C , allows us to evaluate the outward-pointing normal n ^ and the area of the face A in the Equation (3).
{ n = ( x J x I ) × ( x K x I ) x C = x I + x J + x K 3 { n ^ = sign ( n T ( x H x C ) ) n | n | A = 1 2 | n |
Replacing the (3) in (1) and considering that the displacement field u is related to the nodal displacement vector q through the shape function matrix N , the virtual work is written in (4).
δ W p = p m n ^ T (   N   d A ) q = Φ p T q
Defining the vector of nodal pressure p n and considering the surface integral of the shape function matrix as the product of a matrix H with the face’s area A , Equation (5) provides the pressure nodal force vector Φ p , which is related to the nodal pressure vector p n through the matrix J p .
{ p m = [ 1 1 1 ] 3 [ p I p J p K ] = m T p n   N   d A = H A   m T p n n ^ T H A   q = Φ p T q     Φ p = ( H T n ^   m T A )   p n = J p p n
Using Equation (5) for each element j belonging to the analysed structure composed by n e elements, the total nodal force vector Φ p due to the nodal pressure vector p n acting on each node is calculated and included in the force equilibrium in Equation (6), in which the structure’s stiffness matrix K and the other nodal force vector Φ are introduced.
Φ p j = J p j p n j         Φ p = J p p n             K q = Φ + Φ p
When the displacement vector q is available, the normal displacement δ n , i of a node i is evaluated projecting its displacement q i along the mean normal direction n ^ i , which is obtained with the resultant of the outward-pointing normal vectors n ^ k i belonging to the n s pressure-loaded faces surrounding the node i . Then, the matrix J q which connects the nodal normal displacement vector δ n to the nodal displacement vector q is written in (7).
n ^ i = k = 1 n s n ^ k i | k = 1 n s n ^ k i |             δ n , i = n ^ i T q i         δ n = J q q
Partitioning the equation system in (6) by dividing the contributions of the free displacements ( l ) and the constrained ones ( v ) as follows in (8),
[ K l l K l v K v l K v v ] [ q l q v ] = [ Φ l Φ v ] + [ Φ p , l Φ p , v ]
the aim is to put in evidence the surface pressure and surface deformation fields vectors p and δ . In order to reach the goal, several matrices are introduced:
  • The matrix J n which relates the nodal pressure vector p n to the surface pressure field vector p in (9);
    p n = J n p
  • The matrix J δ which relates the surface deformation field vector δ to the nodal normal displacement vector δ n in (10); and
    δ = J δ δ n
  • The matrices J l and J v which collect the free part ( l ) or the constrained one ( v ) of a vector quantity x , useful also to rewrite it in its ordered form x o through the matrix J o in (11).
    { x l = J l x x v = J v x     x o = [ x l x v ] = [ J l J v ] x = J o x
Starting from the Equation (10), the surface deformation field vector δ is written as a function of the free displacement vector q l and the constrained one q v in (12).
δ = J δ J q J o T [ q l q v ] = [ J δ l J δ v ] [ q l q v ] = J δ l q l + J δ v q v
Then, the pressure nodal force vector related to the free displacements Φ p , l in the Equation (8) is written as a function of the pressure field vector p in (13).
Φ p , l = J l J p J n p
Solving the Equation (8) with respect to q l and replacing the latter in (12), the final Equation (14) is obtained.
δ = ( J δ l K l l 1 J l J p J n ) p + [ J δ l K l l 1 ( Φ l K l v q v ) + J δ v q v ]
In Equation (15), the influence matrix C and a vector n are highlighted; the vector n depends on the structure’s boundary conditions—the known forces acting on the free displacements Φ l and the known constrained displacements q v .
{ C = J δ l K l l 1 J l J p J n n = J δ l K l l 1 ( Φ l K l v q v ) + J δ v q v             δ = C p + n
Finally, Equation (15) can be used to evaluate the surface deformation δ due to the surface pressure p by assembling the influence matrix C and the boundary conditions vector n coming from the finite element discretization of the structure.

2.2. Acetabular Cup and Femoral Head Meshes

The regular mesh algorithm allows us to create a triangulation starting from the coordinates x of the nodes. Both for the acetabular cup and for the femoral head, the nodes coordinates, respectively x c and x h , are written in spherical coordinates in (16), considering the cup inner radius R , the cup thickness H and the head radius r . The meshes are shown in the Figure 2.
x c = [ ρ c sin θ c cos φ c ρ c sin θ c sin φ c ρ c cos θ c ] 0 θ c π 0 φ c π R ρ c R + H x h = [ ρ h sin θ h cos φ h ρ h sin θ h sin φ h ρ h cos θ h ] 0 θ h π 0 φ h 2 π 0 ρ h r
The partitioning of the displacements in terms of free or constrained ones, is the following:
  • Regarding the acetabular cup, the nodes on the outer surface are fixed because of the cup fixation with respect to the pelvic bone, namely when ρ c = R + H , while the nodes subjected to the surface pressure are located on the inner surface, namely when ρ c = R ; and
  • Regarding the femoral head, the nodes with the y -coordinate lower than a constant value depending on a chosen angle ϕ 0 are fixed because of the head fixation with respect to the femoral stem, namely when y h < r cos ϕ 0 , while the nodes subjected to the surface pressure are located on the outer surface, namely when ρ h = r .
This information allows us to calculate the influence matrix C and the boundary conditions vector n for both the structures, in order to write the Equation (15) for the two cases in the related reference frames in (17). With the chosen boundary conditions, the boundary conditions vector n is null for both cases.
δ c = C c p c + n c δ h = C h p h + n h
With the aim of implementing the deformation model in an elasto-hydrodynamic lubrication algorithm, it is worth noting that the influence matrices C c and C h are calculated once (before the lubrication algorithm solving), so the related computational expense is not involved in the lubrication algorithm iterative cycles.

2.3. Acetabular Cup and Femoral Head Coupling through Cubic Interpolation

The acetabular cup and femoral head surfaces are subjected to the same pressure, so the total deformation is given by the sum of the single contributions which refer to different grids { θ , φ } ; moreover, the two grids slide with respect to each other because of the hip relative motion between the femoral head and acetabular cup. The way to couple the grids proposed in this work consists of transforming the analysed surfaces’ quantities with a cubic interpolation.
In order to move a surface quantity f from a grid { x 0 , y 0 } to a second grid { x 1 , y 1 } , the f value on the generic point ( x , y ) of the second grid is given in (18).
f ( x , y ) = a T φ a = [ a 0 a 1 a 2 a 3 a 4 a 5 a 6 a 7 a 8 a 9 a 10 a 11 ] T φ = [ 1 x x 2 x 3 y y 2 y 3 x y x 2 y x y 2 x 3 y x y 3 ] T
With reference to the Figure 3, the coefficient vector a is evaluated solving the linear system (19) by imposing that the function f must be equal to its value on the surrounding points ( i 1 , j 1 ) , ( i 1 , j ) , ( i , j 1 ) and ( i , j ) and the same is done for its partial derivatives along the x direction, f x , and the y one, f y .
{ a T φ ( x i 1 , y j 1 ) = f i 1 , j 1   a T φ ( x i 1 , y j ) = f i 1 , j   a T φ ( x i , y j 1 ) = f i , j 1 a T φ ( x i , y j ) = f i j a T φ x ( x i 1 , y j 1 ) = f x i 1 , j 1 a T φ x ( x i 1 , y j ) = f x i 1 , j a T φ x ( x i , y j 1 ) = f x i , j 1 a T φ x ( x i , y j ) = f x i j a T φ y ( x i 1 , y j 1 ) = f y i 1 , j 1 a T φ y ( x i 1 , y j ) = f y i 1 , j a T φ y ( x i , y j 1 ) = f y i , j 1 a T φ y ( x i , y j ) = f y i j A a = b
The vector b is composed of the value of f and its partial derivative along the two directions on the surroundings points—writing the partial derivatives with the central finite difference, the vector b can be expressed as a function of the vector f s through a matrix M defined in (20).
f s = [ f i 1 , j 2 f i , j 2 f i 2 , j 1 f i 2 , j f i 1 , j 1 f i 1 , j f i , j 1 f i j f i 1 , j + 1 f i , j + 1 f i + 1 , j 1 f i + 1 , j ] T b = M f s
Since every point on the second grid can be calculated by the knowledge of the vector f s , through the definition of a matrix J 01 , s , then all the f values on the second grid, included in the vector f 1 , can be calculated through the assembled cubic interpolation matrix J 01 multiplied by the f vector values on the starting grid, f 0 , in (21).
f ( x , y ) = ( φ T A 1 M ) f s = J 01 , s f s           f 1 = J 01 f 0
In order to apply the Equation (21) to the grids belonging, respectively, to the inner surface of the acetabular cup, { θ c , φ c } , and to the outer surface of the femoral head, { θ h , φ h } , the position of the points defined in the grid of arrival has to be rotated in the reference frame of the starting grid because of the hip relative motion. With reference to the Figure 4 and denoting with R x , R y and R z the rotation matrices around the cartesian axes, the global rotation between the two grids is composed of three transformations:
  • Both for the acetabular cup and for the femoral head, the reference frame used in the finite element discretization was rotated by the inclination angle α i n and the anteversion angle β a v with respect to the anatomical reference frame (Antero/Posterior AP, Proximo/Distal PD and Medio/Lateral ML) through a rotation matrix R g defined in (22) [3]; and
    R g = R z ( π 2 β a v ) R x ( α i n )
  • The head anatomical reference frame is rotated with respect to the cup one by the Flexion/Extension angle θ F E , the Adduction/Abduction angle θ A A and the Internal/External Rotation angle θ I E R through the rotation matrix R h i p defined in (23) [3].
    R h i p = R z ( θ F E ) R x ( θ A A ) R y ( θ I E R )
Then, the cup radial unit vector defined in the x c y c z c reference frame, r ^ c , can be written in the x h y h z h reference frame as r ^ c ( h ) in (24), together with the head radial unit vector r ^ h rotated in the cup reference frame as r ^ h ( c ) .
r ^ = [ sin θ cos φ sin θ sin φ cos θ ] T r ^ c ( h ) = R g T R h i p T R g r ^ c r ^ h ( c ) = R g T R h i p R g r ^ h
Since the radial unit vector r ^ depends only on the spherical angles θ and φ , the points of the cup grid { θ c , φ c } can be rotated in the head grid obtaining the points { θ c ( h ) , φ c ( h ) } and the points of the head grid { θ h , φ h } can be rotated in the cup grid obtaining the points { θ h ( c ) , φ h ( c ) } by inverting r ^ in (24) as follows in (25).
θ = arctan ( r ^ 1 2 + r ^ 2 2 r ^ 3 ) φ = arctan ( r ^ 2 r ^ 1 )
Once the rotated grids are obtained, the cubic interpolation matrix that leads from the cup to the head J c h and the one that leads from the head to the cup J h c can be assembled. Generally, a lubrication model algorithm is referred to the cup grid, so the calculated cubic interpolation matrices can be used together with Equation (17) in order to evaluate in (26) the total deformation of both the surfaces δ due to the pressure p acting on them in the cup grid reference frame.
{ δ c = C c p c + n c δ h = C h p h + n h δ = δ c + J h c δ h p c = p p h = J c h p c         δ = ( C c + J h c C h J c h ) p + ( n c + J h c n h )
In Equation (27), the final form of the influence matrix C and of the boundary conditions vector n constituting the definitive deformation model is highlighted.
{ C = C c + J h c C h J c h n = n c + J h c n h           δ = C p + n
Finally, the whole deformation model logic is shown in the Figure 5 through block visualizations, in which the calculation of the boundary conditions vector n is not considered because of the reasons explained in the Section 2.2.

3. Results and Discussion

3.1. Validation

In order to validate Equation (17), two simulations (one for the acetabular cup and one for the femoral head) were performed in the Matlab environment (Matlab R2020b, Mathworks, Natick, Massachusetts, USA). Their results were compared with the ones computed by the software Ansys Mechanical (Ansys 2020 R2, Canonsburg, Pennsylvania, USA) for the same inputs (pressure field, mechanical and geometrical properties) but for discretization performed with quadratic tetrahedral elements instead of the linear ones used in the proposed model. The second order quadratic tetrahedral finite elements used in Ansys allows us to work with coarser mesh with respect to the one used in Matlab; moreover, the Matlab mesh is strictly connected to the interface grid on the acetabular cup surface because it will be used successively in an elasto-hydrodynamic lubrication algorithm.
The input pressure fields, characteristics of elasto-hydrodynamic lubrication shapes [3], were used in the simulations with zero pressure on the domain boundaries, both for the cup simulation and for the head one, as follows in (28).
p c ( θ c , φ c ) = ( p 0 e   ( θ c θ 0 ) 2 + ( φ c φ 0 ) 2 ( α 0 π ) 2 ) sin θ c sin φ c p h ( θ h , φ h ) = { ( p 0 e   ( θ h θ 0 ) 2 + ( φ h φ 0 ) 2 ( α 0 π ) 2 ) sin θ h sin φ h if   0 φ h π 0 if   π < φ h 2 π
The input parameters used in the calculations are listed in the Table 1. Regarding the mechanical properties, they are referred to a plastic acetabular cup (Polyethylene) and a ceramic femoral head (Alumina 88%).
The pressure fields of Equation (28) were described in Matlab code—in order to apply them within the Ansys simulations, the mesh generated by Ansys was imported in the Matlab algorithm, and the pressure field was interpolated through a cubic interpolation matrix on the Ansys mesh nodes and, after, it was imported back in the Ansys model as External Data. In Figure 6, the Ansys meshes related to the acetabular cup and the femoral head together with the surface pressure fields are shown.
The same information shown in Figure 6 is reported in Figure 7 but related to the Matlab model.
In order to compare the surfaces’ deformation, the Ansys solution results are evaluated in terms of directional deformations along the cartesian axes; then, they are projected along the radial direction so that the normal surface deformation is obtained. In Figure 8, the input surface pressure and the output surface deformation of the acetabular cup are shown over the cup grid, both for the Matlab simulation and for the Ansys one. The results show a general satisfactory agreement between the two software in terms of deformation shapes, peak coordinates and magnitude even if a slight underestimation is detected close to deformation peak—this is probably due to different finite element used in Ansys (quadratic tetrahedra) and the adopted choice to distribute uniformly the nodal pressures over the loaded face of a tetrahedron though their arithmetic average. However, also the choice of using of the linear tetrahedra allows us to perform the influence matrix assembly easily and this results in unappreciable underestimation.
In order to validate the model, the Root Mean Square Error (RMSE) and the squared correlation coefficient R 2 between the deformation results are calculated and the related data are plotted against each other in the Figure 9—the plot confirms the underestimation close to the highest values of the deformation but provides good values of the RMSE and the R 2 , respectively, 8.47   ×   10 7   m and 0.99714 .
Then, the same comparison was performed within the femoral head simulation. In Figure 10, the head surface pressure and deformation over the head grid, both in Matlab and in Ansys, are shown. With respect to the acetabular cup case, the results were even more satisfactory—a dominant underestimation or overestimation is not appreciable and the discontinuities generated by the imposed boundary conditions, in terms of null displacement in correspondence of the femoral head fixation, do not result in output instabilities.
In Figure 11, the Matlab deformation data are plotted against the Ansys one for the femoral head simulation, showing two slight deviations—an overestimation for the lowest deformation values and an underestimation for the highest ones. In this case the RMSE and the R 2 obtained are, respectively, 2.67   × 10 8   m and 0.99854 .
It is interesting to notice that, despite the fact that the influence matrix C is referred to the interface between the acetabular cup and the femoral head (in fact it relates the surface deformation δ to the surface pressure p ), it is built by rearranging the global stiffness matrix K and the nodal displacement q and force Φ vectors associated with the two analysed bodies according to the classical Finite Element theory, so other quantities referring instead to the bodies’ bulk continuum are evaluable, such as for example the Von Mises equivalent stress and the strain energy.
Then, these quantities can be calculated also with the Matlab code, so they can be compared. In Figure 12 the above quantities are shown for the cup simulation in Ansys.
The same quantities obtained by the Matlab code are shown in the Figure 13—the comparison is satisfactory in terms of magnitude order and shapes matching.
The same comparison is performed for the femoral head case. In Figure 14, the Von Mises stress and the strain energy in the Ansys environment are shown.
Then, Figure 15 reports the same outputs obtained by Matlab code. In this case, the comparison provides less satisfactory agreements regarding the strain energy—this could be due to the particular boundary conditions imposed on the fixed nodes, which lead to remarkable discontinuities around that zone, as confirmed by Figure 10 and Figure 14.

3.2. Implementation of the Deformation Model in an EHL Lubrication Algorithm

In the context of an elasto-hydrodynamic lubrication algorithm, the proposed deformation model accomplishes the task of calculating the lubricating fluid gap contribution due to the surfaces’ deformation considering also the hip relative motion. In particular, according to Equation (29) the fluid gap field h is composed of the geometrical gap, dependent on the radial clearance c and the dimensionless eccentricity vector n , added to the surfaces’ deformation δ [3,4].
h ( θ , φ ) = c ( 1 n T r ^ ( θ , φ ) ) + δ ( p ( θ , φ ) )
In order to give an example in which the cubic interpolation matrices are used, the surfaces’ deformation over the cup grid is calculated for the same prosthesis characterized by the data listed in the Table 1, but with different parameters related to the gaussian input surface pressure and with a configuration defined by a relative position between the acetabular cup and femoral head deriving from a previously algorithm developed in [3,4].
The new data are listed in the Table 2.
In Figure 16, the prosthesis configuration is shown in the anatomical reference frame of the cup in terms of surfaces’ grids, together with the fluid pressure acting on the two surfaces of the joint and the local reference frames with red x -axis, green y -axis and blue z -axis.
Once the cubic interpolation matrices are built with Equation (21) while the influence matrix and the boundary conditions vector are assembled through Equation (27), the resultant surfaces’ deformation is calculated and shown in Figure 17 together with the surface pressure.

4. Conclusions

In this work, an accurate model for the calculation of the total hip replacement surfaces deformation was proposed with the aim to merge it in EHL lubrication algorithms. The computer code was based on the finite element theory by discretizing the acetabular cup and the femoral head bodies in linear tetrahedra; the elements located on the loaded surfaces are subjected to nodal forces due to the lubricating surface pressure, evaluated by the virtual work principle. The two discretizations of the prosthesis were coupled through a cubic interpolation matrix which connects the two surfaces grids and takes into account the relative motion between the acetabular cup and the femoral head. The model was developed in order to define a global influence matrix and a boundary conditions vector, so that it can be run in the same computational context of an elasto-hydrodynamic lubrication algorithm taking advantage of its accuracy. In order to validate the proposed code, two simulations were conducted and compared with the results elaborated by the software Ansys working with quadratic tetrahedra for the same problem—in particular the surfaces’ deformation of a prosthesis made of a polyethylene acetabular cup and Alumina 88% femoral head was evaluated as a consequence of a surface pressure field modelled with a modified gaussian shape. The comparison provided a satisfactory agreement in terms of surfaces deformation, in particular recording a root mean squared error RMSE and a squared correlation coefficient R 2 equal to, respectively,   8.47 × 10 7   m and 0.99714 for the acetabular cup simulation and equal to 2.67 × 10 8   m and 0.99854 for the femoral head one. Then, the outputs regarding the Von Mises stress and the strain energy were compared resulting in a good matching in terms of magnitude order and shapes.
A third simulation was conducted to evaluate the prosthesis deformation due to a similar pressure field but with an implant configuration which considers a relative rotation between the surfaces due to the hip motion, in order to show the simplicity of implementation of the deformation model in an elasto-hydrodynamic lubrication code.
Finally, the proposed model showed quantitative and qualitative reliability, so that several insights can be conducted and the related future perspectives could regard:
  • The implementation of a routine which elaborates the contact pressure in domain zones characterized by nodes overlapping, in order to consider the model functionality within a mixed elasto-hydrodynamic lubrication algorithm;
  • The adaptation of the finite element model to consider viscoelastic materials, adding deformation contributions which depend on the time history of the pressure field through a viscosity matrix;
  • The usage of the model for other types of joint replacements such as knee, ankle, shoulders, etc.

Author Contributions

Conceptualization, A.R.; methodology, A.R. and A.S.; software, A.S.; formal analysis, A.R. and A.S.; investigation, A.R.; calculation: A.S.; discussion of the results: A.R. and A.S.; resources, A.R.; writing—original draft preparation, A.S.; writing—review and editing, A.R.; supervision, A.R.; funding acquisition, A.R. Both authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by MIUR, PRIN 2017 BIONIC.

Institutional Review Board Statement

Not Applicable.

Informed Consent Statement

Not Applicable.

Data Availability Statement

Not Applicable.

Conflicts of Interest

The authors declare no conflict of interest.

Nomenclature

I ,   J ,   K ,   H Linear tetrahedron nodes
n ^ Tetrahedron face outward-pointing normal
A Tetrahedron face area
N Shape function matrix
p I ,   p J ,   p K Nodal pressures
p n Nodal pressure vector
Φ Nodal force vector
q Nodal displacement vector
K Stiffness matrix
δ n Nodal normal displacement vector
p Surface pressure field vector
δ Surface deformation field vector
C Influence matrix
n Boundary conditions vector or dimensionless eccentricity vector
x Cartesian coordinates vector
ρ ,   θ ,   φ Spherical coordinates
R Acetabular cup inner radius
H Acetabular cup thickness
r Femoral head radius
ϕ 0 Femoral head fixed nodes angle
i , j Surface grid point indices
J 01 Cubic interpolation matrix from the grid 0 to the grid 1
R Rotation matrix
α i n ,   β a v Inclination and anteversion angles
θ F E ,   θ A A ,   θ I E R Hip Flexion/Extension, Adduction/Abduction and Internal/External Rotation angles
r ^ Radial unit vector
E ,   ν Young modulus and Poisson ratio
RMSE ,   R 2 Root Mean Squared Error and squared correlation coefficient
c Radial clearance
h Lubricating fluid gap

References

  1. Fisher, J.; Dowson, D. Tribology of Total Artificial Joints. Proc. Inst. Mech. Eng. Part H J. Eng. Med. 1991, 205, 73–79. [Google Scholar] [CrossRef]
  2. Ruggiero, A.; Sicilia, A. Lubrication modeling and wear calculation in artificial hip joint during the gait. Tribol. Int. 2020, 142, 105993. [Google Scholar] [CrossRef]
  3. Ruggiero, A.; Sicilia, A. A Mixed Elasto-Hydrodynamic Lubrication Model for Wear Calculation in Artificial Hip Joints. Lubricants 2020, 8, 72. [Google Scholar] [CrossRef]
  4. Ruggiero, A.; Sicilia, A.; Affatato, S. In silico total hip replacement wear testing in the framework of ISO 14242-3 accounting for mixed elasto-hydrodynamic lubrication effects. Wear 2020, 460, 203420. [Google Scholar] [CrossRef]
  5. Mattei, L.; Di Puccio, F.; Piccigallo, B.; Ciulli, E. Lubrication and wear modelling of artificial hip joints: A review. Tribol. Int. 2011, 44, 532–549. [Google Scholar] [CrossRef]
  6. Gleghorn, J.P.; Bonassar, L.J. Lubrication mode analysis of articular cartilage using Stribeck surfaces. J. Biomech. 2008, 41, 1910–1918. [Google Scholar] [CrossRef]
  7. Scholes, S.C.; Unsworth, A. Comparison of friction and lubrication of different hip prostheses. Proc. Inst. Mech. Eng. Part H J. Eng. Med. 2000, 214, 49–57. [Google Scholar] [CrossRef]
  8. Ruggiero, A.; Merola, M.; Affatato, S. Finite element simulations of hard-on-soft hip joint prosthesis accounting for dynamic loads calculated from a musculoskeletal model during walking. Materials 2018, 11, 574. [Google Scholar] [CrossRef] [Green Version]
  9. Affatato, S.; Ruggiero, A. A Perspective on Biotribology in Arthroplasty: From In Vitro toward the Accurate In Silico Wear Prediction. Appl. Sci. 2020, 10, 6312. [Google Scholar] [CrossRef]
  10. Jalali-Vahid, D.; Jagatia, M.; Jin, Z.M.; Dowson, D. Prediction of lubricating film thickness in a ball-in-socket model with a soft lining representing human natural and artificial hip joints. Proc. Inst. Mech. Eng. Part J J. Eng. Tribol. 2001, 215, 363–372. [Google Scholar] [CrossRef]
  11. Jalali-Vahid, D.; Jin, Z.M.; Dowson, D. Elastohydrodynamic lubrication analysis of hip implants with ultra high molecular weight polyethylene cups under transient conditions. Proc. Inst. Mech. Eng. Part C J. Mech. Eng. Sci. 2003, 217, 767–777. [Google Scholar] [CrossRef]
  12. Jalali-Vahid, D.; Jagatia, M.; Jin, Z.; Dowson, D. Prediction of lubricating film thickness in UHMWPE hip joint replace-ments. J. Biomech. 2001, 34, 261–266. [Google Scholar] [CrossRef]
  13. Jalali-Vahid, D.; Jin, Z.M. Transient elastohydrodynamic lubrication analysis of ultra-high molecular weight polyethylene hip joint replacements. Proc. Inst. Mech. Eng. Part C J. Mech. Eng. Sci. 2001, 216, 409–420. [Google Scholar] [CrossRef]
  14. Fregly, B.J.; Bei, Y.; Sylvester, M.E. Experimental evaluation of an elastic foundation model to predict contact pressures in knee replacements. J. Biomech. 2003, 36, 1659–1668. [Google Scholar] [CrossRef]
  15. Halloran, J.P.; Easley, S.K.; Petrella, A.J.; Rullkoetter, P.J. Comparison of deformable and elastic foundation finite element simulations for predicting knee replacement mechanics. J. Biomech. Eng. 2005, 127, 813–818. [Google Scholar] [CrossRef]
  16. Askari, E.; Andersen, M.S. A closed-form formulation for the conformal articulation of metal-on-polyethylene hip prostheses: Contact mechanics and sliding distance. Proc. Inst. Mech. Eng. Part H J. Eng. Med. 2018, 232, 1196–1208. [Google Scholar] [CrossRef] [Green Version]
  17. Pérez-González, A.; Fenollosa-Esteve, C.; Sancho-Bru, J.L.; Sánchez-Marín, F.T.; Vergara, M.; Rodríguez-Cervantes, P.J. A modified elastic foundation contact model for application in 3D models of the prosthetic knee. Med. Eng. Phys. 2008, 30, 387–398. [Google Scholar] [CrossRef]
  18. Mukras, S.M.; Kim, N.H.; Schmitz, T.L.; Sawyer, W.G. Evaluation of contact force and elastic foundation models for wear analysis of multibody systems. In Proceedings of the International Design Engineering Technical Conferences and Computers and Information in Engineering Conference, Montreal, QC, Canada, 15–18 August 2010; Volume 44106, pp. 1565–1577. [Google Scholar]
  19. Srivastava, G.; Christian, N.; Higgs, C.F., III. A predictive framework of the tribological impact of physical activities on metal-on-plastic hip implants. Biotribology 2021, 25, 100156. [Google Scholar] [CrossRef]
  20. Wang, F.; Jin, Z. Prediction of elastic deformation of acetabular cups and femoral heads for lubrication analysis of artificial hip joints. Proc. Inst. Mech. Eng. Part J J. Eng. Tribol. 2004, 218, 201–209. [Google Scholar] [CrossRef]
  21. Evans, H.P.; Hughes, T.G. Evaluation of deflection in semi-infinite bodies by a differential method. Proc. Inst. Mech. Eng. Part C J. Mech. Eng. Sci. 2000, 214, 563–584. [Google Scholar] [CrossRef]
  22. Gao, L.; Dowson, D.; Hewson, R.W. Predictive wear modeling of the articulating metal-on-metal hip replacements. J. Biomed. Mater. Res. Part B Appl. Biomater. 2015, 105, 497–506. [Google Scholar] [CrossRef]
  23. Jagatia, M.; Jin, Z. Elastohydrodynamic lubrication analysis of metal-on-metal hip prostheses under steady state entraining motion. Proc. Inst. Mech. Eng. Part H J. Eng. Med. 2001, 215, 531–541. [Google Scholar] [CrossRef]
Figure 1. Linear tetrahedron with I, J, K, and H nodes subjected to nodal pressures acting on the face IJK.
Figure 1. Linear tetrahedron with I, J, K, and H nodes subjected to nodal pressures acting on the face IJK.
Lubricants 09 00041 g001
Figure 2. Acetabular cup (a) and femoral head (b) meshes in the related reference frames.
Figure 2. Acetabular cup (a) and femoral head (b) meshes in the related reference frames.
Lubricants 09 00041 g002
Figure 3. A 2D interpolation scheme with the points surrounding the desired evaluation.
Figure 3. A 2D interpolation scheme with the points surrounding the desired evaluation.
Lubricants 09 00041 g003
Figure 4. Reference frames of the hip joint.
Figure 4. Reference frames of the hip joint.
Lubricants 09 00041 g004
Figure 5. Deformation model blocks logic.
Figure 5. Deformation model blocks logic.
Lubricants 09 00041 g005
Figure 6. Ansys meshes and pressure loaded surfaces—acetabular cup (a) and femoral head (b).
Figure 6. Ansys meshes and pressure loaded surfaces—acetabular cup (a) and femoral head (b).
Lubricants 09 00041 g006
Figure 7. Matlab meshes and pressure loaded surfaces—acetabular cup (a) and femoral head (b).
Figure 7. Matlab meshes and pressure loaded surfaces—acetabular cup (a) and femoral head (b).
Lubricants 09 00041 g007
Figure 8. Acetabular cup surface pressure input in Matlab (surf) and Ansys (black circles) (a); Matlab deformation output (surf) and Ansys deformation output (black circles) related to the acetabular cup (b).
Figure 8. Acetabular cup surface pressure input in Matlab (surf) and Ansys (black circles) (a); Matlab deformation output (surf) and Ansys deformation output (black circles) related to the acetabular cup (b).
Lubricants 09 00041 g008
Figure 9. Matlab deformation results data against Ansys ones regarding the acetabular cup simulation.
Figure 9. Matlab deformation results data against Ansys ones regarding the acetabular cup simulation.
Lubricants 09 00041 g009
Figure 10. Femoral head surface pressure input in Matlab (surf) and Ansys (black circles) (a); Matlab deformation output (surf) and Ansys deformation output (black circles) related to the femoral head (b).
Figure 10. Femoral head surface pressure input in Matlab (surf) and Ansys (black circles) (a); Matlab deformation output (surf) and Ansys deformation output (black circles) related to the femoral head (b).
Lubricants 09 00041 g010
Figure 11. Matlab deformation results data against Ansys ones regarding the femoral head simulation.
Figure 11. Matlab deformation results data against Ansys ones regarding the femoral head simulation.
Lubricants 09 00041 g011
Figure 12. Von Mises stress (a) and strain energy (b) within the acetabular cup simulation in Ansys.
Figure 12. Von Mises stress (a) and strain energy (b) within the acetabular cup simulation in Ansys.
Lubricants 09 00041 g012
Figure 13. Von Mises stress (a) and strain energy (b) within the acetabular cup simulation in Matlab.
Figure 13. Von Mises stress (a) and strain energy (b) within the acetabular cup simulation in Matlab.
Lubricants 09 00041 g013
Figure 14. Von Mises stress (a) and strain energy (b) within the femoral head simulation in Ansys.
Figure 14. Von Mises stress (a) and strain energy (b) within the femoral head simulation in Ansys.
Lubricants 09 00041 g014
Figure 15. Von Mises stress (a) and strain energy (b) within the femoral head simulation in Matlab.
Figure 15. Von Mises stress (a) and strain energy (b) within the femoral head simulation in Matlab.
Lubricants 09 00041 g015
Figure 16. Prosthesis configuration due the chosen relative rotations.
Figure 16. Prosthesis configuration due the chosen relative rotations.
Lubricants 09 00041 g016
Figure 17. Fluid pressure input (a) and total surface deformation results obtained (b).
Figure 17. Fluid pressure input (a) and total surface deformation results obtained (b).
Lubricants 09 00041 g017
Table 1. Input parameters related to the comparative simulations.
Table 1. Input parameters related to the comparative simulations.
ParameterValue
Femoral head radius r14 mm
Acetabular cup inner radius R 14.03 mm
Acetabular cup thickness H 9 mm
Femoral head fixed nodes angle ϕ 0 30°
Acetabular cup Young modulus E c 1.1 GPa
Acetabular cup Poisson ratio ν c 0.42
Femoral head Young modulus E h 245.9 GPa
Femoral head Poisson ratio ν h 0.24
Pressure gaussian peak p 0 107 Pa
Pressure gaussian θ-translation θ 0 2π/3 rad
Pressure gaussian φ-translation φ 0 π/3 rad
Pressure gaussian dimensionless width α 0 0.2
Table 2. Input parameters related to the simulation with relative rotations between the acetabular cup and the femoral head.
Table 2. Input parameters related to the simulation with relative rotations between the acetabular cup and the femoral head.
ParameterValue
Pressure gaussian peak p 0 107 Pa
Pressure gaussian θ-translation θ 0 π/4 rad
Pressure Gaussian φ-translation φ 0 π/4 rad
Pressure gaussian dimensionless width α 0 0.2
Anteversion angle β a v 90°
Inclination angle α i n 45°
Flexion/Extension angle θ F E −40°
Adduction/Abduction angle θ A A 10°
Internal/External Rotation angle θ I E R
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Ruggiero, A.; Sicilia, A. Mathematical Development of a Novel Discrete Hip Deformation Algorithm for the In Silico Elasto-Hydrodynamic Lubrication Modelling of Total Hip Replacements. Lubricants 2021, 9, 41. https://doi.org/10.3390/lubricants9040041

AMA Style

Ruggiero A, Sicilia A. Mathematical Development of a Novel Discrete Hip Deformation Algorithm for the In Silico Elasto-Hydrodynamic Lubrication Modelling of Total Hip Replacements. Lubricants. 2021; 9(4):41. https://doi.org/10.3390/lubricants9040041

Chicago/Turabian Style

Ruggiero, Alessandro, and Alessandro Sicilia. 2021. "Mathematical Development of a Novel Discrete Hip Deformation Algorithm for the In Silico Elasto-Hydrodynamic Lubrication Modelling of Total Hip Replacements" Lubricants 9, no. 4: 41. https://doi.org/10.3390/lubricants9040041

APA Style

Ruggiero, A., & Sicilia, A. (2021). Mathematical Development of a Novel Discrete Hip Deformation Algorithm for the In Silico Elasto-Hydrodynamic Lubrication Modelling of Total Hip Replacements. Lubricants, 9(4), 41. https://doi.org/10.3390/lubricants9040041

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop