Next Article in Journal
Numerical and Experimental Analysis of Nonlinear Vibrational Response due to Pressure-Dependent Interface Stiffness
Next Article in Special Issue
On the Design and Lubrication of Water-Lubricated, Rubber, Cutlass Bearings Operating in the Soft EHL Regime
Previous Article in Journal / Special Issue
Mapping the Micro-Abrasion Mechanisms of CoCrMo: Some Thoughts on Varying Ceramic Counterface Diameter on Transition Boundaries In Vitro
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

by
Alessandro Ruggiero
* and
Alessandro Sicilia
Department of Industrial Engineering, University of Salerno, Via Giovanni Paolo II, nº 132, 84084 Fisciano, Italy
*
Author to whom correspondence should be addressed.
Lubricants 2020, 8(7), 72; https://doi.org/10.3390/lubricants8070072
Submission received: 5 May 2020 / Revised: 24 June 2020 / Accepted: 28 June 2020 / Published: 1 July 2020

Abstract

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

Graphical Abstract

1. 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 analyzed joint coming from the human daily activities (gait, running, climbing and so on) with the purpose to optimize the bio-bearing tribological design especially for minimizing its wear rate.
The lubrication of a natural human synovial joint is a complex mix of elasto-hydrodynamic, full film and boundary lubrication modes and it is still under investigation.
Regarding the artificial joint, the lubrication phenomenon can be analyzed by considering a combination of Boundary Lubrication (BL), Mixed Lubrication (ML), Hydrodynamic Lubrication (HL) and Elasto-Hydrodynamic Lubrication (EHL) depending on load and motion conditions [8,9].
The EHL is a particular case of the hydrodynamic mode characterized by the deformation of the joint surfaces due to the high pressure reached in the gap, and it is typical of tribo-systems with a low geometrical conformity and with low elasticity. The surfaces’ deformation due to the load allows the survival of very low constant fluid film thickness, which contributes to minimizing the surfaces’ wear.
The synovial fluid pressure profile in this zone substantially follows the Hertzian pressure, and the minimum film thickness is located at the exit of this area, causing the well-known pressure spike. The combination of high load and low relative motion could lead to contact between the articulated surfaces despite their deformation, so that in these contact areas the contact pressure rises and in this case the regime is referred to as Mixed Elasto-Hydrodynamic Lubrication (MEHL).
Among the synovial joints, which are the human joints allowing the biggest range of relative motion between the linked bones, one of the most important is the hip, which can be modeled as a spherical joint linking the femur and the pelvis. Due to the hip anatomical position, the hip spherical surfaces, namely the femoral head and the acetabular cup, have to bear a load greater than the body weight during most daily activities.
Moreover, the wear phenomenon analysis in the hip joint is affected by several factors, and it is difficult to identify the wear type that occurs within it [10], especially considering the different biomaterials used in the Total Hip Replacement (THR), mainly polymers, metals and ceramics [11,12,13,14].
Since the MEHL lubrication problem is mostly approached by solving the Reynolds equation, the introduction of surface deformation together with the dependence of some parameters on the pressure and the contact areas brings with it a series of difficulties related to the high non-linearity of the problem and to the transition between the full lubricated areas and the contact ones.
An important role is played by the model that evaluates the surfaces’ deformation starting from the pressure field: some authors, as explained successively, used fast Fourier transform techniques based on the Finite Element Method, FEM, models regarding mainly hard-on-hard pairs; other authors justified the usage of easier models, like the constrained column one, when the prosthesis is a soft-on-hard coupling [15,16,17].
In general, some authors developed models that are able to analyze the tribological performance of artificial hip joints and interesting research studies have been found in the scientific literature.
In [18], the authors studied the effects of the fluid viscosity variation on the contact area regarding a hip implant made of a UHMWPE acetabular cup and a metallic femoral head through a full numerical methodology developed for the mixed lubrication mode, while in [6] a general tribological analysis was conducted for spherical bearings with complex spherical-based geometry by using models to simulate interesting tribological quantities like wear, regarding a mixed EHL regime.
In [7], the authors analyzed the variation of the contact pressure, the contact area and the wear volume loss over several cycles regarding a metal-on-metal hip joint replacement through a finite-element model coupled with Archard’s equation, validating the model with hip simulators’ experimental data.
The analysis of the transition between the full film lubricated areas and the contact ones is a critical issue in mixed lubrication modeling. The simplest idea is represented by the introduction of a contact model activating in the domain portion where the fluid film thickness is lower than a certain limit chosen according to the boundary layer that appears in a tribological system when the two surfaces are excessively pushed against each other, leading to the rising of the contact pressure, which, together with the fluid pressure of the lubricated areas, contributes to the load support.
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.

2. Materials and Methods

2.1. 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 x y z 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.
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.
{ N h i p = [ N A P N P D N M L ] T ω h i p = [ ω A A ω I E R ω F E ] T
The transformation is obtained in (2) by the rotation matrix R c , which takes into account the inclination angle α i n , defined as the angle formed by the cup axis and the sagittal plane, and the anteversion angle β a v , defined as the angle between the cup axis projection on the sagittal plane and the AP axis.
{ R x ( α ) = [ 1 0 0 0 cos ( α ) sin ( α ) 0 sin ( α ) cos ( α ) ] R z ( γ ) = [ cos ( γ ) sin ( γ ) 0 sin ( γ ) cos ( γ ) 0 0 0 1 ] R c = R z ( π 2 β a v ) R x ( α i n )
Finally, the vectors N and ω defined in the cup reference frame are obtained in (3).
{ N = R c T N h i p ω = R c T ω h i p
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.

2.2. 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 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 .
{ · ( f p ) = · ( ρ h v ) + ( ρ h ) t p Ω = p 0 , f = ρ h 3 12 μ
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.
ρ ( p ) = ρ 0 a ρ + b ρ p a ρ + p
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.
μ ( γ ¯ , p ) = μ e f f ( γ ¯ ) e α μ p
The effective viscosity, in (7), depends on the average shear rate γ ¯ , that is the ratio between the sliding velocity v s l 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 .
μ e f f ( γ ¯ ) = μ + μ 0 μ 1 + k γ ¯ n γ ¯ = v s l h
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.
h = h g + δ + u w
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.
{ δ f = D ( p ) δ c = max { 0 Δ b ( h g + δ f + u w ) }   δ = δ f + δ c
The calculation of the contact deformation leads to the determination of the contact pressure p c , in (10), with the same deformation model inversion.
p c = D 1 ( δ c )
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].
D ( p ) = k d p
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.
p t = p + p c
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 w 0 , 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 ).
k w = k w 0 ( h R a ) α w τ w = k w p t v s l
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.
u w = 0 t τ w d t
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).
{ N ( t ) = p t d A V w ( t ) = u w d A

2.3. The Spherical Joint

Regarding a spherical joint, that is the hip case, the Reynolds equation is written in the cup reference frame x y z , 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 r ^ as a function of the chosen spherical angles θ and φ written in (16).
r ^ ( θ , φ ) = [ sin ( θ ) cos ( φ ) sin ( θ ) sin ( φ ) cos ( θ ) ] T 0 θ , φ π
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).
{ sin θ θ ( sin θ f p θ ) + φ ( f p φ ) = R sin θ [ θ ( sin θ ρ h U θ ) + φ ( ρ h U φ ) + R sin θ ( ρ h ) t ] p ( 0 , φ , t ) = p ( π , φ , t ) = p 0 p ( θ , 0 , t ) = p ( θ , π , t ) = p 0
The geometrical gap h g is calculated in (18) with the dimensionless eccentricity vector n and the radial clearance c .
{ c = R r n = e c h g ( θ , φ , t ) = c ( 1 n ( t ) T r ^ ( θ , φ ) )
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).
k d = R [ ( 1 + H R ) 3 1 ] E [ 1 1 2 ν + 2 1 + ν ( 1 + H R ) 3 ]
Thus the fluid film thickness in the Reynolds equation is obtained in (20).
h ( θ , φ , t ) = h g ( θ , φ , t ) + δ ( θ , φ , t ) + u w ( θ , φ , t )
Since the cup is fixed with the reference frame, the velocity v P of a point P located by r ^ 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).
{ v P ( θ , φ , t ) = 0 v Q ( θ , φ , t ) = e ˙ ( t ) + ω ( t ) × [ ( R h ) r ^ ( θ , φ ) e ( t ) ]
Thus the entraining velocity vector v is calculated as the arithmetic average of the P and Q velocities, while the sliding velocity vector v s l is calculated as the difference between the P and Q velocities, as written in (22).
{ v = v P + v Q 2 = v Q 2 v s l = v P v Q = v Q
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).
{ R y ( β ) = [ cos ( β ) 0 sin ( β ) 0 1 0 sin ( β ) 0 cos ( β ) ] R z ( γ ) = [ cos ( γ ) sin ( γ ) 0 sin ( γ ) cos ( γ ) 0 0 0 1 ] R s = R z ( φ ) R y ( θ ) #
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 s l , s , as written in (24).
{ v s = R s T v = 1 2 R s T { e ˙ + ω × [ ( R h ) r ^ e ] } = [ U θ U φ U ρ ] T v s l , s = 2 v s
The norm of the sliding velocity vector in the contact plane v s l , defined in (25), will be used to evaluate the wear rate τ w and the average shear rate γ ¯ .
v s l = v s l , s , θ 2 + v s l , s , φ 2 = 2 U θ 2 + U φ 2
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).
{ N ( t ) = 0   0 π   π p t ( θ , φ , t ) r ^ ( θ , φ ) R 2 sin ( θ ) d θ d φ V w ( t ) = 0   0 π   π u w ( θ , φ , t ) R 2 sin ( θ ) d θ d φ

2.4. 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)
[ sin θ θ ( f sin θ p θ ) + φ ( f p φ ) ] i j n = { R sin θ [ θ ( sin θ ρ h U θ ) + φ ( ρ h U φ ) + R sin θ t ( ρ h ) ] } i j n
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:
h i j n = h g i j n + δ i j n + u w i j n 1
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.
p s = [ p i 1 , j n p i , j 1 n p i j n p i , j + 1 n p i + 1 , j n ] T
Evaluating all of the quantities involved in the equation on the cross stencil, after defining vectors and matrices for the finite differences, like D 2 s , 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.
u θ = sin θ ρ h U θ u φ = ρ h U φ u t = ρ h
f s T D 2 s p s R sin θ i ( D θ s T u θ s + D φ s T u φ s + R sin θ i d t T u t ) = R I ( p s ) = 0
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).
{ h p = δ p = k d v s p = 1 2 R s T [ ω × ( h p r ^ ) ] = [ U θ p U φ p U ρ p ] T v s l p = 4 v s l ( U θ U θ p + U φ U φ p )
The derivatives in (32) are useful to find the ones referred to the average shear rate and to the effective viscosity in (33).
{ γ ¯ p = 1 h 2 ( v s l p h v s l h p ) μ e f f γ ¯ = k n γ ¯ n 1 μ 0 μ ( 1 + k γ ¯ n ) 2
The fluid density and viscosity derivatives are given in (34).
{ ρ p = ρ 0 a ρ b ρ 1 ( a ρ + p ) 2 μ p = μ e f f γ ¯ γ ¯ p e α μ p + α μ μ
Finally the derivatives of the term that directly compare in the Reynolds equation are calculated in (35).
{ f p = ρ p h 3 12 μ + h p ρ h 2 4 μ ρ h 3 12 μ 2 μ p u θ p = sin θ ( ρ p h U θ + ρ h p U θ + ρ h U θ p ) u φ p = ρ p h U φ + ρ h p U φ + ρ h U φ p u t p = ρ p h + ρ h p
Hence the Jacobian row J R , I associated to the scalar equation R I is obtained in (36).
p s T D 2 s f s p s + f s T D 2 s R sin θ i ( D θ s T u θ s p s + D φ s T u φ s p s + R sin θ i d t T u t p s ) = J R , I ( p s ) = R I p s
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).
R ( p ) = 0 J R ( p ) = R p
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 t o l p .
p ( 0 ) w h i l e   r e s ( k ) > t o l p p ( k + 1 ) = max [ 0 p ( k ) λ r J R 1 ( p ( k ) ) R ( p ( k ) ) ] r e s ( k + 1 ) = I | p I ( k + 1 ) p I ( k ) | I p I ( k )
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.
λ r ( | n n | ) = 1 { 0 i f   | n n | 1 ( 1 u r ) ( 1 e | n n | 1 τ r ) i f   | n n | > 1
When the Newton cycle has reached the convergence, the total pressure p t i j n found is used to evaluate the linear wear, in (40), that will be elaborated in the gap calculation at the successive time step.
u w i j n = u w i j n 1 + τ w i j n Δ t
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 r e f by setting to zero the function F defined in (41).
F ( n ) = N ( n ) N r e f = 0
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 ε .
J F ( n ) = F n [ J F ] i j = F i n j F i ( n j ( 1 + ε ) ) F i ( n j ) ε n j
The iterative cycle on the eccentricity, shown in (43), runs until a defined residual becomes lower than a chosen tolerance t o l N .
n ( 0 ) w h i l e   r e s ( h ) > t o l N n ( h + 1 ) = n ( h ) J F 1 ( n ( h ) ) F ( n ( h ) ) r e s ( h + 1 ) = | F ( n ( h ) ) | | N r e f |
The inverse MEHL problem block diagram is explained in Figure 9.

3. 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.
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 ( θ m a x , φ m a x ) , 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.
The same instant is analyzed in Figure 14 with the plots along the θ and φ directions in correspondence with the domain point ( θ m a x , φ m a x ) , 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 m a x ( t ) and the minimum gap h m i n ( 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.
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).
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).
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.
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 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.
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 c y c 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.
V w , r = V w , f n c y c = 27.35   m m 3 y e a r
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.
n = [ 0 1.1 t 0.1 0 ] ω = [ 0 0 500 ] r a d s
With this set of inputs, the results can be analyzed in the rotation plane x y 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.
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.

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

Supplementary Materials

The following are available online at https://www.mdpi.com/2075-4442/8/7/72/s1.

Author Contributions

Conceptualization, A.R.; methodology, A.R., A.S.; software, A.S.; formal analysis, A.R., A.S.; writing A.S.; supervision, A.R.; project administration, A.R.; funding acquisition, A.R. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by MIUR, PRIN 2017- BIONIC.

Conflicts of Interest

The authors declare no conflict of interest.

Nomenclature

FE, AA, IERFlexion–Extension, Adduction–Abduction, Internal–External Rotation
ML, AP, PDMedio–Lateral, Anterior–Posterior, Proximo–Distal hip axes
x , y , z Cartesian axes
N Load vector
ω Angular velocity vector
α i n , β a v Acetabular cup inclination angle and anteversion angle
R x i ( θ ) Rotation matrix for a rotation θ around the x i axis
R c Rotation matrix from the hip reference frame to the cup one
p Fluid pressure
h Fluid film thickness
ρ Fluid density
ρ 0 Nominal fluid density
a ρ , b ρ Dowson–Higginson coefficients
μ Fluid viscosity
μ e f f Effective nominal viscosity
α μ Barus exponential coefficient
μ 0 , μ Cross viscosities in correspondence of 0 and theoretically shear rate
k , n Cross viscosity model parameters
γ ¯ Average shear rate
h g Geometrical gap
δ Total surfaces’ deformation
δ f , δ c Deformation due to fluid pressure and contact deformation
u w Linear wear
Δ b Boundary layer thickness
p c Contact pressure
p t Total pressure
τ w Wear rate
k w , k w 0 Wear factor function and nominal wear factor
α w Modified Archard model exponential coefficient
R a Average roughness
V w Wear volume
θ , φ Spherical angles
t Time
r ^ Radial unit vector
r , R Femoral head and acetabular cup radii
c Radial clearance
H Acetabular cup thickness
E , ν Acetabular cup Young modulus and Poisson coefficient
k d Constraint column model constant
p 0 Boundary pressure
e , n Eccentricity vector and its dimensionless form
v , v s l Entraining and sliding velocity vectors
R s Spherical rotation matrix
i , j , n Finite difference subscripts
p Discretized pressure vector
R , J R Discretized Reynolds equation vector and its analytical Jacobian matrix
λ r , u r , τ r Relaxation factor and its parameters
r e s , t o l Iterative cycles residual and tolerance
F , J F Load difference function vector and its numerical Jacobian matrix
V w , r Volumetric wear rate

References

  1. Pramanik, S.; Agarwal, A.K.; Rai, K. Chronology of total hip joint replacement and materials development. Trends Biomater. Artif. Organs 2005, 19, 15–26. [Google Scholar]
  2. Carr, A.J.; Robertsson, O.; Graves, S.; Price, A.J.; Arden, N.K.; Judge, A.; Beard, D.J. Knee replacement. Lancet 2012, 379, 1331–1340. [Google Scholar] [CrossRef]
  3. Jin, Z.; Yonggang, M.; Yuanzhong, H.; Jianbin, L. In Memoriam: Duncan Dowson (1928–2020). Friction 2020, 8, 1–3. [Google Scholar] [CrossRef] [Green Version]
  4. Fisher, J.; Dowson, D. Tribology of total artificial joints. Proc. Inst. Mech. Eng. Part H J. Eng. Med. 1991, 205, 73–79. [Google Scholar] [CrossRef]
  5. Ruggiero, A.; Sicilia, A. Lubrication modeling and wear calculation in artificial hip joint during the gait. Tribol. Int. 2020, 142, 105993. [Google Scholar] [CrossRef]
  6. Wang, F.; Wang, L.; Sun, M. Tribological modelling of spherical bearings with complex spherical-based geometry and motion. WIT Trans. Eng. Sci. 2010, 66, 3–15. [Google Scholar]
  7. Harun, M.; Wang, F.; Jin, Z.; Fisher, J. Long-term contact-coupled wear prediction for metal-on-metal total hip joint replacement. Proc. Inst. Mech. Eng. Part J J. Eng. Tribol. 2009, 223, 993–1001. [Google Scholar] [CrossRef]
  8. Dumbleton, J.H. Tribology of Natural and Artificial Joints; Elsevier: Amsterdam, The Netherlands, 1981. [Google Scholar]
  9. Nordin, M.; Frankel, V.H. Basic Biomechanics of the Musculoskeletal System; Lippincott Williams & Wilkins: Philadelphia, PA, USA, 2001. [Google Scholar]
  10. Mattei, L.; di Puccio, F.; Piccigallo, B.; Ciulli, E. Lubrication and wear modelling of artificial hip joints: A review. Tribol. Int. 2011, 44, 532–549. [Google Scholar] [CrossRef]
  11. Zivic, F.; Affatato, S.; Trajanovic, M.; Schnabelrauch, M.; Grujovic, N.; Choy, K.L. Biomaterials in Clinical Practice: Advances in Clinical Research and Medical Devices; Springer: Berlin, Germany, 2017. [Google Scholar]
  12. Affatato, S.; Ruggiero, A.; Merola, M. Advanced biomaterials in hip joint arthroplasty. A review on polymer and ceramics composites as alternative bearings. Compos. Part B Eng. 2015, 83, 276–283. [Google Scholar] [CrossRef]
  13. Merola, M.; Affatato, S. Materials for hip prostheses: A review of wear and loading considerations. Materials 2019, 12, 495. [Google Scholar] [CrossRef] [Green Version]
  14. Aherwar, A.; Singh, A.K.; Patnaik, A. Current and future biocompatibility aspects of biomaterials for hip prosthesis. AIMS Bioeng. 2015, 3, 23–43. [Google Scholar] [CrossRef]
  15. Jalali-Vahid, D.; Jagatia, M.; Jin, Z.; Dowson, D. Prediction of lubricating film thickness in a ball-in-socket model with a soft lining representing human natural and artificial hip joints. Proc. Inst. Mech. Eng. Part J J. Eng. Tribol. 2001, 215, 363–372. [Google Scholar] [CrossRef]
  16. Wang, F.; Jin, Z. Prediction of elastic deformation of acetabular cups and femoral heads for lubrication analysis of artificial hip joints. Proc. Inst. Mech. Eng. Part J J. Eng. Tribol. 2004, 218, 201–209. [Google Scholar] [CrossRef]
  17. Jalali-Vahid, D.; Jin, Z.; Dowson, D. Elastohydrodynamic lubrication analysis of hip implants with ultra high molecular weight polyethylene cups under transient conditions. Proc. Inst. Mech. Eng. Part C J. Mech. Eng. Sci. 2003, 217, 767–777. [Google Scholar] [CrossRef]
  18. Wang, F.; Jin, Z. Lubrication modelling of artificial hip joints: From fluid film to boundary lubrication regimes. In Proceedings of the Engineering Systems Design and Analysis, Manchester, UK, 19–22 July 2004; pp. 605–611. [Google Scholar]
  19. Wang, F.; Brockett, C.; Williams, S.; Udofia, I.; Fisher, J.; Jin, Z. Lubrication and friction prediction in metal-on-metal hip implants. Phys. Med. Biol. 2008, 53, 1277. [Google Scholar] [CrossRef]
  20. Zhu, D.; Hu, Y.-Z. The study of transition from elastohydrodynamic to mixed and boundary lubrication. STLE/ASME HS Cheng Tribol. Surveill 1999, 150–156. [Google Scholar]
  21. Zhu, D. On some aspects of numerical solutions of thin-film and mixed elastohydrodynamic lubrication. Proc. Inst. Mech. Eng. Part J J. Eng. Tribol. 2007, 221, 561–579. [Google Scholar] [CrossRef]
  22. Gao, L.; Wang, F.; Yang, P.; Jin, Z. Effect of 3D physiological loading and motion on elastohydrodynamic lubrication of metal-on-metal total hip replacements. Med. Eng. Phys. 2009, 31, 720–729. [Google Scholar] [CrossRef]
  23. Jin, Z.; Dowson, D. A full numerical analysis of hydrodynamic lubrication in artificial hip joint replacements constructed from hard materials. Proc. Inst. Mech. Eng. Part C J. Mech. Eng. Sci. 1999, 213, 355–370. [Google Scholar] [CrossRef]
  24. Askari, E.; Andersen, M. A modification on velocity terms of Reynolds equation in a spherical coordinate system. Tribol. Int. 2019, 131, 15–23. [Google Scholar] [CrossRef] [Green Version]
  25. Jalali-Vahid, D.; Jagatia, M.; Jin, Z.; Dowson, D. Prediction of lubricating film thickness in UHMWPE hip joint replacements. J. Biomech. 2001, 34, 261–266. [Google Scholar] [CrossRef]
  26. Jalali-Vahid, D.; Jin, Z. Transient elastohydrodynamic lubrication analysis of ultra-high molecular weight polyethylene hip joint replacements. Proc. Inst. Mech. Eng. Part C J. Mech. Eng. Sci. 2001, 216, 409–420. [Google Scholar] [CrossRef]
  27. Udofia, I.; Jin, Z. Elastohydrodynamic lubrication analysis of metal-on-metal hip-resurfacing prostheses. J. Biomech. 2003, 36, 537–544. [Google Scholar] [CrossRef]
  28. Wang, F.; Jin, Z. Transient elastohydrodynamic lubrication of hip joint implants. J. Tribol. 2008, 130, 011007. [Google Scholar] [CrossRef]
  29. Kharb, A.; Saini, V.; Jain, Y.; Dhiman, S. A review of gait cycle and its parameters. IJCEM Int. J. Comput. Eng. Manag. 2011, 13, 78–83. [Google Scholar]
  30. Bergmann, G.; Bender, A.; Dymke, J.; Duda, G.; Damm, P. Standardized loads acting in hip implants. PLoS ONE 2016, 11, e0155612. [Google Scholar] [CrossRef]
  31. van Leeuwen, H. The determination of the pressure–viscosity coefficient of a lubricant through an accurate film thickness formula and accurate film thickness measurements. Proc. Inst. Mech. Eng. Part J J. Eng. Tribol. 2009, 223, 1143–1163. [Google Scholar] [CrossRef] [Green Version]
  32. Gao, L.; Dowson, D.; Hewson, R.W. A numerical study of non-Newtonian transient elastohydrodynamic lubrication of metal-on-metal hip prostheses. Tribol. Int. 2016, 93, 486–494. [Google Scholar] [CrossRef] [Green Version]
  33. Archard, J.; Hirst, W. The wear of metals under unlubricated conditions. Proc. R. Soc. Lond. Ser. A Math. Phys. Sci. 1956, 236, 397–410. [Google Scholar]
  34. Gao, L.; Dowson, D.; Hewson, R.W. Predictive wear modeling of the articulating metal-on-metal hip replacements. J. Biomed. Mater. Res. Part B Appl. Biomater. 2017, 105, 497–506. [Google Scholar] [CrossRef]
  35. Gao, L.; Hua, Z.; Hewson, R.; Andersen, M.S.; Jin, Z. Elastohydrodynamic lubrication and wear modelling of the knee joint replacements with surface topography. Biosurface Biotribol. 2018, 4, 18–23. [Google Scholar] [CrossRef]
  36. Jagatia, M.; Jin, Z. Elastohydrodynamic lubrication analysis of metal-on-metal hip prostheses under steady state entraining motion. Proc. Inst. Mech. Eng. Part H J. Eng. Med. 2001, 215, 531–541. [Google Scholar] [CrossRef] [PubMed]
  37. Yan, Y.; Neville, A.; Dowson, D. Biotribocorrosion—An appraisal of the time dependence of wear and corrosion interactions: I. The role of corrosion. J. Phys. D Appl. Phys. 2006, 39, 3200. [Google Scholar] [CrossRef]
  38. di Puccio, F.; Mattei, L. Biotribology of artificial hip joints. World J. Orthop. 2015, 6, 77. [Google Scholar] [CrossRef] [PubMed]
  39. Mattei, L.; di Puccio, F.; Ciulli, E. A comparative study of wear laws for soft-on-hard hip implants using a mathematical wear model. Tribol. Int. 2013, 63, 66–77. [Google Scholar] [CrossRef]
  40. Bergiers, S.; Hothi, H.; Richards, R.; Henckel, J.; Hart, A. Quantifying the bearing surface wear of retrieved hip replacements. Biosurface Biotribol. 2019, 5, 28–33. [Google Scholar] [CrossRef]
Figure 1. Hip joint reference frames. FE: Flexion–Extension rotation; ML: Medio–Lateral axis; AA: Adduction–Abduction rotation; AP: Anterior–Posterior axis; IER: Internal–External Rotation; PD: Proximo–Distal axis.
Figure 1. Hip joint reference frames. FE: Flexion–Extension rotation; ML: Medio–Lateral axis; AA: Adduction–Abduction rotation; AP: Anterior–Posterior axis; IER: Internal–External Rotation; PD: Proximo–Distal axis.
Lubricants 08 00072 g001
Figure 2. (a) Hip loads [30]; (b) hip angular velocities.
Figure 2. (a) Hip loads [30]; (b) hip angular velocities.
Lubricants 08 00072 g002
Figure 3. (a) Hip loads in the cup reference frame; (b) hip angular velocities in the cup reference frame.
Figure 3. (a) Hip loads in the cup reference frame; (b) hip angular velocities in the cup reference frame.
Lubricants 08 00072 g003
Figure 4. Fluid column within the moving surfaces.
Figure 4. Fluid column within the moving surfaces.
Lubricants 08 00072 g004
Figure 5. Acetabular cup reference frame.
Figure 5. Acetabular cup reference frame.
Lubricants 08 00072 g005
Figure 6. Spherical grid domain.
Figure 6. Spherical grid domain.
Lubricants 08 00072 g006
Figure 7. Relaxation factor function.
Figure 7. Relaxation factor function.
Lubricants 08 00072 g007
Figure 8. Mixed Elasto-Hydrodynamic Lubrication (MEHL) problem block diagram.
Figure 8. Mixed Elasto-Hydrodynamic Lubrication (MEHL) problem block diagram.
Lubricants 08 00072 g008
Figure 9. Inverse MEHL problem block diagram.
Figure 9. Inverse MEHL problem block diagram.
Lubricants 08 00072 g009
Figure 10. Eccentricity vector during the gait cycle.
Figure 10. Eccentricity vector during the gait cycle.
Lubricants 08 00072 g010
Figure 11. (a) Pressure field in the full film lubrication phase; (b) gap field in the full film lubrication phase.
Figure 11. (a) Pressure field in the full film lubrication phase; (b) gap field in the full film lubrication phase.
Lubricants 08 00072 g011
Figure 12. (a) Pressure and gap in correspondence of the pressure peak point in the full film 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.
Figure 12. (a) Pressure and gap in correspondence of the pressure peak point in the full film 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.
Lubricants 08 00072 g012
Figure 13. (a) Pressure field in the mixed lubrication phase; (b) gap field in the mixed lubrication phase.
Figure 13. (a) Pressure field in the mixed lubrication phase; (b) gap field in the mixed lubrication phase.
Lubricants 08 00072 g013
Figure 14. (a) Pressure and gap in correspondence with the pressure peak point in the mixed lubrication phase along the θ direction; (b) pressure and gap in correspondence with the pressure peak point in the mixed lubrication phase along the φ direction.
Figure 14. (a) Pressure and gap in correspondence with the pressure peak point in the mixed lubrication phase along the θ direction; (b) pressure and gap in correspondence with the pressure peak point in the mixed lubrication phase along the φ direction.
Lubricants 08 00072 g014
Figure 15. Maximum pressure and minimum gap evolution during the gait cycle.
Figure 15. Maximum pressure and minimum gap evolution during the gait cycle.
Lubricants 08 00072 g015
Figure 16. (a–d) Fluid pressure transition from the full film lubrication phase to the mixed one; (e–h) contact pressure transition from the full film lubrication phase to the mixed one.
Figure 16. (a–d) Fluid pressure transition from the full film lubrication phase to the mixed one; (e–h) contact pressure transition from the full film lubrication phase to the mixed one.
Lubricants 08 00072 g016
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%.
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%.
Lubricants 08 00072 g017
Figure 18. Worn cup at the end of the gait cycle with contact point track.
Figure 18. Worn cup at the end of the gait cycle with contact point track.
Lubricants 08 00072 g018
Figure 19. Contact point track in the spherical angles plane.
Figure 19. Contact point track in the spherical angles plane.
Lubricants 08 00072 g019
Figure 20. Cumulated wear volume during the gait.
Figure 20. Cumulated wear volume during the gait.
Lubricants 08 00072 g020
Figure 21. (a) Pressure, (b) gap and (c) linear wear profiles in the rotation plane in correspondence with different radial clearance values.
Figure 21. (a) Pressure, (b) gap and (c) linear wear profiles in the rotation plane in correspondence with different radial clearance values.
Lubricants 08 00072 g021
Figure 22. Load components in the rotation plane against different radial clearance values.
Figure 22. Load components in the rotation plane against different radial clearance values.
Lubricants 08 00072 g022
Table 1. Simulation input parameters.
Table 1. Simulation input parameters.
ParameterValue
Acetabular cup inclination angle α i n 45 ° [10]
Acetabular cup anteversion angle β a v 90 °
Femoral head radius r 14   m m [10]
Radial clearance c 100   μ m [10]
Acetabular cup thickness H 9.5   m m [15]
Acetabular cup Young modulus E 1   G P a [10]
Acetabular cup Poisson ratio ν 0.4 [10]
Barus model exponential factor α μ 19.8 · 10 9   P a 1 [31]
Cross model upper limit viscosity μ 0 40   P a   s [32]
Cross model lower limit viscosity μ 9   m P a   s [32]
Cross model parameter k 9.54 [32]
Cross model parameter n 0.73 [32]
Synovial fluid nominal density ρ 0 850   k g / m 3
Density model parameter a ρ 5.9 · 10 8   P a [31]
Density model parameter b ρ 1.34 [31]
Boundary pressure p 0 0   P a
Boundary layer thickness Δ b 30   n m [19,35]
Roughness R a 0.1   μ m [10]
Nominal wear factor k w 0 10 7   m m 3 / N / m [10]
Archard model exponential factor α w 2.24 [34]
Spherical domain mesh density 100 × 100
Time domain steps number 70
Pressure tolerance t o l p 10 5
Load tolerance t o l N 10 2
Relaxation parameter u r 0.2
Relaxation parameter τ r 0.2
Increment for the numerical Jacobian ε 10 4

Share and Cite

MDPI and ACS Style

Ruggiero, A.; Sicilia, A. A Mixed Elasto-Hydrodynamic Lubrication Model for Wear Calculation in Artificial Hip Joints. Lubricants 2020, 8, 72. https://doi.org/10.3390/lubricants8070072

AMA Style

Ruggiero A, Sicilia A. A Mixed Elasto-Hydrodynamic Lubrication Model for Wear Calculation in Artificial Hip Joints. Lubricants. 2020; 8(7):72. https://doi.org/10.3390/lubricants8070072

Chicago/Turabian Style

Ruggiero, Alessandro, and Alessandro Sicilia. 2020. "A Mixed Elasto-Hydrodynamic Lubrication Model for Wear Calculation in Artificial Hip Joints" Lubricants 8, no. 7: 72. https://doi.org/10.3390/lubricants8070072

APA Style

Ruggiero, A., & Sicilia, A. (2020). A Mixed Elasto-Hydrodynamic Lubrication Model for Wear Calculation in Artificial Hip Joints. Lubricants, 8(7), 72. https://doi.org/10.3390/lubricants8070072

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

Article Metrics

Back to TopTop