Lubrication modeling and wear calculation in artificial hip joint during the gait

Tribology International 142 (2020) 105993

Contents lists available at ScienceDirect

Tribology International journal homepage: http://www.elsevier.com/locate/triboint

Lubrication modeling and wear calculation in artificial hip joint during the gait Alessandro Ruggiero *, Alessandro Sicilia Department of Industrial Engineering, University of Salerno, Italy

A R T I C L E I N F O

A B S T R A C T

Keywords: Total hip replacement Lubrication Gait Wear

In the framework of Total Hip Replacement (THR), the aim of this paper was to propose a novel 3D contactlubrication model for the in-silico prosthetic wear calculation during the gait. The approach is based on the Reynolds equation and in-vivo loads were selected for the simulations. The hip relative angular velocities were calculated by a commercial musculoskeletal software. In the considered conditions, by using the Archard wear theory, the model allows the determination of the unsteady fluid film pressure and the volumetric wear rate prediction. With reference to a Ceramic on Polyethylene hip prosthesis, the volumetric wear rate was calculated showing a good agreement with the values obtained by other Authors as a result of FEM approaches or in-vitro tests.

1. Introduction Arthroplasty is a surgical procedure consisting in the physical replacement of an unhealthy natural human synovial joint with an artificial one. The artificial joints have to guarantee biocompatibility, fixation, mobility, load capacity, stability and minimal friction and wear of the tribo-system [1]. Nowadays the development of more and more accurate computa­ tional methods for the pre-clinical wear prediction of artificial implants are very actual issue [2], such as for example in Total Knee Replacement (TKR) or in Total Hip Replacement (THR). Accurate models could avoid the standard in-vitro time-consuming investigation procedures (simu­ lators) and should contribute as tool for the optimization of the human prostheses tribological design. Main target, regarding the prosthesis design, is its tribological characterisation since, of course, a high volu­ metric wear rate is undesired. Obviously, wear is a complex and multiscale phenomenon, and its accurate in silico prediction needs sufficiently accurate and validated theoretical models which require deep theoretical and experimental investigations in many fields of knowledge, such as contact mechanics, stress-strain analysis, musculoskeletal multibody modelling, synovial lubrication modelling (boundary/mixed, hydrodynamic and EHD), tribo-corrosion, bio-materials characterizations, etc. The wear phenomenon in a joint is affected by several factors, like

the coupled materials [3], the geometry of the tribopair, the roughness of the contact surfaces, the load conditions [4] and the lubrication of the tribosystem [5–7]. Regarding Total Hip Replacement (THR), main biomaterials actually used are polymers, metals and ceramics. The polymers are mainly the Ultra High Molecular Weight PolyEthylene (UHMWPE) and Ultra High Cross Linked PolyEthylene (UHXLPE) which, despite of their high wear resistance, produce wear debris detected as source of infection by the synovial joint, which tries to eject them causing the osteolysis [8]. The metallic biomaterials, like stainless steel, titanium (Ti) or cobalt-chromium-molybdenum (CoCrMo) alloys, are characterised by high wear resistance, high hardness, toughness and surface finishing, but still now there is no sufficient scientific demonstration about the metal wear debris non-toxicity [9]. The ceramics used in arthroplasty, like zirconia and alumina or composites starting from their combination, are the most performing materials thanks to their high wear resistance, low friction and high wettability, in addition to their excellent biocompatibility [8]. Interesting biomaterials reviews regarding the Total Hip Arthroplasty (THA) are available in literature [8–11]. From tribological point of view, in order to reduce the implant failure and the number of revisions, the researchers studied the various combinations of materials testing the pairings, looking for the optimal configuration in terms of volumetric wear rate, which is an important pre-clinical issue [12–14]. Since one of the main goal of a prosthesis design is its duration,

* Corresponding author. Department of Industrial Engineering (DIIn), University of Salerno, Via Giovanni Paolo II, 132, Fisciano, SA, Italy. E-mail address: [email protected] (A. Ruggiero). URL: https://docenti.unisa.it/004299/home (A. Ruggiero). https://doi.org/10.1016/j.triboint.2019.105993 Received 8 July 2019; Received in revised form 27 September 2019; Accepted 30 September 2019 Available online 7 October 2019 0301-679X/© 2019 Elsevier Ltd. All rights reserved.

A. Ruggiero and A. Sicilia

Tribology International 142 (2020) 105993

summed up in Ref. [20]. It is interesting to analyse the Stokes couple-stress fluid model [21] because of the non-Newtonian behaviour induced by the interaction between the fluid and the macromolecules within it, through the couple-stress parameter. This non-Newtonian fluid model has already been adopted in Ref. [22] to describe the nat­ ural lubrication in the ankle joint. Some authors proposed THA lubri­ cation models based on the Reynolds equation in various conditions: - in Ref. [23] is developed a full numerical solution to the hydrody­ namic lubrication problem of metal-on-metal or ceramic-on-ceramic hip joint replacements under general load and speed conditions during steady walking, obtaining results about the survival of a small but finite film thickness, the increase of the minimum film thickness due to the presence of a dimple on the acetabular cup, especially when it is deep and close to the direction of the resultant load and, finally, the necessity of an elasto-hydrodynamic analysis; - in Ref. [24] is conducted a parametric analysis on the hydrodynamic lubrication of the artificial hip joint and used the Reynolds equation with the Bergmann cycle inputs, showing that the angular motion provides a negligible sliding effect compared to the squeeze one; in particular his results show that the joint reaches early the contact phase, due to the low viscosity of the Newtonian synovial fluid, and that a better lubrication is provided with the increase of the viscosity, the increase of the femoral head radius and the decrease of the radial clearance; - in Ref. [25] is presented an elasto-hydrodynamic lubrication simu­ lation of metal on metal total hip implant, considering both steady state and transient physiological loading and motion conditions of the gait cycle; they focused on the position of the contact area which was located within the cup dimensions, thanks to the anatomical position of the cup; moreover their results showed that the signifi­ cant squeeze action is due to the horizontal anterior-posterior and medial-lateral load components; - in Ref. [26] the Authors performed a multibody modelling to investigate the hip squeaking and the path of the contact point be­ tween the head and the cup of ceramic on ceramic hip implant, considering three-dimensional physiological gait cycle loading and motion conditions and neglecting the lubrication; the contact force is studied with the Hertzian theory, regarding the normal component, and with the modified Coulomb friction law, regarding the tangential component; they concluded that the three-dimensional vibrations of the head within the cup play an important role regarding the squeaking effect and the wear of the artificial hip joints.

Fig. 1. Hip joint reference frame.

it is necessary to achieve detailed models of the lubrication phenomenon occurring in the synovial joint in order to optimize design parameters toward a better lubrication and then a better wear resistance. The lubrication modes which take place in natural human synovial joint arise from complex mechanisms studied by several researchers and described in Ref. [15] and in Ref. [16]. The complexity of the natural synovial joint lubrication is overcome when the joint is replaced with a prosthesis. In this case we can refer to the three classical lubrication modes (boundary, mixed, full film) in dependence of the loading and the motion conditions [17]. The full fluid film lubrication mode is charac­ terised by the complete separation of the articular surfaces by the mean of the synovial fluid, which produces a pressure field able to sustain loads due to surfaces sliding or squeeze motion [15]. The latter mode is able to carry high loads for short time and is typical of fast and intense physical activities like jumping or running [16]. When the joint un­ dergoes both sliding and squeeze actions, the high pressure reached in the cavity becomes able to deform the surfaces increasing the coupling conformity and ensuring the continuity of the surface separation [15]. The boundary lubrication mode occurs when a high load doesn’t allow to keep the surfaces separated. In this case the responsible of the lubrication is a glycoprotein constituent of the synovial fluid, named lubricin, which is absorbed by both the articular surfaces, forming a macromolecular monolayer that is effective in the friction reduction [15]. This lubrication mode is strictly affected by the chemical proper­ ties of the synovial fluid, while the rheological properties, such as the viscosity, are responsible for the full fluid film mode [15]. The hydro­ dynamic lubrication modelling is commonly approached by the Rey­ nolds equation under its classical hypotheses. In the cases of natural lubrication the equation can be sophisticated considering the porosity of the surfaces with the Darcy law, in order to describe the synovial fluid filtration in the articular cartilage [18], the Elasto-Hydrodynamic lubrication mode taking into account the surface deflection field [19], the non-Newtonian behaviour of the lubricant fluid [15,16]. Regarding the rheological behaviour of the synovial fluid, a number of non-Newtonian fluid models have been proposed in literature, like Power Law Model, Cross Model, Herschel-Bulkley Model and others

The final aim of this research is to give a contribution in the tribo­ logical modelling of THR considering the presence of the lubricant (synovial fluid) proposing a novel lubrication model based on the so­ lution of the Reynolds equation during the gait, which takes in to ac­ count the rise of the Hertzian pressure field in the contact area, occurring when the eccentricity gets over the radial clearance. The advantage is the possibility to consider, during the gait, both the boundary lubrication mode and the hydrodynamic one in the calculation of the cup wear rate during one gait. 2. Theoretical background 2.1. Musculoskeletal simulation In this work the hip kinematics data were calculated by using OpenSim®, which is a musculoskeletal multibody software based on the inverse dynamics approach. The aim of a musculoskeletal simulation is to calculate the relative motion and the internal forces between bodies connected by joints. Although the software is able to provide the hip loads too, for the lubrication analysis the more realistic Bergmann in vivo measured hip loads [27] are chosen to be used. In the software data base several human models are available and 2

A. Ruggiero and A. Sicilia

Tribology International 142 (2020) 105993

Fig. 4. Hip loads during the gait cycle.

strikes the ground again, during which there is no contact between foot and ground; the gait cycle period is about 1.2 s. The right hip angles trajectories during the gait cycle were calculated by OpenSim® software with reference to a musculoskeletal lower limb system characterised by 23� of freedom and moved by 92 muscles: the lower limb system was developed by researchers and the experimental data available with the model, like the kinematic ones, are collected from different subjects. The procedure that leads to the system degrees of freedom evolution starting from measurements is called Inverse Ki­ nematics and it is performed by the means of markers applied on a subject: recording the marker position, the calculation of the degrees of freedom that match the measures is achieved through a minimization. The results regarding the right hip during walking are shown in the Fig. 2. The output of the kinematic analysis provides also the angular ve­ locity vector of the femur with respect to the pelvis during the time, shown in the Fig. 3. Regarding the hip joint internal forces during the gait cycle the Bergmann’s hip-on-femur in vivo measured loads [27] are used and reproduced in the Fig. 4. As will be introduced in the next paragraph, the forces acting on the acetabular cub are calculated; in order to adapt the Reynolds equation to a spherical joint by the means of the reference frame xyz shown in the Fig. 1, they are rotated by the femoral neck angle equal to αc . The Bergmann’s hip-on-femur loads are represented by the vector NB , while the femoral head angular velocity vector with respect to the acetabular cup is ωg . The two vectors are defined in the equation (1). � NB ¼ ½ NAP NPD NML �T (1) ωg ¼ ½ ωAA ωIER ωFE �T

Fig. 2. Hip angles trajectories during the gait cycle.

Fig. 3. Hip angular velocity components during the gait cycle.

usable for inverse kinematics, inverse dynamics, static optimization calculations and so on. Each model is characterised by default pro­ portions between the bodies which can be modified and adapted to a particular subject through the scaling tool. The model here used consists in a lower limb biomechanical system composed by three open chains, the back and the two legs, characterised by average weight and proportions. Regarding the hip joint, shown in the Fig. 1, it is a spherical joint, so it allows three relative rotations between the linked bodies, namely pelvis and femur. The three rotations are called Flexion/Extension (FE) in the sagittal plane, Adduction/Abduction (AA) in the frontal plane and Internal/External Rotation (IER) in the horizontal plane; these rotations are performed around three orthogonal axes called respectively Medial/ Lateral (ML), Anterior/Posterior (AP) and Proximal/Distal (PD). The gait cycle is composed by the stance phase, from 0%, when the right heel strikes the ground, to about 60%, when the toes leave the ground, and by the swing phase, from 60% until 100% when the heel

The hip-on-cup loads Nc;g are the opposite of the Bergmann loads NB : Nc;g ¼

NB

(2)

After defining the rotation matrix Rc , in (3), in which αc represents the neck angle, 2 3 1 0 0 Rc ¼ 4 0 cosð αc Þ sinð αc Þ 5 (3) 0 sinð αc Þ cosð αc Þ the rotation of the hip-on-cup loads Nc;g and the relative angular velocity vector ωg can be performed in (4), obtaining Nc and ωc .

3

A. Ruggiero and A. Sicilia

Tribology International 142 (2020) 105993

Fig. 5. Loads and relative motion in the acetabular cup reference frame during the gait cycle.

Fig. 6. Couple-stress flow for different macromolecules size.

8 < Nc ¼ RT Nc;g c : ωc ¼ RT ωg c

Fig. 7. Couple-stress function.

stress parameter η is obtained writing the equation (5) in an alterna­ tive form, (6), in which the couple-stress length l appears: it could be assumed as representative of the characteristic dimension of the hyal­ uronic acid macromolecules within the synovial fluid [22]. � � rffiffiffi ∂ ∂2 u η (6) τ¼μ u l2 2 l ¼ ∂y μ ∂y

(4)

Considering a normal femoral neck angle, αc ¼ 54� [28], the relative motion and load conditions of the hip joint, referred to xyz frame, are shown in the Fig. 5. (see Fig. 6) These vectors will be the definitive inputs of the Reynolds equation.

Executing the dynamical equilibrium of the infinitesimal fluid vol­ ume within the surface gap and by integrating them with appropriate boundary conditions [22] the fluid velocity profile is obtained in (7), it is shown in 2–6:

2.2. Lubrication modelling A modification of the Reynolds equation is performed taking advantage of the non-Newtonian behaviour of the Stokes couple-stress fluid model [21], that is characterised by the relation between the shear stress τ and the shear strain ∂u=∂y shown in the equation (5) [22].

τ¼μ

∂u ∂y

3

η

∂u ∂y3

ui ðyÞ ¼

(5)

The fluid behaviour is affected by the viscosity μ and the couplestress parameter η. A more meaningful interpretation of the couple-

� 1 ∂p hy y2 2μ ∂xi � � �y� 2 3 y h þ cosh cosh 2 l ∂p 6 l l 7 � � 15 4 h μ ∂xi 1 þ cosh l

Ui;1 þ

Ui;2

Ui;1

h

y

(7)

Integrating once again the velocity profile ui , the specific flow qi is 4

A. Ruggiero and A. Sicilia

Tribology International 142 (2020) 105993

� 8 ∂ > > r¼ > > ∂ x > > > > > < V ¼ ½U > > > > > > > > > : f ¼g

∂ ∂ ∂y ∂z 0

W�

�

→r⋅½f rp� ¼ r⋅ðhVÞ þ

∂h ∂t

(12)

h3 12μ

Differentiating the velocity profile according to the couple-stress model, the shear stress τi is obtained in the equation (13).

τi ðyÞ ¼ μ

Zh qi ¼

ui ðyÞ dy ¼

Ui;1 þ Ui;2 h 2

0

6 6 6 3 h ∂p 6 6 61 12μ ∂xi 6 6 6 4

h l 12

Ui;1 h

1 ∂p ðh 2 ∂xi

� � h 3 sinh l � �7 2 h 7 7 1 þ cosh l 7 7 7 � �3 7 h 7 7 l 5 (8)

� � � � � � 8 ∂ ∂p ∂ ∂p ∂ ∂ ∂h > Þ þ Þ þ R sinðθÞ ðsinðθÞhU ðhU sinðθÞ sinðθÞf þ f ¼ R sinðθÞ > θ ϕ > > ∂θ ∂ϕ ∂ϕ ∂θ ∂ϕ ∂θ ∂t > > > > > < ∂p ∂p pð0; ϕÞ ¼ pðπ; ϕÞ ¼ p0 or ð0; ϕÞ ¼ ðπ; ϕÞ ¼ 0 > ∂θ ∂θ > > > > > > > ∂p ∂p > : pðθ; 0Þ ¼ pðθ; πÞ ¼ p0 or ðθ; 0Þ ¼ ðθ; πÞ ¼ 0 ∂ϕ ∂ϕ

g

h3 ∂p 12μ ∂xi

(9)

The dimensionless function g (10) is characteristic of the couplestress fluid non-Newtonian behaviour: it depends on the ratio h= l defined as ξ. ξ gðξÞ ¼ 1

12

sinhðξÞ 2 1 þ coshðξÞ ξ3

(13)

(15)

Transforming the differential operators with the spherical co­ ordinates, is possible to write the modified Reynolds equation, for the spherical joint, in (15), followed by the Dirichlet or Neumann boundary conditions on the cup edges [23]. The fluid film thickness h depends on the radial clearance c and the dimensionless eccentricity n ¼ e=c by the equation (16) [23].

Denoting the arithmetic average surface velocity with Ui and incor­ porating the term in the square bracket in a function g, the specific flow can be written qi ¼ Ui h

2yÞ

Assuming the femoral head as a sphere, with radius r, and the acetabular cup as an hemispherical cavity, with radius R, a reference frame fixed with the cup and located in its centre, with axes x and z lying in a plane tangent to the edge of the cup and axis y forming a right-hand orthonormal reference frame is assumed. The spherical joint system, shown in Fig. 8, demonstrates that the relative motion is due to the eccentricity e time variation and to the angular velocity vector ω of the femoral head with respect to the acetabular cup. In this configuration a point P on the cup is located by the unit vector br , that depends on the spherical angles θ and ϕ, used in (14). 8 < R sinðθÞcosðϕÞ ðP Os Þ ¼ Rbr ¼ R sinðθÞsinðϕÞ (14) : R cosðθÞ

Fig. 8. Spherical joint reference frame.

2

Ui;2

hðθ; ϕ; tÞ ¼ cð1

nðtÞ⋅br ðθ; ϕÞÞ

(16)

The velocity of a point Q on the femoral head, located by the unit vector br , is due to the eccentricity variation with respect to the time e_ and to the angular velocity vector ω ½29�, while the velocity of the point P located by the same unit vector br is null. The two velocities are shown in (17). � vP ¼ 0 (17) _ þ ωðtÞ � ½ðR hðθ; ϕ; tÞÞ ​ br eðtÞ� vQ ðθ; ϕ; tÞ ¼ eðtÞ

(10)

A typical couple-stress function is represented in Fig. 7. Note that for l→0, g→1, so that the fluid returns to behave as a Newtonian fluid. The continuity equation for the infinitesimal fluid volume, allows to obtain the general modified Reynolds equation (11). �� 3 � � �� 3 � � ∂ h ∂p ∂ h ∂p ∂ðUhÞ ∂ðWhÞ ∂h (11) g þ g ¼ þ þ ∂x ∂z 12μ ∂x 12μ ∂z ∂x ∂z ∂t

In order to evaluate the velocity component in the spherical frame, the vector vQ is rotated by the means of a rotation matrix Rs , in (18), that depends on the spherical angles θ and ϕ. Then, performing the arith­ metic average of the velocities of P and Q points, the velocity terms Uθ and Uϕ which appear in the Reynolds equation are obtained, written in the equation (19).

Defining a new function f, the equation (11) takes a more compact form in (12), which easily explains that the pressure field, included in the term on the equation left-hand, arises as a consequence of the sliding effect and the squeeze one, respectively the first term and the second one on the right-hand of the equation. 5

Tribology International 142 (2020) 105993

A. Ruggiero and A. Sicilia

v¼

RTs vQ ¼ ½ Uθ 2

Uϕ

(19)

U ρ �T

At this point the hydrodynamic problem is entirely characterised. 2.3. Contact modelling In order to consider the occurring of the boundary lubrication mode, defined by the contact area, the Hertz contact theory is taken into ac­ count [30,31]. In the particular case of a ball and socket joint, the Hertz theory provides a constitutive non-linear elastic relation between the contact force Fc and the penetration depth δ of the two surfaces [26]. The stiffness k of the relation depends on the relative curvature between ball and socket and their mechanical properties. The Hertz model is shown in the equation (20). 8 1 1 1 > > ¼ > > R r R > eq > > > > > > > < 3 1 1 ν2b 1 ν2s (20) →Fc ¼ kδ2 ¼ þ * > E E E > b s > > > > > > pffiffiffiffiffiffi > > > 4E* Req > : k¼ 3

Fig. 9. Contact area reference frame.

Since the aim is to develop a hydrodynamic-boundary lubrication model, it is necessary to know the location of the contact area centre rC on the spherical joint reference frame, in terms of spherical angles θc and ϕc : the contact area rises where the eccentricity vector norm gets over the radial clearance, so the dimensionless eccentricity vector norm jnj gets over the unity. Moreover it is assumed that, in this position, the penetration depth δ is represented by the negative of the film thickness. The contact reference frame is shown in the Fig. 9. The calculation of the contact area centre location and the penetration depth are described in the equation (21).

Table 1 Input parameters regarding a CoP pair.

Fig. 10. Spherical angles grid domain.

2 8 cosðθÞ > > > Ry ðθÞ ¼ 4 0 > > < 2 sinðθÞ cosðϕÞ > > > > Rz ðϕÞ ¼ 4 sinðϕÞ > : 0

3 0 sinðθÞ 1 0 5 0 cosðθÞ 3 →Rs ðθ; ϕÞ ¼ Rz ðϕÞRy ðθÞ sinðϕÞ 0 cosðϕÞ 0 5 0 1

Head radius

r

14 ​ mm ½25�

Radial clearance

c

50 ​ μm ½26�

Fluid viscosity

μ

1:5⋅10

Dimensionless couple-stress length

l=c

0:1 ​ ½22�

Cup Poisson coefficient

νc

0:47 ​ ½38�

νh

0:23 ½26�

kw

10

μf

0:03 [14]

Cup Young modulus Head Poisson coefficient

(18)

Head Young modulus Wear factor Friction coefficient Boundary pressure

Fig. 11. From loads to eccentricities workflow. 6

Ec Eh

p0

3

​ Pa ​ s ½24�

1 GPa ½6� 358 ​ GPa ½26� 7

​ mm3 =Nm ½38�

0 ​ Pa ½24�

A. Ruggiero and A. Sicilia

Tribology International 142 (2020) 105993

Fig. 12. Eccentricity vector and its time derivative during the gait cycle.

0qffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 8 n2x þ n2y > > A θc ¼ arctan@ > > ( > nz < rC ¼ R br ðθc ; ϕc Þ jnj > 1→ → > > δ ¼ cðjnj 1Þ > � � > > : ny ϕc ¼ arctan nx

(21)

Once the penetration depth δ is known, the radius of the circular contact area a can be calculated, followed by the mean contact pressure pm and the maximum one pmax , (22). 8 > 3 > > kδ2 > > ¼ p > m < pffiffiffiffiffiffiffiffi Ac a ¼ Req δ→Ac ¼ πa2 → (22) > > > 3 > > > : pmax ¼ 2pm A point P in the contact area is located by the vector rP , calculated starting from the centre location rC and its orientation, resulting in the rotation matrix Rs;c , that is Rs with angles θc and ϕc . The vector rP exists in the domain’s area characterised by negative film thickness h (23). hðθ; ϕÞ < 0→rP ðθ; ϕÞ ¼ rC þ Rs;c uC

Fig. 13. Dimensionless eccentricity norm during the gait cycle.

(23)

The ellipsoidal Hertz pressure field can be applied over the contact area once the transfer in the contact reference frame Cxc yc is done. So the local position vector of the point P, uC , is needed, in (24). uC ¼ RTs;c ðrP

rC Þ ¼ ½ uc;θ

uc;ϕ

uc;ρ �T

(24)

The uC vector entries are used to evaluate the Hertzian pressure field by the equation (25). sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u2c;θ þ u2c;ϕ (25) pHz ðθ; ϕÞ ¼ pmax 1 a2 At this point the hydrodynamic-boundary problem is entirely char­ acterised, so the calculation of the load and the friction torque vectors is possible. The load vector is easily calculated by integrating the pressure field on the socket surface, as shown in (26). ZZ NðtÞ ¼ pðθ; ϕ; tÞbr R2 sinðθÞdθdϕ#ð26Þ (26) Regarding the friction torque, it is necessary to distinguish the fullfilm lubricated areas by the contact ones. In the first case, the shear stresses are calculated by the equation (27) and organised in the shear stress vector defined in the spherical frame, ðτ w Þs (28).

Fig. 14. [23] hip reference frame.

7

A. Ruggiero and A. Sicilia

Tribology International 142 (2020) 105993

Fig. 15. (a) Calculated eccentricities in hip reference frame; (b) [25] MoM EHL eccentricities in hip reference frame.

8 > > > <

τw;θ ¼ τθ ðy ¼ hÞ ¼ μ

2.4. Wear modelling

vθ h ∂p þ h 2R ∂θ

(27)

> > v h ∂p > : τw;ϕ ¼ τϕ ðy ¼ hÞ ¼ μ ϕ þ h 2R sinðθÞ ∂ϕ ðτ w Þs ¼ ½ τw;θ

Wear calculation was approached by means of Archard wear model [31]. It assumes that the wear volume V is proportional to the normal contact load Fc and to the sliding distance s through a constant kw , wear factor, which depends on the coupling materials, as shown in (33).

(28)

τw;ϕ 0 �T

_ þ ωðtÞ � ðR ​ br vQ ðθ; ϕ; tÞ ¼ eðtÞ

The sliding distance is calculated starting from the sliding velocity vc of the contact area centre, in (34). The vector vc is rotated by the means of Rs;c to evaluate its components in the spherical frame.

(29)

eðtÞÞ

_ þ ωðtÞ � ðrC vc ðtÞ ¼ RTs;c ½eðtÞ

vQ;t eðtÞÞ→bt ¼ �� �� vQ;t

dV ¼ kw Fc ðtÞds ¼ kw Fc ðtÞvsd ðtÞdt

vc;ρ �T

(34)

(36)

Integrating over the time the equation (36), the time trend of the wear volume VðtÞ is obtained, in (37). Zt

(31)

μf pHz RTsbt

vc;ϕ

The infinitesimal wear volume during a time dt is written in the equation (36).

(30)

The shear stress vector in the spherical frame in the contact area is given by the equation (31), rotating the unit vector by the way of the rotation matrix Rs . The friction coefficient μf is determined by the ma­ terials coupling. ðτ w Þs ¼

eðtÞÞ� ¼ ½ vc;θ

The sliding velocity vsd is assumed to be the norm of vc in the contact plane, as shown in the equation (35). qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ds (35) vsd ðtÞ ¼ v2c;θ þ v2c;ϕ ¼ dt

Since the time derivative of the eccentricity vector e_ represents a normal velocity in the contact area, while the other term refers to a tangential velocity vQ;t , the tangential unit vector bt is calculated in (30) by the means of the norm of the latter term. vQ;t ¼ ωðtÞ � ðR ​ br

(33)

V ¼ k w Fc s

In the second case, a simple Coulomb model is proposed, based on the knowledge of the tangential unit vector bt in every point of the domain. A generic point Q in the contact area slides with velocity vQ , in the equation (29).

VðtÞ ¼ kw

s

The knowledge of the shear stress vector ðτ w Þ over the whole socket surface domain leads to the calculation of the friction torque vector, given in the equation (32). ZZ Mf ðtÞ ¼ ðRbr � Rs ðτ w Þs ÞR2 sinðθÞdθdϕ (32)

Fc ðtÞvsd ðtÞdt

(37)

0

3. Fluid pressure field calculation method The modified Reynolds equation (15) is a Partial Differential Equa­ tion (PDE). The pressure doesn’t depend strictly on the time, but is related to it through the time dependent eccentricity and angular ve­ locity. The pressure can be found solving the elliptic problem for each time step using the finite difference method [32]. In this work, the model is developed in Matlab® environment. The method consists into discretize the spatial domain, represented in this case by the femoral cup spherical angles, and by writing the equation for each grid point (Fig. 10).

The Reynolds equation furnish the pressure field pðθ; ϕ; tÞ, over the time t, due to the relative motion, identified by ωðtÞ, and the eccen­ tricity, identified by the dimensionless eccentricity nðtÞ. Then one can calculate the load NðtÞ and the friction torque Mf ðtÞ.

8

A. Ruggiero and A. Sicilia

Tribology International 142 (2020) 105993

Fig. 16. Pressure (a) and film thickness (b) fields followed by pressure spherical plot (c) at 59.7% of gait.

In the equation (38) i and j refer to the spherical angles spatial domain, respectively θ and ϕ, while n refers to the time domain t. The Reynolds equation can be written as:

The aim is to write a linear invertible system to found the pressure pnij

in each grid point, using the finite difference: in particular the central

� � � � � � �� �� ∂ ∂p ∂ ∂p n ∂ ∂ ∂h n ¼ R sinðθÞ sinðθÞ sinðθÞf þ f ðsinðθÞhUθ Þ þ ðhUϕ Þ þ R sinðθÞ ∂θ ∂ϕ ∂ϕ ij ∂θ ∂ϕ ∂θ ∂t ij

9

(38)

A. Ruggiero and A. Sicilia

Tribology International 142 (2020) 105993

Fig. 17. Pressure and film thickness fields followed by pressure spherical plot at 84.5% of gait.

differences, for the first and second spatial derivatives, and the forward difference, for the time derivative, are chosen. Starting from the first term of the left-hand of the equation, it assumes the form in (39).

� �� � � � � þ � þ � � ∂ ∂p n dþ dm ⋅sinðθÞ dm ⋅f θ dθ ⋅pθ m ⋅sinðθÞ dm ⋅f θ dθ ⋅pθ sinðθÞf sinðθÞ ¼ sinðθ iÞ ∂θ ∂θ ij Δθ

10

(39)

A. Ruggiero and A. Sicilia

Tribology International 142 (2020) 105993

The vectors in (39) are: � � 8 pθ ¼ pni 1;j pnij pniþ1;j T > > > > > > > > > θ ¼ ½ θi 1 θi θiþ1 �T > > > > > > > > > > ½ 0 1 1 �T > þ > > d ¼ > m > 2 > > > > > > > > > < ½ 1 1 0 �T dm ¼ 2 > > > > > > > > > ½0 1 1 �T > þ > > > dθ ¼ > Δθ > > > > > > > > > > ½ 1 1 0 �T > > d ¼ > θ > Δθ > > > > > : �T � n n n f θ ¼ f i 1;j f ij f iþ1;j

Fig. 18. Friction torque vector during the gait cycle.

(40)

Combining the vectors in (40), the matrices D2θ and Fθ are built, in (42) and (41). 8 þ þT > > Dþ > m ¼ d m dm < (41) Dm ¼ d m d m T > > > : F ¼ sinðθ ÞsinðθÞf T θ

i

θ

8� � � þ� þ þ > < Dθ :;:;k ¼ Dm dθ k � � � � →D2θ ¼ > : D θ :;:;k ¼ Dm dθ k

Dþ θ

Dθ Δθ

(42)

The vector aθ defined as in (43) allows to write the first term of the left-hand side of the Reynolds equation as in (44). aθ ¼ D2θ Fθ

(43)

� � �� ∂ ∂p n sinðθÞ ¼ aθ ⋅pθ sinðθÞf ∂θ ∂θ ij

(44)

Following the same procedure, the discretization of the second term of the left-hand side of the equation is performed, in (45), defining new vectors reported in (46). � þ � � � � � �� dþ dm ⋅f ϕ dϕ ⋅pϕ ∂ ∂p n m ⋅f ϕ dϕ ⋅pϕ (45) ¼ f ∂ϕ ∂ϕ ij Δϕ

Fig. 19. Maximum pressure and minimum film thickness during the gait cycle.

8 p ¼ � pn i;j 1 ϕ > > > > > > > > > > ½0 > > dþ > ϕ ¼ > > < > > > > > > > > > > > > > > > :

dϕ ¼

½ 1

� n f ϕ ¼ f i;j

1

pnij

pni;jþ1

�T

1 1 �T Δϕ 1 0� Δϕ f nij

(46)

T

f ni;jþ1

�T

The combination of the vectors in (46) leads to the matrix D2ϕ , in (47). 8 þ þT þ > < Dϕ ¼ d ϕ d m Dþ Dϕ ϕ →D2ϕ ¼ (47) > Δϕ :D ¼ d d T ϕ m ϕ Fig. 20. Wear volume during the gait cycle. 11

A. Ruggiero and A. Sicilia

Tribology International 142 (2020) 105993

Fig. 21. Eccentricities characterised by high squeeze velocitiy.

Fig. 22. Fluid pressure field due to significant squeeze action.

Then the second term is writable as a scalar product, reported in (49), between the vector aϕ , defined in (48), and the pressure vector along ϕ direction pϕ . (48)

aϕ ¼ D2ϕ f ϕ �

�

��n

∂ ∂p f ∂ϕ ∂ϕ

��n � � � � ∂ ∂ ðsinðθÞhUθ Þ þ ðhUϕ Þ ¼ R sinðθi Þ dcθ ⋅uθ þ dcϕ ⋅uϕ R sinðθÞ ∂θ ∂ϕ ij

¼ aϕ ⋅pϕ ij

(49)

The sliding term on right-hand side of the equation is written as in (50), with vectors defined in (51).

12

(50)

A. Ruggiero and A. Sicilia

Tribology International 142 (2020) 105993

8 <

� ps ¼ pni

pni;j

1;j

: a ¼ ½ ½a � θ 1 I

½aϕ �1

1

pnij

pni;jþ1

½aθ �2 þ ½aϕ �2

pniþ1;j ½aϕ �3

�T ½aθ �3 �

→aI ⋅ps ¼ bI

(56)

In order to achieve a well stated problem, it remains to add the Dirichlet or Neumann discretized boundary conditions, reaching the linear system in (57). The column vector p is obtained scanning for each i the entire j domain. 8 n pn2j pn1j pnnθ þ1;j pnnθ ;j > p1j ¼ pnnθ þ1;j ¼ p0 or ¼ ¼0 > > Δθ Δθ > > > > > < aI ⋅ps ¼ bI > > > > > > n n > > n pn pni1 pi;nϕ þ1 pi;nϕ : pi1 ¼ pni;nϕ þ1 ¼ p0 or i2 ¼ ¼0 Δϕ Δϕ

The linear system in its matrix form is written in (58), and can be inverted for each time step in order to evaluate the pressure field evo­ lution over time.

Fig. 23. Load capacity due to significant squeeze action.

8 > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > <

dcθ ¼

dþ θ þ dθ 2

dcϕ ¼

dþ ϕ þ dϕ 2

� n hθ ¼ hi � n hϕ ¼ hi;j

1;j

hnij

hniþ1;j

1

hnij

hni;jþ1

�T �T

The solving of the Reynolds equation requires the cavitation constraint, that doesn’t allow the negative pressure rise. Then an itera­ tive method is used, instead of the classical system inversion, to find the vector p. The iterative methods, regarding the elliptic equations, consist in building the iteration matrix G, obtained by splitting A into two other matrices, M and N in (60), in order to find the solution vector p at the ðk þ1Þ iteration starting from the one at the ðkÞ iteration. The iterative cycle stops when the convergence is reached, as described in (61).

(51)

A¼M �

uϕ ¼ hϕ V ϕ The squeeze term on right-hand side of the equation is written as in (52), with vectors defined in (53). � � �� ∂h n R sinðθÞ R sinðθÞ ¼ R sinðθi Þ½R sinðθi Þðdt ⋅ht Þ� (52) ∂t ij

> > > :

h ht ¼

hnij 1

hnij

iT

N→Mpðkþ1Þ ¼ NpðkÞ þ b

�� � G ¼ M 1 N while ​ ​ max �pðkþ1Þ pðkÞ � < tol ​ → 1 ~ b¼M b pðkþ1Þ ¼ GpðkÞ þ ~ b

(60) (61)

The iterative method allows to add the cavitation constraint, described in (62), at each iteration step. � � � � (62) if pðkþ1Þ I < 0 then pðkþ1Þ I ¼ 0 The iterative method used in this work is the Successive Over Relaxation one, SOR. It is based on the splitting of A in its triangular parts, lower L and upper U, and diagonal part D, as shown in (63), used to build M and N matrices in (64).

(53)

(63)

A¼LþDþU

Finally the whole right-hand side of the equation can be denoted as a constant bI , in (54). � � R sinðθi Þ dcθ ⋅uθ þ dcϕ ⋅uϕ þ R sinðθi Þðdt ⋅ht Þ ¼ bI (54)

8 > > > <

MSOR ¼

> > > : NSOR ¼ ð1

At this point, the whole equation has assumed the form in (55). aθ ⋅pθ þ aϕ ⋅pϕ ¼ bI

(58)

Ap ¼ b

The computational model allows to evaluate the contact phase modifying the rows of A and the components of b in the grid points where the film thickness is negative, as described in (59). 8 < aI ¼ ½ 0 0 1 0 0 � n if hij < 0→ (59) : bI ¼ pnHzij

> > > > > �T � > > > V θ ¼ Uθ ni 1;j Uθ nij U nθiþ1;j > > > > > > > > > �T � n > n n > > > V ϕ ¼ U ϕi;j 1 Uϕ ij U ϕi;jþ1 > > > > > > > > > uθ ¼ sinðθÞhθ V θ > > > > :

8 ½ 1 1 �T > > d ¼ > < t Δt

(57)

D þ ωL

ω ωÞD ω

ωU

(64)

The SOR method takes advantage of the convergence parameter ω tuned on its optimal value, described in (65), that depends on the grid resolution. vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u 2 2 u →ωopt ¼ h¼u 1 (65) t 1 1 þ sinðπhÞ þ Δθ2 Δϕ2

(55)

Analysing the vectors pθ and pϕ , the cross stencil vector ps is built (Fig. 10) then the vectors aθ and aϕ are rearranged in a new vector aI , defining a new scalar product, in (56).

13

A. Ruggiero and A. Sicilia

Tribology International 142 (2020) 105993

Once the equation is solved, the load and the friction torque vectors, in (66), are calculated with the trapezoidal integrals, while the wear volume, in (67), is evaluated through a summation. ZZ ZZ 8 NðtÞ ¼ pðθ; ϕ; tÞbr R2 sinðθÞdθdϕ ¼ f N dθdϕ < ZZ ZZ (66) : Mf ðtÞ ¼ ðRbr � Rs ðτ w Þs ÞR2 sinðθÞdθdϕ ¼ f Mf dθdϕ Zt

Zt Fc ðtÞvsd ðtÞdt ¼

VðtÞ ¼ kw

The geometrical, rheological and mechanical parameters, together with the coupling constants and the Dirichlet boundary pressure, listed in the Table 1, are taken from available literature; the load and relative motion conditions are shown in Fig. 5; the SOR iterative method is used to solve the modified Reynolds equation. In the Fig. 12 the plot of the dimensionless eccentricity vector, fol­ lowed by its time derivatives, and its norm in 4–2, are shown (see Fig. 13). The dimensionless eccentricity vector norm plot points out that the joint is always in contact, since it is greater than the unity for the entire cycle. In particular the y component of the dimensionless eccentricity vector is the closest to the unity and, together with the z component, concurs to the contact occurring, while the x component has the lowest absolute value, as expected. A comparison on the eccentricity calculation with [25] is shown in Fig. 15, despite the term of comparison is based on a Metal-on-Metal (MoM) prosthesis elasto-hydrodynamic lubrication model, with c ¼ 30 μm . The results are interesting for the qualitative match of the curve trends, written in Ref. [25] reference frame shown in Fig. 14, consid­ ering that the radial clearance actually used in this paper is c ¼ 50 μm . Then, the surf plots of the pressure and film thickness fields at about 59,7% of the cycle, followed by the spherical plot of the pressure, are reported in Fig. 16. In this phase, the contact area is located approximately in θ ¼ π=3 and ϕ ¼ π=2, with a certain extension and the maximum pressure reached is about 4 MPa. As the spherical plot demonstrates, the contact area location points towards the vertical direction of the hip joint, namely the PD direction, where the load is higher. In the surf plot, the transition zone, from the full film lubricated area to the contact one, is visualised: around the Hertzian pressure field the fluid pressure in­ creases because of the smaller film thickness. The same plots during the least critical swing instant, about 84.5% of the cycle, in Fig. 17, are shown. Going towards the cycle, during the swing phase, the contact area location moves to the cup centre, decreasing its dimension, approxi­ mately in θ ¼ π=2 and ϕ ¼ π =2, with a lower pressure peak, of about 2,6 MPa, due to the load decrease. The friction torque, calculated with the resulting pressure field evolution, is shown in Fig. 18. The time trend of the minimum film thickness and the maximum pressure experienced in the pair are reported in the Fig. 19. The maximum pressure experienced in the coupling during the cycle follows the vertical load profile, as confirmed in Refs. [23,29,39], while the minimum film thickness is set to a threshold value during the contact phase, chosen in order to not run into numerical indeterminacy resulting in an undetermined pressure value, characteristic of the modified Rey­ nolds equation, occurring when the film thickness goes to zero, namely the pressure in the joint is undefined if the gap doesn’t exist, regarding a lubricated environment. The threshold value is chosen by imposing the lambda ratio, char­ acteristic of the lubrication mode, lower than the hydrodynamic value, for example equal to 2 [24].

0

fVw ðtÞdt

(67)

0

The explained resolution method is able to calculate the pressure field evolution starting from the eccentricity time variation and the angular velocity vector. In order to evaluate the eccentricity evolution due to a certain load condition, with the same relative motion, an additional iterative cycle on the eccentricity is necessary, based on the Newton method, that sets to zero the difference FðnÞ between the loads calculated with the trial eccentricity, NðnÞ, and the reference loads Nref , shown in (68). FðnÞ ¼ NðnÞ

Nref ¼ 0

(68)

The Newton method is based on the building of the Jacobian matrix, obtained by developing FðnÞ through the Taylor’s series stopped at the first order, like in (69). ~ þ Δn→FðnÞ ¼ Fð~ n¼n n þ ΔnÞ ffi Fð~ nÞ þ JF ð~ nÞΔn ¼ 0

(69)

The Jacobian matrix entries are calculated choosing the small value of ε for the finite differences associated to the derivatives of the F components with respect to the n ones. It is shown in (70). � � Fi nj ð1 þ εÞ Fi n j ∂Fi (70) ½JF ð~ nÞ�ij ¼ ffi ∂nj εnj At this point the vector quantity Δn is calculated, with which the ~ is done, until the maximum of the ab­ update of the trial eccentricity n solute value of F is dropped below a known tolerance quantity, as described in (71). while maxðjFð~ nÞjÞ < tol Δn ¼ ½JF ð~ nÞ� 1 Fð~ nÞ ~¼n ~ þ Δn n

(71)

The workflow of the Newton iterative method is shown in the Fig. 11. 4. Results and discussion In this paragraph, the simulation results regarding the lubrication of a Ceramic-on-Polyethylene prosthesis (CoP), in particular made of Alumina ceramic femoral head against UHXLPE acetabular cup, during the gait cycle, are shown. Moreover, some comparisons with respect to other authors are performed. A CoP pair is chosen because of its good agreement between dura­ tion, fixation, stability and inexpensiveness with respect to others implant choices, regarding the THA. The Soft-on-Hard main problem is represented by the wear of the soft material (polyethylene) because of the contact with the hard material (ceramic). In fact, despite of the good quality of the implant, the CoP pair is subject to a high number of re­ visions, as shown in Ref. [33], so the modelling of its tribological properties is interesting from a research point of view, in order to optimize its performance. Since the polymer is considered to have a viscoelastic behaviour, an adjustment of the Hertz elastic theory should be done, for example considering a damping contribution in the calcu­ lation of the contact force proportional to the penetration velocity in order to take into account the energy dissipation characteristic of the material viscoelastic behaviour, as performed in Ref. [26]: anyway many authors used the Hertz non-linear constitutive relation to analyse the polymers deformation [34–37].

λ¼

h ¼2 Ra

(72)

The λ chosen value coupled with roughness information about the materials pair [24], calculated in (73), leads to a threshold value on the ratio film thickness over radial clearance, determined in (74). � qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Ra;head ¼ 0:01 μm (73) →Ra ¼ R2a;head þ R2a;cup ¼ 2:5 μm Ra;cup ¼ 2:5 μm � � h � � c lim h Ra λ λ¼ ¼ 0:1 → ¼ c lim c Ra c

14

(74)

A. Ruggiero and A. Sicilia

Tribology International 142 (2020) 105993

The model set the film thickness to the threshold value with logic explained in (75). � � � � h h if hnij < c then hnij ¼ c (75) c lim c lim

5. Conclusions The aim of this work was to accurately model the tribological behaviour of the artificial hip joint prosthesis during walking accounting for unsteady synovial lubrication phenomena for wear determination. In order to calculate the volumetric wear rate of a hip prosthesis during the gait cycle, an hydrodynamic-boundary lubrication model of the spherical joint was proposed and developed in Matlab® environ­ ment. The simulations were performed considering a walking physical activity described by Bergmann in vivo load measurements and by the hip angles trajectories obtained by OpenSim® musculoskeletal multi­ body software in correspondence of a normal gait cycle. The developed model allows the tribological simulation of a generic THR pair by varying the geometrical, rheological and mechanical parameters of the tribo-system; moreover the proposed model allows to optimize the prosthesis tribological design by selecting the configuration charac­ terised by the minimum volumetric wear rate during the time for improving its longevity, also by considering several kind of human ki­ nematics. The model was validated with reference to a 28 mm CoP (Alumina vs UHXLPE) pair allowing to evaluate faithfully the volumetric wear rate in the joint and achieving a good agreement with in vitro measurements performed by other authors: in particular, a volumetric

The wear volume during the gait cycle is shown in Fig. 20: is char­ acterised by an initial increase due to the contact occurring, followed by a rapid growth, namely a higher slope, during the stance phase, while the swing phase shows a lower slope, as expected. Moreover, the wear plot shows that the volume of lost material in one gait cycle is Vwear;f , in (76). Vwear;f ¼ 2⋅10

15

m3 cyc

(76)

With the hypothesis of linearity between the wear volume and the number of cycles, considering that about 106 cycles are executed in one year by an average man and assuming that the gap geometry doesn’t change considerably during the cycles despite the loss of material, the volumetric wear rate can be evaluated as in (77). V_ wear ¼ Vwear;f ncyc ¼ 2⋅10

15

m3 cyc mm3 ¼2 106 year cyc year

(77)

3

wear rate of 2 mm year was obtained, of the same magnitude order of those of

The volumetric wear rate calculated in (77) is comparable, from a magnitude order point of view, with the ones found in literature reports 3

aware of 6,59 mm [38] and 3,35 � 0,29 year

mm3 year

3

3

mm 6; 59 mm year ½38� and 3; 35 year ½9� obtained, the first by using a FEM model

and the second one in vitro by using simulators. Of course the proposed model could be more and more improved for a more accurate in silico wear assessment. In particular future re­ searches are needed devoted mainly to:

½9�.

The resulting squeeze velocities in Figs. 4–1, obtained by the using of the finite difference for the dimensionless eccentricity time derivatives, don’t produce the pressure required to keep the surfaces totally sepa­ rated, during the gait cycle load and motion conditions: the full sepa­ ration could be established during others physical activities, like running or jumping, characterised by higher values of squeeze velocities [16]. However, since the developed model is able to distinguish over the socket surface domain the lubricated areas by the contact ones, if a significant squeeze action occurs, the fluid pressure field can reach higher value in order to keep the surfaces separated and bearing very high loads even if for short times. A demonstration of this capacity is given considering, with the same input parameters, the eccentricities shown in the Fig. 21, characterised by a significant time slope without contact, coupled with a constant angular velocity on the z axis, ωz ¼ 2rad=s. The resulting pressure field, in pure full-film lubrication, shown in the Fig. 22, reaches a peak value of about 80 MPa, leading to a high load capacity in Fig. 23, demonstrated by the y component load value of about 2,7 kN.

- improve the mathematical description of the contact in a mixed lubrication regime; - remove some classical Reynolds hypotheses, as rigid surfaces, in order to take into account the elasto-hydrodynamic lubrication mode, with various deflection models switchable with the boundary contact phase, or smooth surfaces, in order to take into account the roughness together with its role in the wear phenomenon and its variation related to the progressive wear; - establish a linking between the lubrication model and a FEM one to take advantage of the calculated deformation field by taking into account the more detailed mechanical behaviour of UHXLPE and by simulating other types of THR like CoC and MoM pairs; - use other types of non-Newtonian fluid model, to take into account the synovial fluid rheological behaviour in a more accurate way. Acknowledgements This research was supported by MIUR PRIN 2017 BIONIC.

List of Symbols N

Load Angular velocity R Rotation matrix τ Shear stress μ Viscosity η Couple-stress parameter l Couple-stress length h Fluid film thickness p Pressure x; y; z Cartesian coordinates U; V; W Velocity components q Specific flow gðξÞ Couple-stress function t Time

ω

15

A. Ruggiero and A. Sicilia

r R c br θ; ϕ e n

ν

E δ a Mf s Fc kw V V_ wear

i; j; n d; D A; b M; N JF

Tribology International 142 (2020) 105993

Head radius Cup radius Radial clearance Radial unit vector Spherical angles Eccentricity Dimensionless eccentricity Poisson coefficient Young modulus Penetration depth Contact area radius Friction torque Sliding distance Contact force Wear factor Wear volume Volumetric Wear rate Discretization indices Finite difference vectors and matrices Linear system matrix and known term vector Matrices forming the iteration matrix Jacobian matrix for Newton cycle

References

[19] Wang W, Wang F, Jin Z, Dowson D, Hu Y. Numerical lubrication simulation of metal-on-metal artificial hip joint replacements: ball-in-socket model and ball-onplane model. Proc Inst Mech Eng J J Eng Tribol 2009;223(7):1073–82. [20] Chhabra R. Non-Newtonian fluids: an introduction. In: Rheology of complex fluids. New York: Springer; 2010. p. 3–34. [21] Stokes VK. Couple stresses in fluids. Phys Fluids 1966;9(9):1709–15. [22] Ruggiero A, G� omez E, D’Amato R. Approximate closed-form solution of the synovial fluid film force in the human ankle joint with non-Newtonian lubricant. Tribol Int 2013;57:156–61. [23] Jin ZM, Dowson D. A full numerical analysis of hydrodynamic lubrication in artificial hip joint replacements constructed from hard materials. Proc Inst Mech Eng C J Mech Eng Sci 1999;213(4):355–70. [24] Ramjee S. Numerical analysis of lubrication in an artificial hip joint. 2008. [25] 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(6):720–9. [26] Askari E, Flores P, Dabirrahmani D, Appleyard R. Nonlinear vibration and dynamics of ceramic on ceramic artiflcial hip joints: a spatial multibody modelling. Nonlinear Dyn 2014;76(2):1365–77. [27] Bergmann G, Bender A, Dymke J, Duda G, Damm P. Standardized loads acting in hip implants. PLoS One 2016;11(5):e0155612. [28] Byrne DP, Mulhall KJ, Baker JF. Anatomy & biomechanics of the hip. Open Sport Med J 2010;4(1). [29] Askari E, Andersen M. A modification on velocity terms of Reynolds equation in a spherical coordinate system. Tribol Int 2019;131:15–23. [30] Fischer-Cripps AC, Gloyna EF, Hart WH. Introduction to contact mechanics. New York: Springer; 2000. [31] Popov VL, Hefl M, Willert E. Handbook of contact mechanics. Springer; 2019. [32] LeVeque RJ. Finite difference methods for ordinary and partial differential equations: steady-state and time-dependent problems. Washington: Siam: Seattle; 2007. [33] RIAP. Progetto registro italiano ArtroProtesi, quarto report. 2017. [34] Wang A, Sun DC, Stark C, Dumbleton JH. Wear mechanisms of UHMWPE in total joint replacements. Wear 1995;181:241–9. [35] Van Citters DW, Kennedy FE, Collier JP. Rolling sliding wear of UHMWPE for knee bearing applications. Wear 2007;263(7–12):1087–94. [36] Guezmil M, Bensalah W, Mezlini S. Tribological behavior of UHMWPE against TiAl6V4 and CoCr28Mo alloys under dry and lubricated conditions. J Mech Behav Biomed Mater 2016;63:375–85. [37] Gevaert MR, LaBerge M, Gordon JM, DesJardins JD. The quantification of physiologically relevant cross-shear wear phenomena on orthopaedic bearing materials using the MAX-shear wear testing system. J Tribol 2005;127(4):740–9. [38] Hung J-P, Wu JS-S. A comparative study on wear behavior of hip prosthesis by finite element simulation. Biomed Eng: Appl Basis Commun 2002;14(04):139–48. [39] Harun MN, Wang FC, Jin ZM, Fisher J. “Long-term contact-coupled wear prediction for metal-on-metal total hip joint replacement,” Proceedings of the Institution of Mechanical Engineers. Participant J: J Eng Tribol 2009;223(7): 993–1001.

[1] Fisher J, Dowson D. Tribology of total artificial joints. Proc Inst Mech Eng H J Eng Med 1991;205(2):73–9. [2] Strickland MA, Taylor M. “"In-silico wear prediction for knee replacements—methodology and corroboration. J Biomech 2009;42(10):1469–74. [3] Affatato S, Ruggiero A, De Mattia JS, Taddei P. Does metal transfer affect the tribological behaviour of femoral heads? Roughness and phase transformation analyses on retrieved zirconia and Biolox® Delta composites. Compos B Eng 2016; 92:290–8. [4] Ruggiero A, Merola M, Affatato S. Finite element simulations of hard-on-soft hip joint prosthesis accounting for dynamic loads calculated from a musculoskeletal model during walking. Materials 2018;11(4):574. [5] Wang A, Essner A, Polineni VK, Stark C, Dumbleton JH. Lubrication and wear of ultra-high molecular weight polyethylene in total joint replacements. Tribol Int 1998;31(1–3):17–33. [6] Mattei L, Di Puccio F, Piccigallo B, Ciulli E. Lubrication and wear modelling of artificial hip joints: a review. Tribol Int 2011;44(5):532–49. [7] Unsworth A. Recent developments in the tribology of artificial joints. Tribol Int 1995;28(7):485–95. [8] Affatato S, Ruggiero A, Merola M. Advanced biomaterials in hip joint arthroplasty. A review on polymer and ceramics composites as alternative bearings. Compos B Eng 2015;83:276–83. [9] Merola M, Affatato S. Materials for hip prostheses: a review of wear and loading considerations. Materials 2019;12(3):495. [10] Zivic F, Affatato S, Trajanovic M, Schnabelrauch M, Grujovic N, Choy KL. Biomaterials in clinical practice: advances in clinical research and medical devices. Springer; 2017. [11] Aherwar A, Singh A, Patnaik A. Current and future biocompatibility aspects of biomaterials for hip prosthesis. AIMS Bioeng 2015;3(1):23–43. [12] Affatato S, Ruggiero A, Jaber S, Merola M, Bracco P. Wear behaviours and oxidation effects on different uhmwpe acetabular cups using a hip joint simulator. Materials 2018;11(3):433. [13] Ruggiero A, D’Amato R, G� omez E. Experimental analysis of tribological behavior of UHMWPE against AISI420C and against TiAl6V4 alloy under dry and lubricated conditions. Tribol Int 2015;92:154–61. [14] Ruggiero A, G� omez E, Merola M. Experimental comparison on tribological pairs UHMWPE/TIAL6V4 alloy, UHMWPE/AISI316L austenitic stainless and UHMWPE/ AL2O3 ceramic, under dry and lubricated conditions. Tribol Int 2016;96:349–60. [15] Nordin M, Frankel VH. Basic biomechanics of the musculoskeletal system. Lippincott Williams & Wilkins; 2001. [16] Dumbleton JH. Tribology of natural and artificial joints. Elsevier; 1981. [17] Jin ZM, Medley JB, Dowson D. Fluid film lubrication in artificial hip joints. Tribol Ser 2003;41:237–56. [18] Ruggiero A, G� omez E, D’Amato R. Approximate analytical model for the squeezefilm lubrication of the human ankle joint with synovial fluid filtrated by articular cartilage. Tribol Lett 2011;41(2):337–43.

16