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

: 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 deﬁnition of the inﬂuence 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.


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].

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 normaln located on the face's centre C.

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 , , and , is subjected to the pressure action represented by the nodal pressures , and on the respective face with outward-pointing normal located on the face's centre . Starting from the nodal pressures, the relative nodal force vector was calculated with the virtual work principle-the virtual work done by the pressure acting on the tetrahedron's face which produces a virtual displacement is given in (1), considering the analysed face area and the nodal displacement vector .
The pressure acting on the face is considered to be equal to the arithmetic average of the nodal pressures in (2).
By knowing the coordinates of the tetrahedron's nodes, the calculation of a vector normal to the face, together with the coordinates of the face centre , allows us to evaluate the outward-pointing normal and the area of the face in the Equation (3).
Replacing the (3) in (1) and considering that the displacement field is related to the nodal displacement vector through the shape function matrix , the virtual work is written in (4). 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.
The pressure acting on the face is considered to be equal to the arithmetic average of the nodal pressures p m in (2).
By knowing the coordinates x of the tetrahedron's nodes, the calculation of a vector n normal to the I JK face, together with the coordinates of the face centre x C , allows us to evaluate the outward-pointing normaln and the area of the face A in the Equation (3).
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).
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 .
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.
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 directionn i , which is obtained with the resultant of the outward-pointing normal vectorsn ki 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).
Partitioning the equation system in (6) by dividing the contributions of the free displacements (l) and the constrained ones (v) as follows in (8), 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 (9) • The matrix J δ which relates the surface deformation field vector δ to the nodal normal displacement vector δ n in (10); and • 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).
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).
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).
Solving the Equation (8) with respect to q l and replacing the latter in (12), the final Equation (14) is obtained.
Lubricants 2021, 9, 41 5 of 19 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 .
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.

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.
Solving the Equation (8) with respect to and replacing the latter in (12), the final Equation (14) is obtained.
In Equation (15), the influence matrix and a vector are highlighted; the vector depends on the structure's boundary conditions-the known forces acting on the free displacements and the known constrained displacements .
Finally, Equation (15) can be used to evaluate the surface deformation due to the surface pressure by assembling the influence matrix and the boundary conditions vector coming from the finite element discretization of the structure.

Acetabular Cup and Femoral Head Meshes
The regular mesh algorithm allows us to create a triangulation starting from the coordinates of the nodes. Both for the acetabular cup and for the femoral head, the nodes coordinates, respectively and , are written in spherical coordinates in (16)   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 = + , while the 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.
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.

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).
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 .

ditions vector
for both the structures, in order to write the Equation (1 cases in the related reference frames in (17). With the chosen boundary c boundary conditions vector is null for both cases.
With the aim of implementing the deformation model in an elasto-h lubrication algorithm, it is worth noting that the influence matrices and lated once (before the lubrication algorithm solving), so the related comp pense is not involved in the lubrication algorithm iterative cycles.

Acetabular Cup and Femoral Head Coupling through Cubic Interpolation
The acetabular cup and femoral head surfaces are subjected to the sam the total deformation is given by the sum of the single contributions which ent grids , ; moreover, the two grids slide with respect to each other b hip relative motion between the femoral head and acetabular cup. The way grids proposed in this work consists of transforming the analysed surfac with a cubic interpolation.
In order to move a surface quantity from a grid , to a second the value on the generic point , of the second grid is given in (18).
With reference to the Figure 3, the coefficient vector is evaluated solv system (19) by imposing that the function must be equal to its value on the points − 1, − 1 , − 1, , , − 1 and , and the same is done for rivatives along the direction, , and the one, .
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).
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).
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 α in and the anteversion angle β av 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 • The head anatomical reference frame is rotated with respect to the cup one by the Flexion/Extension angle θ FE , the Adduction/Abduction angle θ AA and the Internal/External Rotation angle θ IER through the rotation matrix R hip defined in (23) [3].  Then, the cup radial unit vector defined in the reference frame, , can be written in the reference frame as in (24), together with the head radial unit vector rotated in the cup reference frame as .
Once the rotated grids are obtained, the cubic interpolation matrix that leads from the cup to the head and the one that leads from the head to the cup 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 acting on them in the cup grid reference frame.
In Equation (27), the final form of the influence matrix and of the boundary conditions vector constituting the definitive deformation model is highlighted.
Finally, the whole deformation model logic is shown in the Figure 5 through block visualizations, in which the calculation of the boundary conditions vector is not considered because of the reasons explained in the Section 2.2.
Since the radial unit vectorr 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 θ Once the rotated grids are obtained, the cubic interpolation matrix that leads from the cup to the head J ch and the one that leads from the head to the cup J hc 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.
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. 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.

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).
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%).

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, PA, 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).
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 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.  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. The same information shown in Figure 6 is reported in Figure 7 but related to the Matlab model.  Lubricants 2021, 9, x 11 of 19 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 between the deformation results are calculated and the related data are plotted against each other in the Figure 9  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.  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. 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.
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 obtained are, respectively, 2.67 × 10 m and 0.99854. It is interesting to notice that, despite the fact that the influence matrix is referred to the interface between the acetabular cup and the femoral head (in fact it relates the surface deformation to the surface pressure ), it is built by rearranging the global stiffness matrix and the nodal displacement 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.  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.
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 obtained are, respectively, 2.67 × 10 m and 0.99854. It is interesting to notice that, despite the fact that the influence matrix is referred to the interface between the acetabular cup and the femoral head (in fact it relates the surface deformation to the surface pressure ), it is built by rearranging the global stiffness matrix and the nodal displacement 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 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 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. 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 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 Figures 10 and 14.
Lubricants 2021, 9, x 15 of 19 be due to the particular boundary conditions imposed on the fixed nodes, which lead to remarkable discontinuities around that zone, as confirmed by Figures 10 and 14.

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 ℎ is composed of the geometrical gap, dependent on the radial clearance and the dimensionless eccentricity vector , added to the surfaces' deformation [3,4].
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.

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].
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. Table 2. Input parameters related to the simulation with relative rotations between the acetabular cup and the femoral head.

Parameter Value
Pressure gaussian peak p 0 10 7 Pa Pressure gaussian θ-translation θ 0 π/4 rad Pressure Gaussian ϕ-translation ϕ 0 π/4 rad Pressure gaussian dimensionless width α 0 0.2 Anteversion angle β av 90 • Inclination angle α in 45 • Flexion/Extension angle θ FE −40 • Adduction/Abduction angle θ AA 10 • Internal/External Rotation angle θ IER 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.

Adduction/Abduction angle 10°
Internal/External Rotation angle 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 -axis, green -axis and blue -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.  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.

Adduction/Abduction angle 10°
Internal/External Rotation angle 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 -axis, green -axis and blue -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.

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.

Conflicts of Interest:
The authors declare no conflict of interest.