A Mixed Elasto-Hydrodynamic Lubrication Model for Wear Calculation in Artificial Hip Joints

The aim of this paper was to propose a novel in silico mixed elasto-hydrodynamic lubrication model with the purpose of wear prediction in Total Hip Replacements (THRs). The model considers the progressive wear contribution in the calculation of the meatus filled by the non-Newtonian synovial fluid. The results were referred to the gait cycle kinematics, calculated by using musculoskeletal multibody software, while the loading was assumed by literature in vivo measurements. The simulations allow evaluating the fluid and the contact pressure fields and the acetabular cup wear over the time. The results were obtained considering a Ultra High Molecular Weight PolyEthylene, UHMWPE, cup and were compared with results from the literature, showing a good agreement in terms of total volume wear of the cup.


Introduction
Nowadays many people suffer of osteoarthritis, a joint disease caused by several factors, such as aging, trauma or intense sport activities. The disturbance increases when the joint surfaces slide excessively on each other, causing cartilage deterioration and, consequently, direct contact between the bones, which could produce deformation, wear and pain.
A solution currently adopted is the replacement of the unhealthy joint with a prosthesis, an implant made of biomaterials designed to improve stability, load capacity and mobility, and to guarantee minimal friction and wear.
Examples of actual implants are represented by hip and knee artificial joints, which are the main human joints replaced: several surgical strategies and techniques are adopted in terms of total or partial replacements [1,2].
The implants can be viewed as complex tribological systems and for this reason the implant activity produces wear. The duration is strictly connected to the wear rate of the joint, hence researchers are currently focused on measuring or estimating it through experiments (in vitro), mainly by using simulators. Nowadays, the scientific interest is rising around the possibility to obtain accurate theoretical models in order to move toward the in silico approach for the analysis of the prosthesis behavior in terms of wear [3][4][5][6][7].
Since the in vitro approach requires so much time and resources, even if it is a very accurate modality in the framework of the phenomena prediction, a deep knowledge of the bio-tribology of the implants is needed, so that numerical models could be built to simulate in a more and more accurate way, in general and in a faster way with respect to the in vitro approach in particular, the best combinations of materials and geometry in correspondence of generic loading conditions of the Some techniques are based on the reduced Reynolds equation, FEM models or asperity contact pressure models, and they are adopted by authors and available in the scientific literature [6,[18][19][20][21].
In [22], an EHL simulation model was developed to analyze the tribological behavior of a metal-on-metal hip implant subjected to steady-state and physiological loading and motion gait cycle conditions by using the multi-grid method and the fast Fourier transform. The authors investigated the contact area position within the cup and the distribution of the fluid film thickness in dependence on the load components.
In [23], a full numerical analysis of the hydrodynamic lubrication problem was reported regarding artificial hip joints made of hard materials, such as metal-on-metal or ceramic-on-ceramic, through the Reynolds equation solved during walking conditions, in order to study the effects of the design parameters like radii and clearance on the magnitude of the predicted film thickness and considering a machined dimple on the acetabular cup in several positions.
Askari and Andersen, in [24], studied the hydrodynamic lubrication problem of hip implants during gait, and developed a computational model based on multibody dynamics and the Reynolds equation, considering the translational and rotational relative motion between femoral head and acetabular cup, and investigating the main differences with respect to the case in which only rotational motion is considered.
Jalali-Vahid et al., in [25], conducted a wide parametric analysis of the EHL problem of hip joints made of a UHMWPE acetabular cup against metallic or ceramic femoral heads by solving the Reynolds equation.
In [26], the authors considered the cyclic load and speed walking conditions in an EHL analysis for the artificial hip joint replacements made of a UHMWPE cup against a hard femoral head, solving the non-linear Reynolds equation and the elasticity equation through the Newton-Raphson method and finding that the predicted minimum fluid film thickness stayed constant despite large angular velocity and load changes.
In [27], an interesting EHL analysis was reported for metal-on-metal hip-resurfacing prostheses under simple steady-state rotation, by solving the Reynolds equation with the finite difference method and the elasticity equation with a finite element model, comparing the predicted minimum fluid film thickness with the one calculated through the Hamrock and Dowson formula.
In [28], a general transient EHL model was developed for metal-on-metal artificial hip joints during walking conditions using an equivalent discrete spherical convolution model and the corresponding fast Fourier transform to evaluate the elastic deformation of the spherical bearing surfaces, investigating the effects of the cup inclination angle and the lubricant viscosity on the prosthesis lubrication.
The aim of this paper is to build a novel numerical model able to solve a modified Reynolds equation for a soft-on-hard artificial hip joint with a UHMWPE acetabular cup against a hard femoral head (ceramic or metallic) in the case of a MEHL mode, taking into account: - The progressive wear phenomenon in the lubricating gap calculation through a modified Archard law; - The non-Newtonian synovial fluid behavior through the cross-viscosity model; - The cup surface deflection using a constrained column model.
In this paper, the model was used in the case of normal gait cycle, but it is adaptable to any hip kinematics and loads conditions.

The Hip Joint during the Gait
The hip joint can be viewed as a spherical joint that links the pelvis and femur bodies, so it allows three relative rotations defined in three axes: The Flexion-Extension rotation (FE) around the Medio-Lateral axis (ML) perpendicular to the sagittal plane; - The Adduction-Abduction rotation (AA) around the Anterior-Posterior axis (AP) perpendicular to the frontal plane; - The Internal-External Rotation (IER) around the Proximo-Distal axis (PD) perpendicular to the horizontal plane.
The lubrication analysis is performed in a reference frame xyz fixed at the cup center, with the x and z axes lying at the cup edges and the y axis lying on the cup axis to form a right-hand reference frame as shown in Figure 1 together with the three hip axes. The lubrication analysis is performed in a reference frame fixed at the cup center, with the and axes lying at the cup edges and the axis lying on the cup axis to form a right-hand reference frame as shown in Figure 1 together with the three hip axes. During a certain relative motion, the hip is subjected to a contact force that is composed by its three projections on the three hip axes.
Once the relative angular velocity vector and the contact force vector time data in the hip joint reference frame are available [5] and defined in (1), it is necessary to rotate them in order to study the lubrication in the cup reference frame.
The transformation is obtained in (2) by the rotation matrix , which takes into account the inclination angle , defined as the angle formed by the cup axis and the sagittal plane, and the anteversion angle , defined as the angle between the cup axis projection on the sagittal plane and the AP axis. During a certain relative motion, the hip is subjected to a contact force that is composed by its three projections on the three hip axes.
Once the relative angular velocity vector ω hip and the contact force vector N hip time data in the hip joint reference frame are available [5] and defined in (1), it is necessary to rotate them in order to study the lubrication in the cup reference frame.
The transformation is obtained in (2) by the rotation matrix R c , which takes into account the inclination angle α in , defined as the angle formed by the cup axis and the sagittal plane, and the anteversion angle β av , defined as the angle between the cup axis projection on the sagittal plane and the AP axis.
Finally, the vectors N and ω defined in the cup reference frame are obtained in (3).
The gait cycle is composed by the stance phase, starting when the heel strikes the ground and ending when the toes leaves it, and by the swing phase, until the heel strikes the ground again [29].
Since the interest is to perform a tribological analysis during walking, the time data about the angular velocity vector and the contact force vector during the gait cycle are needed.
OpenSim ® software was chosen for the simulations allowing the calculation of the three hip angles evolution during the gait cycle and the angular velocities, which are directly involved in the lubrication analysis.
Regarding the contact force, the experimental measurements done by Bergmann through instrumented implants [30] were used, considering them more realistic than the simulated ones.
Both the angular velocity vector and the Bergmann loads are referred to a single normal gait cycle with a period of 1.1 s and to an average man with average body proportions, as shown in Figure 2. The gait cycle is composed by the stance phase, starting when the heel strikes the ground and ending when the toes leaves it, and by the swing phase, until the heel strikes the ground again [29].
Since the interest is to perform a tribological analysis during walking, the time data about the angular velocity vector and the contact force vector during the gait cycle are needed.
OpenSim® software was chosen for the simulations allowing the calculation of the three hip angles evolution during the gait cycle and the angular velocities, which are directly involved in the lubrication analysis.
Regarding the contact force, the experimental measurements done by Bergmann through instrumented implants [30] were used, considering them more realistic than the simulated ones.
Both the angular velocity vector and the Bergmann loads are referred to a single normal gait cycle with a period of 1.1 s and to an average man with average body proportions, as shown in Figure  2. After the rotation, the definitive load and motion input needed for the lubrication analysis are shown in Figure 3. After the rotation, the definitive load and motion input needed for the lubrication analysis are shown in Figure 3. After the rotation, the definitive load and motion input needed for the lubrication analysis are shown in Figure 3.

The Modified Reynolds Lubrication Equation
Under the classical hypothesis and referring to the Figure 4, the dynamical equilibrium and the mass one will lead to the Reynolds equation.

The Modified Reynolds Lubrication Equation
Under the classical hypothesis and referring to the Figure 4, the dynamical equilibrium and the mass one will lead to the Reynolds equation. The general Reynolds lubrication equation for a compressible fluid is given in (4), in which the fluid pressure depends on some parameters and inputs, such as the sliding effect and the squeeze one on the right-hand side of the equation. The equation is integrated over a domain Ω closed by the boundary pressure .
The dependence of the fluid density on the fluid pressure [31] is written in (5) through the Dowson-Higginson relationship. The density-pressure relation is typical of mechanical applications, while the synovial fluid is generally considered uncompressible: at any rate the introduction of this dependence in the model leads to taking into account the eventuality of this particular behavior, which could be The general Reynolds lubrication equation for a compressible fluid is given in (4), in which the fluid pressure p depends on some parameters and inputs, such as the sliding effect and the squeeze one on the right-hand side of the equation. The equation is integrated over a domain Ω closed by the boundary pressure p 0 .
The dependence of the fluid density ρ on the fluid pressure [31] is written in (5) through the Dowson-Higginson relationship. The density-pressure relation is typical of mechanical applications, while the synovial fluid is generally considered uncompressible: at any rate the introduction of this dependence in the model leads to taking into account the eventuality of this particular behavior, which could be considered in case of faster physical activities; moreover, even if the density keeps a constant value with a density-pressure relationship it will not affect the Reynolds equation nature and then its solution.
The fluid viscosity µ is based on the Barus model [31], as in (6), with a modification that considers an effective nominal viscosity µ e f f in order to take into account the non-Newtonian behaviour of the synovial fluid. Despite the fact that the Barus model was introduced for mechanical bearings application, is interesting to implement it in the lubrication model because it can significantly affect the synovial fluid viscosity during faster human physical activities (jumping or running): during the gait cycle and for a soft-on-hard hip prosthesis the pressure levels reached in the gap are not able to make the viscosity variation appreciable with respect to the pressure one.
The effective viscosity, in (7), depends on the average shear rate γ, that is the ratio between the sliding velocity v sl and the gap h through a Cross model, in order to consider the viscosity increase due the synovial fluid protein aggregation effect in correspondence with low sliding velocity [10,32]. The introduction of the Cross model allows us to examine a wide range of biomechanical activities, because it considers the synovial fluid shear thinning effect at low sliding velocity, up to a viscosity µ 0 , and, on the other hand, returns a constant asymptotic viscosity µ ∞ in correspondence of very high sliding velocity, through the parameters k and n.
The gap h is composed, in (8), of the geometrical gap h g , depending on the particular pair undeformed geometry and configuration, the surfaces' elastic deformation δ and the linear wear u w cumulated instant by instant. The geometrical gap h g together with the entraining velocity v are characteristic of the particular joint, and their determination modifies the Reynolds equation adapting it to the particular geometry and relative motion.
The deformation δ, in (9), is composed of the fluid pressure deformation δ f , obtained by a deformation model D, and the contact deformation δ c , which is calculated by avoiding the overlap between the surfaces. The overlap is considered if the fluid film thickness decreases below a boundary layer limit ∆ b , which depends on the protein layer adsorbed by the surfaces [10]. In particular, if, in a specific zone of the domain, the fluid pressure p produces a surface deformation δ f not big enough to guarantee the complete surfaces' separation, the resulting overlapping is used to evaluate another amount of deformation δ c , namely the contact deformation, which rises in order to avoid the surfaces' interpenetration. The calculation of the contact deformation leads to the determination of the contact pressure p c , in (10), with the same deformation model inversion.
The deformation model adopted is the constrained column one, in (11), and it is justified, in case of soft-on-hard pairs, by the low cup Young modulus with respect to the head one, which is considered rigid, and it consists of a proportional dependence between the deformation and the pressure through a spring constant k d defined by the cup geometrical and mechanical properties. Despite the polymeric mechanical behavior being visco-elastic, the constrained column model provides a good approximation of the soft acetabular cup deformation in this framework [15,17,25,26].
Then, the domain will be divided into lubricated areas and contact ones, so that the total pressure p t , expressed as the sum of the two contributions in (12), will support the load on the joint.
Moreover, the total pressure produces wear, together with the sliding velocity, through a modified Archard model in which the wear rate τ w is calculated, in (13), with a wear factor k w that, in dependence on the λ ratio (the ratio between the gap h and the combined average roughness R a ), modifies the intensity of the nominal wear factor k w0 , so that the wear intensity takes into account the contact severity [33][34][35]. In particular, the wear coefficient function k w is expressed as a negative power α w of the lambda ratio, resulting in severe wear in those zones characterized by low film thickness, down to the boundary layer thickness (contact areas, λ = ∆ b /R a ), and in negligible wear in the other zones characterised by film thickness that is much larger than the average roughness (lubricated areas, λ > 3).
The wear rate obtained is integrated, in (14), instant-by-instant over the time t to calculate the contribution of the linear wear u w to the gap.
Once the total pressure p t is obtained, it is integrated into the domain's area to find the load N, and the same is done with the linear wear u w to evaluate the wear volume V w , as written in (15).

The Spherical Joint
Regarding a spherical joint, that is the hip case, the Reynolds equation is written in the cup reference frame xyz, shown in Figure 5, and fixed to the cup center O c with axes x and z lying on the cup edges and axis y lying on the cup axis to form a right-hand reference frame. Then, the equation is written in spherical coordinates through the radial unit vector as a function of the chosen spherical angles and written in (16).
The Reynolds equation written in spherical coordinates, coupled with the boundary pressure on the cup edges, leads to the spherical Reynolds problem in (17).
The geometrical gap ℎ is calculated in (18) with the dimensionless eccentricity vector and the radial clearance .
The spring constant used for the constrained column model depends on the cup radius , the cup thickness , the cup Young modulus and Poisson ratio and is shown in (19).
Thus the fluid film thickness in the Reynolds equation is obtained in (20).
Since the cup is fixed with the reference frame, the velocity of a point P located by on the cup is null, while the velocity of a point Q on the head is composed of a translational part, due to the time variation of the eccentricity , and a rotational part, due to the head angular velocity vector [24]. The two velocity vectors are defined in (21).
Thus the entraining velocity vector is calculated as the arithmetic average of the P and Q velocities, while the sliding velocity vector is calculated as the difference between the P and Q velocities, as written in (22). Then, the equation is written in spherical coordinates through the radial unit vectorr as a function of the chosen spherical angles θ and ϕ written in (16).
The Reynolds equation written in spherical coordinates, coupled with the boundary pressure p 0 on the cup edges, leads to the spherical Reynolds problem in (17).
The geometrical gap h g is calculated in (18) with the dimensionless eccentricity vector n and the radial clearance c.
The spring constant k d used for the constrained column model depends on the cup radius R, the cup thickness H, the cup Young modulus E and Poisson ratio ν and is shown in (19).
Thus the fluid film thickness in the Reynolds equation is obtained in (20).
Since the cup is fixed with the reference frame, the velocity v P of a point P located byr on the cup is null, while the velocity v Q of a point Q on the head is composed of a translational part, due to the time variation of the eccentricity e, and a rotational part, due to the head angular velocity vector ω [24]. The two velocity vectors are defined in (21).
Thus the entraining velocity vector v is calculated as the arithmetic average of the P and Q velocities, while the sliding velocity vector v sl is calculated as the difference between the P and Q velocities, as written in (22).
In order to write the velocity vectors in spherical coordinates, they must be rotated by the mean of the rotation matrix R s , which depends on the spherical angles θ and ϕ shown in (23).
The rotated entraining velocity vector in spherical coordinates v s is directly involved in the Reynolds equation, and it is also useful to calculate the sliding velocity vector in spherical coordinates v sl,s , as written in (24).
The norm of the sliding velocity vector in the contact plane v sl , defined in (25), will be used to evaluate the wear rate τ w and the average shear rate γ.
Then, the Reynolds equation is numerically solved to find the total pressure field evolution p t (θ, ϕ, t) for a given eccentricity n(t) and angular velocity ω(t) time functions. Once the total pressure p t (θ, ϕ, t) is calculated, it can be integrated into the domain's area to find load vector N(t), and the same is done with the linear wear field u w (θ, ϕ, t) to evaluate the volume wear V w (t), in (26).

Discrete Reynolds Equation
The spherical Reynolds problem in a MEHL regime is a non-linear Partial Differential Equation (PDE), and it can be discretized on a grid domain, represented by the discrete spherical angles θ i and ϕ j , and solved for each time step t n by using the finite difference method [5,36], building the problem in a Matlab ® environment.
A grid domain is defined and shown in Figure 6. Once the grid domain is defined, the Reynolds equation can be written through its quantities in each grid point in (27).
In order to take into account the gap geometry update due to the progressive linear wear, the film thickness is calculated in (28) considering the linear wear accumulated until the instant , so that the gap depends on the pressure only through the deformation: ℎ = ℎ + + 28 A quantity involved in the problem, for example the fluid pressure , evaluated in its discrete form on the cross stencil is written in (29)  Evaluating all of the quantities involved in the equation on the cross stencil, after defining vectors and matrices for the finite differences, like , , and , including, in (30), the products inside the derivatives in the vectors , and , the Reynolds equation is writeable, in (31), in each grid point as a scalar non-linear equation made of products between vectors and matrices.
= sin ℎ = ℎ = ℎ 30 − sin Executing the derivatives with respect to the pressure cross stencil vector of the equation , the row , of the Reynolds equation's Jacobian is obtained analytically.
Firstly the velocity derivatives with respect to pressure together with the gap derivative are given in (32).
The derivatives in (32) are useful to find the ones referred to the average shear rate and to the effective viscosity in (33). Once the grid domain is defined, the Reynolds equation can be written through its quantities in each grid point in (27) sin In order to take into account the gap geometry update due to the progressive linear wear, the film thickness is calculated in (28) considering the linear wear accumulated until the instant t n−1 , so that the gap depends on the pressure only through the deformation: A quantity involved in the problem, for example the fluid pressure p, evaluated in its discrete form on the cross stencil is written in (29) as a vector p s of five elements.
Evaluating all of the quantities involved in the equation on the cross stencil, after defining vectors and matrices for the finite differences, like D 2s , D θs , D ϕs and d t , including, in (30), the products inside the derivatives in the vectors u θs , u ϕs and u t , the Reynolds equation is writeable, in (31), in each grid point as a scalar non-linear equation R I p s made of products between vectors and matrices.
Executing the derivatives with respect to the pressure cross stencil vector p s of the equation R I , the row J R,I of the Reynolds equation's Jacobian is obtained analytically.
Firstly the velocity derivatives with respect to pressure together with the gap derivative are given in (32).
The derivatives in (32) are useful to find the ones referred to the average shear rate and to the effective viscosity in (33).
The fluid density and viscosity derivatives are given in (34).
Finally the derivatives of the term that directly compare in the Reynolds equation are calculated in (35).
Hence the Jacobian row J R,I associated to the scalar equation R I is obtained in (36).
Closing the equation set with the ones associated to the boundary conditions, the Reynolds equation set R and its analytic Jacobian J R are obtained in (37).
Then, the Newton method with relaxation is used to find the fluid pressure vector p. Starting from an initial trial pressure p (0) , in correspondence with each k iteration, the updated pressure p (k) is subjected to the cavitation constraint so that it cannot be negative. The iterative method, explained in (38), runs until the defined residual becomes lower than a chosen tolerance tol p .
The relaxation factor λ r is chosen as a function of the current eccentricity vector norm |n n |, because it has been seen that the problem requires a greater under-relaxation for a higher eccentricity vector norm for numerical convergence purposes. The function used in (39) enforces a gradual exponential under-relaxation over the unitary norm eccentricity vector through chosen parameters u r and τ r , and it is shown in Figure 7. When the Newton cycle has reached the convergence, the total pressure found is used to evaluate the linear wear, in (40), that will be elaborated in the gap calculation at the successive time step. = + Δ 40 Finally, the MEHL problem block diagram is clarified in Figure 8. When the Newton cycle has reached the convergence, the total pressure p t n ij found is used to evaluate the linear wear, in (40), that will be elaborated in the gap calculation at the successive time step.
Finally, the MEHL problem block diagram is clarified in Figure 8. The MEHL Reynolds equation solution needs eccentricity and angular velocity to find pressure, while in general the eccentricity is unknown and the aim is to find it while taking advantage of the known load.
In order to solve the inverse problem, a further Newton iterative cycle is used: starting from a trial eccentricity n (0) , the updated eccentricity n (h) is evaluated so that the calculated loads N n (h) converge to the reference ones N re f by setting to zero the function F defined in (41).
The Jacobian J F of the function F is built in the numerical way, in (42), by means of a small arbitrary increment defined by ε.
When the Newton cycle has reached the convergence, the total pressure found is used to evaluate the linear wear, in (40), that will be elaborated in the gap calculation at the successive time step. = + Δ 40 Finally, the MEHL problem block diagram is clarified in Figure 8.  The iterative cycle on the eccentricity, shown in (43), runs until a defined residual becomes lower than a chosen tolerance tol N .
The inverse MEHL problem block diagram is explained in Figure 9.
The inverse MEHL problem block diagram is explained in Figure 9.

Results and Discussion
In this paragraph the results of the simulation regarding the gait cycle are shown. The input data are represented by the Bergmann loads and the angular velocity vector is shown in Figure 3, while the parameters used for the simulation are listed in Table 1 and are comparable with the ones found in the literature and referenced in the table.

Results and Discussion
In this paragraph the results of the simulation regarding the gait cycle are shown. The input data are represented by the Bergmann loads and the angular velocity vector is shown in Figure 3, while the parameters used for the simulation are listed in Table 1 and are comparable with the ones found in the literature and referenced in the table.
The first result shown in Figure 10 is the time evolution of the dimensionless eccentricity vector n(t) components during the gait. In Figure 10 it is important to note that the y component of the eccentricity is the main component responsible for the contact between the surfaces, especially in the stance phase of the cycle, because of the combined action with the other two components, while the squeeze action is mostly due the x and z components.
The total pressure p t (θ, ϕ, t) and the gap h(θ, ϕ, t) fields are shown in two different time instants in order to analyze a full film lubrication phase of the cycle, at 1.43% of the gait, and a mixed one, at 92.9%, in which the coexistence and the transition between the lubricated areas and the contact ones are visualized.
Firstly, the surface plots for the full lubrication phase are shown in Figure 11. As expected, the peak pressure is located in the domain portion characterized by small gap supported by the deformation.

Increment for the numerical Jacobian 10
The first result shown in Figure 10 is the time evolution of the dimensionless eccentricity vector components during the gait. In Figure 10 it is important to note that the component of the eccentricity is the main component responsible for the contact between the surfaces, especially in the stance phase of the cycle, because of the combined action with the other two components, while the squeeze action is mostly due the and components.  in order to analyze a full film lubrication phase of the cycle, at 1.43% of the gait, and a mixed one, at 92.9%, in which the coexistence and the transition between the lubricated areas and the contact ones are visualized. Firstly, the surface plots for the full lubrication phase are shown in Figure 11. As expected, the peak pressure is located in the domain portion characterized by small gap supported by the deformation. The same instant is then represented in the following plots in Figure 12, where the pressure and the gap ℎ are shown in the and directions in the correspondence of the domain point , , where the pressure peak is reached. The same instant is then represented in the following plots in Figure 12, where the pressure p t and the gap h are shown in the θ and ϕ directions in the correspondence of the domain point (θ max , ϕ max ), where the pressure peak is reached.
in order to analyze a full film lubrication phase of the cycle, at 1.43% of the gait, and a mixed one, at 92.9%, in which the coexistence and the transition between the lubricated areas and the contact ones are visualized.
Firstly, the surface plots for the full lubrication phase are shown in Figure 11. As expected, the peak pressure is located in the domain portion characterized by small gap supported by the deformation. The same instant is then represented in the following plots in Figure 12, where the pressure and the gap ℎ are shown in the and directions in the correspondence of the domain point , , where the pressure peak is reached. Regarding the instant characterized by a mixed lubrication phase, the same surface plots are shown in Figure 13, in which a sudden interruption of the film thickness h(θ, ϕ) decrease is seen where it reaches the boundary layer thickness ∆ b : in this zone the pressure p t (θ, ϕ) follows the rules of the deformation model depending on the amount of contact deformation. lubrication phase along the θ direction; (b) pressure and gap in correspondence of the pressure peak point in the full film lubrication phase along the φ direction.
Regarding the instant characterized by a mixed lubrication phase, the same surface plots are shown in Figure 13, in which a sudden interruption of the film thickness ℎ , decrease is seen where it reaches the boundary layer thickness Δ : in this zone the pressure , follows the rules of the deformation model depending on the amount of contact deformation. The same instant is analyzed in Figure 14 with the plots along the and directions in correspondence with the domain point , , which is characterised by the peak pressure: in these plots the coexistence of the lubricated area and the contact one can be seen in a better way. The same instant is analyzed in Figure 14 with the plots along the θ and ϕ directions in correspondence with the domain point (θ max , ϕ max ), which is characterised by the peak pressure: in these plots the coexistence of the lubricated area and the contact one can be seen in a better way. point in the full film lubrication phase along the φ direction.
Regarding the instant characterized by a mixed lubrication phase, the same surface plots are shown in Figure 13, in which a sudden interruption of the film thickness ℎ , decrease is seen where it reaches the boundary layer thickness Δ : in this zone the pressure , follows the rules of the deformation model depending on the amount of contact deformation. The same instant is analyzed in Figure 14 with the plots along the and directions in correspondence with the domain point , , which is characterised by the peak pressure: in these plots the coexistence of the lubricated area and the contact one can be seen in a better way. It is interesting to follow over time the maximum pressure p max (t) and the minimum gap h min (t) evolution in Figure 15, in order to evaluate the magnitude order of the pressure reached in the coupling and to analyze the phases characterized by the contact: during the gait cycle, the analyzed pair is almost always in a mixed lubrication phase, while the pure full film lubrication (full film existence on the whole surface) is experienced only in the first instants of the gait characterized by a fast decrease of the gap towards the boundary layer thickness ∆ b ; hence after about the 7% of the gait cycle there was a coexistence of lubricated areas and contact ones and the load was supported by both fluid pressure and contact pressure in some locations of the domain, which can be clearly visualized in subsequent figures; the maximum pressure qualitatively follows the most critical load components' time trend; it reaches a peak value of 10.46 MPa at about 16% of the cycle and it almost reaches this value again in the second phase of the stance period at about 47% of the gait.
peak point in the mixed lubrication phase along the φ direction.
It is interesting to follow over time the maximum pressure and the minimum gap ℎ evolution in Figure 15, in order to evaluate the magnitude order of the pressure reached in the coupling and to analyze the phases characterized by the contact: during the gait cycle, the analyzed pair is almost always in a mixed lubrication phase, while the pure full film lubrication (full film existence on the whole surface) is experienced only in the first instants of the gait characterized by a fast decrease of the gap towards the boundary layer thickness Δ ; hence after about the 7% of the gait cycle there was a coexistence of lubricated areas and contact ones and the load was supported by both fluid pressure and contact pressure in some locations of the domain, which can be clearly visualized in subsequent figures; the maximum pressure qualitatively follows the most critical load components' time trend; it reaches a peak value of 10.46 MPa at about 16% of the cycle and it almost reaches this value again in the second phase of the stance period at about 47% of the gait.  Figure 16 shows the fields of fluid pressure and contact pressure on the acetabular cup surface. The instants analyzed are referred to the time interval from about 12% of the cycle to 24%, during the stance phase, when the minimum film thickness has Figure 15. Maximum pressure and minimum gap evolution during the gait cycle. Figure 16 shows the fields of fluid pressure p and contact pressure p c on the acetabular cup surface. The instants analyzed are referred to the time interval from about 12% of the cycle to 24%, during the stance phase, when the minimum film thickness has just reached the boundary thickness and the high fluid pressure is gradually replaced by the contact one during the transition (Supplementary Video S1). It is interesting to follow over time the maximum pressure and the minimum gap ℎ evolution in Figure 15, in order to evaluate the magnitude order of the pressure reached in the coupling and to analyze the phases characterized by the contact: during the gait cycle, the analyzed pair is almost always in a mixed lubrication phase, while the pure full film lubrication (full film existence on the whole surface) is experienced only in the first instants of the gait characterized by a fast decrease of the gap towards the boundary layer thickness Δ ; hence after about the 7% of the gait cycle there was a coexistence of lubricated areas and contact ones and the load was supported by both fluid pressure and contact pressure in some locations of the domain, which can be clearly visualized in subsequent figures; the maximum pressure qualitatively follows the most critical load components' time trend; it reaches a peak value of 10.46 MPa at about 16% of the cycle and it almost reaches this value again in the second phase of the stance period at about 47% of the gait.  Figure 16 shows the fields of fluid pressure and contact pressure on the acetabular cup surface. The instants analyzed are referred to the time interval from about 12% of the cycle to 24%, during the stance phase, when the minimum film thickness has The resulting total pressure p t , together with the sliding of the surfaces, causes the wear of the coupling and the gradual removal of material. In Figure 17 the linear wear field u w is shown, together with the total pressure, along the gait cycle in four instances between about 12% and 98% of the gait, in order to visualise the growth of the material removed by the surfaces (Supplementary Video S2). contact pressure transition from the full film lubrication phase to the mixed one.
The resulting total pressure , together with the sliding of the surfaces, causes the wear of the coupling and the gradual removal of material. In Figure 17 the linear wear field is shown, together with the total pressure, along the gait cycle in four instances between about 12% and 98% of the gait, in order to visualise the growth of the material removed by the surfaces (Video 2). Figure 17. (a-d) Total pressure field evolution from 12% of the cycle to 98%; (e-h) linear wear field evolution from 12% of the cycle to 98%.
The linear wear at the end of the cycle is then shown in Figure 18 in order to give an idea of the amount of worn material during a gait cycle. Together with the linear wear, which turned out to have a magnitude order comparable with the nanometers, the trajectory of the pressure peak is shown, from the red cross to the blue square, which could be evaluated as the contact point track on the acetabular cup surface. It can be seen that the contact point track initially keeps its position around the Proximo-Distal axis, which is the most loaded direction during the stance phase; it then moves to the cup center during the swing phase and finally returns to its initial position. The linear wear u w at the end of the cycle is then shown in Figure 18 in order to give an idea of the amount of worn material during a gait cycle. Together with the linear wear, which turned out to have a magnitude order comparable with the nanometers, the trajectory of the pressure peak is shown, from the red cross to the blue square, which could be evaluated as the contact point track on the acetabular cup surface. It can be seen that the contact point track initially keeps its position around the Proximo-Distal axis, which is the most loaded direction during the stance phase; it then moves to the cup center during the swing phase and finally returns to its initial position.
The resulting total pressure , together with the sliding of the surfaces, causes the wear of the coupling and the gradual removal of material. In Figure 17 the linear wear field is shown, together with the total pressure, along the gait cycle in four instances between about 12% and 98% of the gait, in order to visualise the growth of the material removed by the surfaces (Video 2). Figure 17. (a-d) Total pressure field evolution from 12% of the cycle to 98%; (e-h) linear wear field evolution from 12% of the cycle to 98%.
The linear wear at the end of the cycle is then shown in Figure 18 in order to give an idea of the amount of worn material during a gait cycle. Together with the linear wear, which turned out to have a magnitude order comparable with the nanometers, the trajectory of the pressure peak is shown, from the red cross to the blue square, which could be evaluated as the contact point track on the acetabular cup surface. It can be seen that the contact point track initially keeps its position around the Proximo-Distal axis, which is the most loaded direction during the stance phase; it then moves to the cup center during the swing phase and finally returns to its initial position. A focus on the contact point track is shown in the spherical angles θ and ϕ plane in Figure 19. A focus on the contact point track is shown in the spherical angles and plane in Figure 19. The surface integration of the linear wear , , field evolution allows evaluating the wear volume trend during the gait, as shown in Figure 20. The resulting wear volume is practically null in the initial full film lubrication phase; it then grows with a big slope during the stance phase and then the slope decreases toward the swing phase. Since the linear wear and thus the worn volume turned out to be small, it can be assumed that the geometry of the surface varies proportionally with time at least during the first million cycles, which is considered to be the amount of cycles executed by an average man in a year: thus a volumetric wear rate , can be obtained in (44) by multiplying the last value of the wear volume , at the end of the cycle by the number of cycles done in a year. Obviously, this assumption is valid for the first million cycles: the wear phenomenon over a million cycles was studied by other researchers, for example [37,38], through experimental in vitro wear tests, and it can be seen that the worn material grows proportionally with the cycles for 1-1.5 million cycles with a certain slope (running-in phase) and then the slope decreases reaching a constant value (steady-state phase). Since the obtained volumetric wear rate , is referred to as the The surface integration of the linear wear u w (θ, ϕ, t) field evolution allows evaluating the wear volume V w (t) trend during the gait, as shown in Figure 20. The resulting wear volume is practically null in the initial full film lubrication phase; it then grows with a big slope during the stance phase and then the slope decreases toward the swing phase. A focus on the contact point track is shown in the spherical angles and plane in Figure 19. The surface integration of the linear wear , , field evolution allows evaluating the wear volume trend during the gait, as shown in Figure 20. The resulting wear volume is practically null in the initial full film lubrication phase; it then grows with a big slope during the stance phase and then the slope decreases toward the swing phase. Since the linear wear and thus the worn volume turned out to be small, it can be assumed that the geometry of the surface varies proportionally with time at least during the first million cycles, which is considered to be the amount of cycles executed by an average man in a year: thus a volumetric wear rate , can be obtained in (44) by multiplying the last value of the wear volume , at the end of the cycle by the number of cycles done in a year. Obviously, this assumption is valid for the first million cycles: the wear phenomenon over a million cycles was studied by other researchers, for example [37,38], through experimental in vitro wear tests, and it can be seen that the worn material grows proportionally with the cycles for 1-1.5 million cycles with a certain slope (running-in phase) and then the slope decreases reaching a constant value (steady-state phase). Since the obtained volumetric wear rate , is referred to as the Since the linear wear and thus the worn volume turned out to be small, it can be assumed that the geometry of the surface varies proportionally with time at least during the first million cycles, which is considered to be the amount of cycles executed by an average man in a year: thus a volumetric wear rate V w,r can be obtained in (44) by multiplying the last value of the wear volume V w, f at the end of the cycle by the number of cycles n cyc done in a year.
Obviously, this assumption is valid for the first million cycles: the wear phenomenon over a million cycles was studied by other researchers, for example [37,38], through experimental in vitro wear tests, and it can be seen that the worn material grows proportionally with the cycles for 1-1.5 million cycles with a certain slope (running-in phase) and then the slope decreases reaching a constant value (steady-state phase). Since the obtained volumetric wear rate V w,r is referred to as the running-in phase, it can be intended as an indicator of the prosthesis performance in terms of duration. The volumetric wear rate value obtained is comparable with the ones found in the literature and scientific reviews. In particular, the value is close to the ones found by [10,[38][39][40] even if the latter refers to retrieved metal on metal hip replacements.
With reference to the pressure and gap profiles of Figure 12, the results showed that the classical asymmetrical shape between the inlet zone and the cavitation one was not pronounced for the used set of the prosthesis parameters and for the particular kinematical and loading data: this is probably due to the low entity of the relative sliding motion with respect to the squeezing one.
In fact, the mentioned asymmetrical shape is characteristic of purely mechanical application in which the sliding effect is amplified by the much higher angular velocity between the components.
Then, in order to see the EHL features, another simulation was conducted with simplified inputs, described in Equation (45). In particular, the following were considered: - Only the time evolution of the y component of the dimensionless eccentricity n from 0 to 1.1 over 0.1 s; - Only the time evolution of the z component of the angular velocity ω as a constant value equal to 500 rad/s, so that a faster relative sliding motion was considered.
With this set of inputs, the results can be analyzed in the rotation plane xy characterized by θ = π/2 and 0 ≤ ϕ ≤ π. Moreover, while all the parameters used are the same reported in the Table 1, a sensitivity analysis on the radial clearance c was performed in the range from 80 to 140 µm, because it represents one of the parameters that most affect the performance of a lubricated joint.
The results are shown in Figure 21 in terms of pressure, gap and linear wear profiles in the rotation plane. The volumetric wear rate value obtained is comparable with the ones found in the literature and scientific reviews. In particular, the value is close to the ones found by [10,[38][39][40] even if the latter refers to retrieved metal on metal hip replacements.
With reference to the pressure and gap profiles of Figure 12, the results showed that the classical asymmetrical shape between the inlet zone and the cavitation one was not pronounced for the used set of the prosthesis parameters and for the particular kinematical and loading data: this is probably due to the low entity of the relative sliding motion with respect to the squeezing one.
In fact, the mentioned asymmetrical shape is characteristic of purely mechanical application in which the sliding effect is amplified by the much higher angular velocity between the components.
Then, in order to see the EHL features, another simulation was conducted with simplified inputs, described in equation (45). In particular, the following were considered: -Only the time evolution of the component of the dimensionless eccentricity from 0 to 1.1 over 0.1 seconds; -Only the time evolution of the component of the angular velocity as a constant value equal to 500 rad/s, so that a faster relative sliding motion was considered.  Table 1, a sensitivity analysis on the radial clearance was performed in the range from 80 to 140 µm, because it represents one of the parameters that most affect the performance of a lubricated joint.
The results are shown in Figure 21 in terms of pressure, gap and linear wear profiles in the rotation plane. Firstly, the classical EHL shapes were visualized, in particular the growth of the pressure along the inlet zone, which takes place over a greater area than the one covered by the pressure fall-off towards the cavitation zone; then the minimum film thickness was localized at the exit of the deformation zone.
The effect of the radial clearance increase led the pressure profile to increase and to become more localized around the deformation zone, and this is due to the decrease of geometrical conformity between the surfaces' curvatures. The pressure rise kept the gap thickness almost constant, while the resulting wear became greater, as expected. Firstly, the classical EHL shapes were visualized, in particular the growth of the pressure along the inlet zone, which takes place over a greater area than the one covered by the pressure fall-off towards the cavitation zone; then the minimum film thickness was localized at the exit of the deformation zone.
The effect of the radial clearance increase led the pressure profile to increase and to become more localized around the deformation zone, and this is due to the decrease of geometrical conformity between the surfaces' curvatures. The pressure rise kept the gap thickness almost constant, while the resulting wear became greater, as expected. In Figure 22 the load x and y components are shown against the radial clearance: the y component increased with the radial clearance as expected, because of the pressure growth, while the x component decreased because the pressure bell moved towards the cup center in correspondence with the radial clearance increase, so that the area covered by the high pressure is centered around the y axis. In Figure 22 the load and components are shown against the radial clearance: the component increased with the radial clearance as expected, because of the pressure growth, while the component decreased because the pressure bell moved towards the cup center in correspondence with the radial clearance increase, so that the area covered by the high pressure is centered around the axis.

Conclusions
This paper proposed a novel numerical model with the purpose of analyzing in silico the MEHL of an artificial hip joint in the case of soft-on-hard prosthesis and to predict the implant wear during the gait.
The model was used to analyze the tribological configuration of an artificial hip joint, made of a UHMWPE acetabular cup against a rigid (metallic or ceramic) hard femoral head, subjected to the gait cycle loading and motion conditions, represented by the in vivo load measurements achieved by Bergmann [30] and the angular velocities calculated by using the musculoskeletal software OpenSim® during a normal gait cycle.
The proposed model is based on MEHL conditions during the gait, also considering the synovial fluid compressibility through the dependence of its density on the pressure, the viscosity dependence on the pressure by using the Barus model, and the shear-thinning non-Newtonian behavior in the Cross formulation. The solution was achieved by using a novel convergence algorithm based on two iterative cycles solved by the relaxed Newton iterative method, which includes the surfaces' progressive wear phenomena during the time.
The performed simulations allowed the calculation of gap thickness, fluid/contact pressure and acetabular cup stress/strain during the gait, and were able to predict the volumetric wear rate of the prosthesis. In particular: -The profiles and the shapes of the analyzed quantities are in good agreement with the scientific literature [10,18,26]; -The predicted volumetric wear is in good agreement with the values found in other researches in the same framework ( [10,[38][39][40]); -The simulation conducted in the framework of the radial clearance sensitivity analysis showed the expected tribological behavior in terms of classical EHL shapes. Even if the results obtained are satisfactory, of course the proposed model is improvable in order to characterize the prosthetic tribological performance in a more and more accurate way. Future researches will be devoted mainly to:

Conclusions
This paper proposed a novel numerical model with the purpose of analyzing in silico the MEHL of an artificial hip joint in the case of soft-on-hard prosthesis and to predict the implant wear during the gait.
The model was used to analyze the tribological configuration of an artificial hip joint, made of a UHMWPE acetabular cup against a rigid (metallic or ceramic) hard femoral head, subjected to the gait cycle loading and motion conditions, represented by the in vivo load measurements achieved by Bergmann [30] and the angular velocities calculated by using the musculoskeletal software OpenSim ® during a normal gait cycle.
The proposed model is based on MEHL conditions during the gait, also considering the synovial fluid compressibility through the dependence of its density on the pressure, the viscosity dependence on the pressure by using the Barus model, and the shear-thinning non-Newtonian behavior in the Cross formulation. The solution was achieved by using a novel convergence algorithm based on two iterative cycles solved by the relaxed Newton iterative method, which includes the surfaces' progressive wear phenomena during the time.
The performed simulations allowed the calculation of gap thickness, fluid/contact pressure and acetabular cup stress/strain during the gait, and were able to predict the volumetric wear rate of the prosthesis. In particular: - The profiles and the shapes of the analyzed quantities are in good agreement with the scientific literature [10,18,26]; - The predicted volumetric wear is in good agreement with the values found in other researches in the same framework ( [10,[38][39][40]); - The simulation conducted in the framework of the radial clearance sensitivity analysis showed the expected tribological behavior in terms of classical EHL shapes.
Even if the results obtained are satisfactory, of course the proposed model is improvable in order to characterize the prosthetic tribological performance in a more and more accurate way. Future researches will be devoted mainly to: -Build a finite element model to calculate the deformation field of both acetabular cup and femoral head, in order to analyze the tribological behavior of other types of THRs in which both contact bodies must be considered deformable; -Consider the surfaces' topography in the gap calculation in order to analyze a more realistic surface in the mixed lubrication mode and also to consider material transfer phenomena; -Improve the wear modeling in order to consider more specific wear modes, such as, for example, delamination effects.