Jtam.dvi


JOURNAL OF THEORETICAL

AND APPLIED MECHANICS

3, 39, 2001

CONTACT PROBLEMS WITH FRICTION, ADHESION AND

WEAR IN ORTHOPAEDIC BIOMECHANICS.

PART II – NUMERICAL IMPLEMENTATION AND

APPLICATION TO IMPLANTED KNEE JOINTS

Jerzy Rojek

Józef Joachim Telega

Stanisław Stupkiewicz

Institute of Fundamental Technological Research, Warsaw

e-mail: jrojek@ippt.gov.pl; jtelega@ippt.gov.pl; sstupkie@ippt.gov.pl

The present paper is the second part of the contribution by Rojek and
Telega (2001). An alternative adhesion law was used to the study of
bone-implant interface. Numerical schemewas developed and applied to
theknee jointafter arthroplasty. Influence ofweardebris on this interface
and currently used wearmodels were investigated.

Key words: unilateral contact, adhesion, friction, wear, knee joint after
arthroplasty, FEM

1. Introduction

In the first part of the paper, a general model of adhesion with friction
and generalized unilateral contact conditions was developed for two deforma-
ble bodies, within the range of small deformations. Such a model applies to
not necessarily linear elastic bodies. However, having in mind application to
orthopaedic biomechanics, only linear elastic solids have been considered.

In the second part of the paper we develop a numerical implementation
for unilateral contact problems with friction and adhesion, see Sections 2-4.
Particularly, in Section 2 a general form of finite element equations is given.
Section 3 is concerned with regularization of contact constraints. Numerical
contact algorithm is elaborated in Section 4. In Section 5 we first test the
developed numerical algorithm and next, perform the analysis of proximal
tibia after total knee replacement (TKR). The last section deals with wear in



680 J.Rojek et al.

joints after arthroplasty. Currently used phenomenological models of wear are
also discussed.

2. Finite element equations

Introducing space discretization of displacements

u(t,x)=N(x)a(t) (2.1)

to the set of differential equations from Part I (Rojek and Telega, 2001) and
applying the standard finite element procedurebased on theGalerkinmethod,
cf Zienkiewicz and Taylor (1998), we obtain the following set of semi-discrete
equations of motions

Mä=Fext+Fc−F int (2.2)

whichare thebasis of the explicit dynamic formulation. In the above equations
N is the matrix of interpolation functions, a and ä are the vectors of nodal
displacements and accelerations, respectively, M is the mass matrix, Fext,
F int and Fc are the vectors of external loads, internal and contact forces,
respectively. The vectors andmatrices in Eq. (2.2) are defined as follows

M=

∫

Ω

ρN⊤N dΩ

Fext =

∫

Ω

N
⊤
f dΩ+

∫

Γ1

N
⊤
g dΓ (2.3)

F int =

∫

Ω

B
⊤
σ dΩ

where σ is the Cauchy stress tensor, B is the strain-displacement operator
matrix, g is the vector of prescribed traction, f is the vector of given body
forces. Calculation of nodal contact forces Fc will be explained later on.
Equations (2.2) are integrated in time using an explicit scheme, in which

the displacements ak+1 at time tk+1 are obtained from the equations for the
known configuration at time tk

äk =M−1D

(
Fkext+F

k
c −F

k
int

)
MD = diagM

ȧ
k+1/2 = ȧk−1/2+ äk∆tk+1/2 ∆tk+1/2 =

1

2
(∆tk+∆tk+1)

a
k+1 =ak+ ȧk+1/2∆tk+1

(2.4)



Contact problems with friction... 681

Now the superscript k denotes that a variable is taken at instant tk. The
effectiveness of the explicit dynamic formulation is based on the use of a dia-
gonal mass matrix MD, so there is no need to solve a system of equations.
The main disadvantage is the limitation of the time step size ∆t due to the
conditional stability of the integration scheme.
The new contact formulation has been implemented within the dynamic

explicit formulation in thefinite element programSimpact (Rojek et al., 1998).
Because of its efficiency in the analysis of large-scale systems the dynamic
explicit approach has become very popular in many applications, cf Rojek et
al. (1998), Chenot and Bay (1998).
The explicit dynamic codes are best suited to carry out simulation of dy-

namic processes, nevertheless they can be applied to quasi-static simulation
as well. We can obtain approximate solutions of linear and non-linear quasi-
static problems introducing adequate damping which would damp out the
oscillations. Adding the damping term to Eq. (2.2) we have

Mä+Cȧ=Fext+F c−F int (2.5)

where C is the dampingmatrix and ȧ is the vector of nodal velocities.
Usually it is sufficient to damp out the lowest vibration modes. This can

be achieved by applying global viscous damping close to the critical value for
the lowest frequencies. The loading velocity should correspond to the dynamic
properties of the system as well, and it should be introduced for a time suf-
ficiently long in comparison to the period of vibrations in the lowest natural
mode. Using this methodology we have performed the analysis of total knee
replacement subjected to quasi-static loading, presented in Section 5.
Dynamic relaxation, amethodology elaborated to analyse quasi-static pro-

blems with dynamicmodel is presented byUnderwood (1983). For survey pa-
pers on finite element method in biomechanics the reader is referred to Mac-
kerle (1998) andPrendergast (1997), cf also Joshi et al. (2000), Rakotomanana
(2000).

3. Regularization of contact constraints

The complete set of contact conditions that must be taken into account
consists of the conditions for normal contact with adhesion, cf Rojek and
Telega (2001)

un ­ 0 −R
e
n+knunβ

2 ­ 0 (−Ren+knunβ
2)un =0 (3.1)



682 J.Rojek et al.

relations for tangential contact with adhesion and friction

R
e
T = kTuTβ

2

φ¬ 0 λ­ 0 φλ=0 (3.2)

u̇T +λR
i
T =0

where
φ= ‖RiT‖−µ|Rn−knunβ

2| (3.3)

as well as the evolution law for the adhesion intensity β given by relations

β̇=−
{
−
1

b
[w− (knu

2
n+kT‖uT‖

2)β]−
}1/p

if β∈ [0,1[

β̇¬−
{
−
1

b
[w− (knu

2
n+kT‖uT‖

2)]−
}1/p

if β=1

(3.4)

Regularization of contact constraints simplifies the variational formulation
of the contact problem.Herewe shall apply the penaltymethod to enforce the
constraints for normal and tangential constraints.
The regularization of the normal contact conditions (3.1) is accomplished

by introducing the penalty coefficient k−n . The impenetrability condition is
enforced by

R−n = k
−

nu
−

n (3.5)

where R−n and u
−
n stand for the negative parts of Rn and un, respectively

(zero irreversible normal traction is taken into account here). The relation
(3.5) applies to compression only, the relations for adhesive traction

Rn = knunβ
2 if un > 0 (3.6)

remainsunchanged.For consistency in notation, however,we shall rewrite it as

R+n = k
+
nu
+
nβ
2 (3.7)

where R+n and u
+
n stand for thenonnegativeparts of Rn and un, respectively.

The above modifications should be taken into account in Eq. (3.6), which
now takes the form

φ= ‖RiT‖−µ|R
−

n | (3.8)

The regularization of frictional constraints is carried out by introducing into
Eq. (3.2)3 the friction tangential penalty coefficient εT

u̇T +λR
i
T =

1

εT
Ṙ
i
T (3.9)



Contact problems with friction... 683

where Ṙ
i
T is the time derivative of R

i
T . For large displacement formula-

tion the Lie derivative should be used here, cf Laursen and Simo (1993). The
regularization becomes exact as εT →∞.

Penalization of normal and tangential contact constraints is equivalent to
specifying additional constitutive relations on the interface. Thus the elastic
behaviour has been introduced here for compression in normal contact, and
the frictional contact has been reformulated as a problem analogous to that
of elastoplasticity, the elastic component of tangential velocity being specified

by (1/εT)Ṙ
i
T .

Regularized contact conditions constitute the basis of the numerical algo-
rithm for the calculation of contact interaction forces.

4. Numerical contact algorithm

At every step of the time integration the contact algorithm performs two
basic tasks:

• contact search to establish the contact area,

• evaluation of contact forces.

Fig. 1. Discretized contact surfaces

The contact search in the implemented finite element algorithm is carried
out by checking the location of the points discretizing one of the two surfaces,
called ”slave”,with respect to theother one, called ”master” (seeFig.1). In the
standard contact algorithm (without adhesion), contact is established for the
slave nodes penetrating through the master surface (un < 0). In the present
algorithm the contact search has been modified to include the possibility of
adhesive contact (un > 0).



684 J.Rojek et al.

Once the contact is established between a certain slave node S and a
master segment definedwithnodes 1 through 4 (Fig.1) at time tk, the relative
displacements ukn and u

k
T , the tangential relative velocity u̇

k
T are obtained.

Given these variables and interface constitutive parameters, k−n , k
+
n , kT ,w,µ,

b and p, the contact interactions Rkn and R
k
T , and intensity of adhesion β

k

are calculated. Previous values of irreversible tangential force (RiT)
k−1 and

intensity of adhesion βk−1 are required to calculate new state at the contact
interface.
The contact force applied to the node S is obtained by integrating the

contact traction Rkn and R
k
T over a part of the boundary ΓS associated with

the node S (see Fig.1)

(Fkc)S =(R
k
nn+R

k
T)|ΓS| (4.1)

|ΓS| – surface measure of ΓS.
The reaction applied to the master surface is next replaced by equivalent

forces applied to the master segment nodes 1 through 4

(Fkc)i =−Ni(F
k
c)S i=1, . . . ,4 (4.2)

where Ni is the shape functioncorrespondingto the ithnodeat theprojection
of the slave node S on the master segment.
Thenumerical algorithmused to calculate the state of the contact interface

is summarized below as follows:

1. Evolution of intensity of adhesion

If βk−1 > 0 (if adhesion bonds exist), a trial value of β̇k is calculated
using relation (3.4)1

β̇ktrial =−
{
−
1

b

[
w−
(
k+n (u

(+)k
n )

2
+kT‖u

k
T‖
2
)
βk−1

]}1/p
(4.3)

Normal contact force contribution is taken into account in Eq. (4.3) for
tensile interface traction only. If β̇ktrial < 0,we have a decohesion process
and a new value of the intensity of adhesion is calculated using forward
Euler scheme

βk =βk−1+ β̇ktrial∆t
k (4.4)

otherwise we set
βk =βk−1 (4.5)



Contact problems with friction... 685

2. Normal contact force

If ukn ¬ 0 the compressive traction is calculated from Eq. (3.5)

Rkn = k
−

nu
k
n (4.6)

otherwise the tensile tractiondue to adhesion is calculated fromEq. (3.7)

Rkn = k
+
nu
k
n(β
k)2 if βk > 0 (4.7)

or Rkn is set to zero

Rkn =0 if β
k =0 (4.8)

3. Adhesive tangential contact force

Adhesive tangential contact traction is calculated from Eq. (3.2)1

(ReT)
k = kTu

k
T(β
k)2 if βk > 0 (4.9)

or (ReT)
k is set to zero, if adhesion is absent

(ReT)
k =0 if βk =0 (4.10)

4. Frictional tangential contact force

For tensile normal contact traction (Rn > 0) the friction force is set to
zero

(RiT)
k =0 (4.11)

If the normal contact traction is compressive, the friction contact trac-
tion is calculated using the radial return algorithm analogous to that
used in elastoplasticity. First we calculate a trial state employing Eqs.
(3.9)

(RiT)
k
trial =(R

i
T)
k−1+εTu̇

k
T∆t

k (4.12)

and check the slip condition (3.8)

φtrial = ‖(R
i
T)
k
trial‖−µ|R

k
n| (4.13)

If φtrial ¬ 0, we have the case of stick contact and the frictional traction
is assigned the trial value

(RiT)
k =(RiT)

k
trial (4.14)

otherwise (slip contact) we perform a return map

(RiT)
k =µ|Rkn|

(RiT)
k
trial

‖(RiT)
k
trial‖

(4.15)



686 J.Rojek et al.

Thus we have the state of the contact interface completely defined. The
above algorithm takes into account all possible loading, unloading and relo-
ading conditionswhichwill bedemonstrated in an example in thenext section.
The intensity of adhesion in the present formulations tends to zero (β → 0)
as the process of decohesion continues, but never achieves zero. For practical
handling of cases with adhesive bondsbrokenwehave introduced a small limit
value βmin, below which we assume that the adhesion disappears, and then
we set β=0.

5. Numerical examples

The aim of this section is twofold. First, we will test only the elaborated
numerical algorithm for contact problem with adhesion and friction. Second,
the discretized interfacemodelwill be applied to the analysis of proximal tibia
after TKR.

5.1. Study of a contact interface under different loading conditions

We analyze a contact interface in 2D with different properties and under
different loading conditions. Using the developed routines we calculate the
contact forces under prescribed relative displacements. In the first case we
analyze the normal behaviour of a contact interface with adhesion with the
following properties: kn = 20kPa/mm, w = 2kJ/mm

2, p = 1. The loading
was imposed by the constant velocity u̇n = 2mm/s, t ∈ [0,0.2]s. We have
analyzed the interface for different viscosity: b = 1 · 10−4, 1 · 10−2, 5 · 10−2

and 1 · 10−1kN·s/mm. We start from the total adhesion (β = 1) and zero
displacement (uT = un =0). Figure 2 depicts the normal force history obta-
ined for these conditions.We can see that until the elasticity energy is smaller
than the limit of adhesion energy, the force increases linearly. After passing
this limit adhesion starts to decrease, although for higher values of viscosity
the adhesion force can be higher than the limit resulting from the limit of ad-
hesion energy. The higher the value of viscosity b, the longer will be the delay
in the adhesion force decrease. Finally the interfacial force decreases nearly to
zero, what means that the adhesion bonds are nearly broken.

The same interface was loaded with prescribed velocity u̇n =2mm/s for
t∈ [0,0.1]s, then unloaded −u̇n =−2mm/s for t∈ [0.1,0.15]s, and loaded
−u̇n = 2mm/s for t ∈ [0.15,0.3] s. The problem was studied for two values
of viscosity: b = 1 · 10−4 and 1 · 10−1kN·s/mm. The results in the form



Contact problems with friction... 687

Fig. 2. Adhesion force for different viscosity

Fig. 3. Adhesion force for different viscosity under loading and unloading



688 J.Rojek et al.

Fig. 4. Tangential force for adhesion with friction

Fig. 5. Tangential force for adhesion with friction under loading and unloading



Contact problems with friction... 689

of relation adhesion force vs. relative displacement are shown in Fig.3. We
can see the unloading for the stiffness of the adhesion decreased.We can also
see the viscous effect in unloading for higher value of the viscosity. It is thus
shown that the adhesive bonds can keep breaking even during unloading if the
viscosity is present.

We analyze the tangential behaviour of a contact interface with adhe-
sion and friction The interface has the following properties: k−n = kT =
20kPa/mm, w = 2kJ/mm2, b = 1 · 10−1kN·s/mm, p = 1, µ = 0.3. The
loading was imposed by the constant velocity u̇T = 2mm/s, t ∈ [0,0.2]s.
Constant compression was applied by setting constant normal displacement
un = −0.5mm. Figure 4 presents the tangential force history obtained. We
can see that tangential force has the frictional component equal to the sliding
limit µRn. The tangential force after breaking adhesive bonds tends to the
sliding limit now.

The same interface was loaded with the imposed velocity according to the
following program: u̇n =2mm/s for t∈ [0,0.1]s, u̇n =0 for t∈ [0.1,0.11] s,
u̇n = −2mm/s for t ∈ [0.11,0.16] s, u̇n = 2mm/s for t ∈ [0.16,0.31] s.
Figure 5 represents the tangential force historyobtained.Wecan see thatwhen
the loading is constant, the adhesion keeps decreasing because of relaxation.
Under unloading conditions the adhesive and friction forces have opposite
orientation. Under reloading the sign of friction force is changed and again is
the same as that of the adhesion component. When adhesion disappears, the
tangential force is equal to the sliding limit.

The examples analyzed here demonstrate correctness of the formulation
andnumerical algorithm.The frictionandadhesion forces havebeencalculated
correctly underdifferent complex loading.This verification allows us to use the
numerical algorithm for real problems of contact with adhesion and friction.

5.2. Analysis of adhesion in the total knee replacement

Weapply the interfacemodelwith adhesion to ananalysis of proximal tibia
after TKR. Fixation of prosthesis is one of the most important mechanical
aspects of the performance of the bone-implant system. As we already know,
loosening can be caused by variety of mechanical and biological phenomena,
cf Sections 2, 3 in Rojek and Telega (2001) and Section 6 which follows. In
the present analysis we shall study loosening of a cemented stem initiated
by debonding at the stem-cement interface, and we shall focus on the role of
adhesion in the fixation of implants.

The knee experiences an oscillating varus-valgus loading.Off-centered load
can cause tensile stresses on the knee implant interfaces (Fig.6). Cement has



690 J.Rojek et al.

poor tensile properties, therefore debonding can occur under tension. One of
the design objective of a knee implant is avoidance of the interface tension.

Fig. 6. Interface tension and compression in the knee implant under off-centred
loading

Fig. 7. Finite elementmodel of tibia with implant

The proximal tibia with mounted tibial component of the condylar re-
placement has been discretized with tetrahedral finite elements (Fig.7). The
tibial component of the endoprosthesis has a plastic bearing surface placed on
a metal tray with central cemented post. Isotropic elastic properties for the
materials in the model have been assumed as follows:

• cortical bone: Young’s modulus 1.7 ·1010Pa, Poisson’s ratio 0.32,

• cement: Young’s modulus 2.27 ·109Pa, Poisson’s ratio 0.335,



Contact problems with friction... 691

• steel: Young’s modulus 2.14 ·1011Pa, Poisson’s ratio 0.30,

• polyethylene: Young’s modulus 1 ·109Pa, Poisson’s ratio 0.30.

Contactwithadhesionand frictionhasbeen taken intoaccount at the steel-
bone and steel-cement interfaces, whilematerials at other interfaces have been
considered to be fully bonded.The following properties have been assumed for
all the interfaces with adhesion:

• interface normal stiffness kn = k
−
n =1 ·10

12Pa/m,

• interface tangential stiffness kT =1 ·10
11Pa/m,

• energy of adhesion w=0.36J/m2,

• viscosity parameter b=0,

• Coulomb friction coefficient µ=0.1.

The above interface properties yield the limit tensile stress 0.6MPa at the
separation 0.6 ·10−3mm, and the limit shear interface strength 0.19MPa at
the tangential relative displacement 1.9 · 10−3mm (friction, however, raises
total tangential resistance). These properties correspond to a relatively weak
interface strength, weaker than that of real bonding. Such an assumption,
however, was necessary to allow the interface to debond.

Fig. 8. Finite element model with loading condition

The tibia has been fixed distally at the base of the bone and the loading
has been applied non-symmetrically on the bearing surface (Fig.8). The load



692 J.Rojek et al.

condition approximated a tibio-femoral contact forces in one leg standing po-
sition. The compressive force of 400% BW (body weight) has been assumed.
For BW = 700N we have Rmax = 2600N. The force has been distributed
over the supposed contact area between tibial and femoral comoponents.

Quasi-static loading has been analyzed using dynamic relaxation method.
Loading has been growing linearly from 0 to the final value of 400%BW.The
quasi-static behaviour has been obtained by introducing adequate (close to
critical) damping. Distribution of stresses in the tibia with implant under the
loading of 320% BW is shown in Fig.9b. Figure9a shows the displacement
vectors for the same level of loading; it can be clearly seen that the side of the
tray opposite to the loading is raised. This causes tensile contact tractions.

Evolution of interface normal stresses on the tray and stem surfaces is
shown in Fig.10. Evolution of tangential contact stresses is shown in Fig.11.
Thedecohesion process can beanalyzed studying the evolution of the intensity
of adhesion β, cf Fig.12.

The results inFigs.10-12 arepresented for thebondedarea only, the reduc-
tion of this area corresponds to the decohesion process. Initially we have bon-
ding at the whole surface. Until certain level of loading the bonding strength
is not reached and the bonding is not broken. As we can see in Fig.12, the
bonded area up to 320% BW is practically not changed. For higher loads de-
cohesion takes place.We can see that once decohesion has started, increasing
loading leads quickly to to complete loosening of the prothesis.

6. Wear in implanted joints, wear models

Looking back, one should first mention ferrography, which was invented
in the early 1970’s and has gained use in the analysis of wear in machines,
cf Evans et al. (1980). The ferrographic technique is based on the magnetic
precipitation of particles from suspensions. In the papers by Evans (1983),
Evans et al. (1980, 1981,1982), Mears et al. (1978), possibilities of employing
the ferrography for the analysis of wear particles present in the synovial fluid
aspirated from surgical joint replacements and human arthritic joints were
examined.

Passing to recent results, a collection of papers contained inDowson (1998)
presents a good survey of experimental results on the subject of friction and
wear of joint replacements and implant materials, cf also Saikko and Ahlroos
(2000). Sophisticated single station machine simulators capable of measuring





 



Contact problems with friction... 695

both the wear and friction were also described, cf also Saikko (1993), Walker
et al. (1997).

6.1. Basic wear effects

In the context of joint prostheses, there are two main groups of wear ef-
fects. Thefirst group of problems is related towear phenomena in the artificial
joint and canbe clasified as amostly tribological problemwhen the contacting
bodies are traditionalmaterials, usuallymetal or ceramic andultra-highmole-
cular weight polyethylene (UHMWPE). However, even the implanted joint is
lubricatedby the synovial fluid, thoughnot so effectively as natural healthy jo-
ints, cfWang et al. (1998). The second group of problems is related to contact
phenomena at bone-implant interfaces, cf the previous sections. In this case,
one of the contacting bodies is a living tissue, thus several non-mechanical
issues occur. The two tribological systems are strongly coupled, namely the
UHMWPE wear debris is commonly agreed to be a major cause of long term
failure of total joint replacement through aseptic loosening, cf Laffargue et al.
(1998), Revell et al. (1997). Aseptic loosening is a process of prosthesis loose-
ning due to bone loss in the absence of infection or mechanical failure. The
bone loss process is induced by the inflamatory response to UHMWPE wear
particles. A coupling effect in an opposite direction is also observed, namely
the bone cement andboneparticles originating from thebone-prosthesis inter-
face may penetrate the joint and strongly accelerate the wear process due to
the third body abrasion. Extensive experimental studies on friction and wear
of ceramics were carried out by Sasaki (1992). However, a situation typical for
ceramic bone implants was not investigated.

For these reasons, the study of the two tribological systems and their re-
lative coupling is an important area of biomechanics. Many efforts have been
concentrated on experimental studies of wear performance of UHMWPE by
means of traditional pin-on-disk devices (Klapperich et al., 1999; Saikko and
Ahlroos, 2000), special joint simulators (Saikko, 1993) and analysis of explan-
ted joints after in vivo operation (Atkinson et al., 1985). Dowson (1998) and
Wanget al (1998) provide extensive reviewsof the recent advances in the tribo-
logy of artificial joints. For instance, topics such as themajorwearmechanisms
and lubrication regimes, effects of surface roughness and joint kinematics on
wear mechanisms are discussed in the paper byWang et al. (1998).

Wear particles of different materials are found in the cells next to the
loosened artificial joint prostheses, but also in distant locations, such as liver,
spleen and kidney (Revell et al., 1997). The transport of particles by the
host defense mechanisms is affected by the particle material, size and shape,



696 J.Rojek et al.

therefore it is not only the wear volume that is important, but also the size
and shape of wear particles generated at contact interfaces, cf Dowson (1998).
Clearly, with decreasing size the number of particles per unit wear volume
increases. In fact, the particle sizemay vary from0.2µmto 100µm, cf Ingham
andFischer (1999), whichmay lead to a several orders ofmagnitude difference
in the number of particles for a fixed wear volume.

6.2. Models of plasticity-determined wear mechanisms

There are threemain basicwearmechanisms occuring at artificial joint in-
terfaces, namely adhesive wear, abrasive wear and surface fatigue wear (Wang
et al., 1998). According to Lim and Ashby (1987), the adhesive and abrasive
(and also to some extent fatigue) wear mechanisms are plasticity-determined.
Adhesive wear is explained by welding at asperity contacts followed by remo-
val of shearing (due to sliding) of softer asperity fragments which can further
create wear particles, cf Archard (1953). Abrasive wear is caused by hard
asperities or particles (two- or three-body abrasive wear) ploughing the softer
surface. This is accompanied by severe plastic deformation of the subsurface
layer. In some cases (e.g. hard and sharp particles) ploughing can cause chip
cutting, thus leading to direct formation of wear particles butmore oftenwear
results from repeated plastic deformation andmicro-crack propagation.

Wear due to fatigue is caused by crack propagation caused by repeated
cyclic loading of the surface. This can occur on the macro-level, as in rolling
contact, or on micro-level due to multiple asperity contacts during sliding.
The latter case is often recalled when discussing adhesive or abrasive wear
mechanisms since wear rates resulting from assumption that every asperity
contact produces wear particle are several orders of magnitude higher than
those observed experimentally.

Despite the fundamental differences betwen the differentwearmechanisms
(and respective friction mechanisms) they have the plastic deformation origin
in common. Since the microscopic plastic deformation in the subsurface layer
is a manifestation of the macroscopic phenomenon of friction, the correspon-
dingdissipation rates are directly related.Therefore, theplasticity-determined
wear mechanisms are well described by the Archard wear law, which predicts
the wear rate to be proportional to the frictional dissipation rate. Of course,
depending on the contact conditions (surface roughness, range of contact pres-
sures, lubrication regime, etc.) different friction andwearmechanismsmay be
active for the same pair of contacting materials and each mechanism may be
characterized by a different proportionality factor of the wear law.

Below, the classical model of adhesive wear of Archard (1953) is outlined



Contact problems with friction... 697

and the wear law is derived. Next, the delamination theory of wear is briefly
presented, as it seemsmuch relevant for description of wear phenomena at the
contact interfaces in joint prostheses.

6.2.1. Archard law of adhesive wear

A simple model of adhesive wear has been proposed by Archard (1953).
Assuming that wear particles created at single asperity contacts have similar
shapes and their size is characterized by the diameter a of single real contact
areas, the volume of thewear particle is proportional to a3. Further, assuming
that the sliding distance over which the wear particle is formed equals the
contact area diameter a (after sliding this distance, the next asperity takes
the load) the volumetric wear rate W , i.e. the volume worn per unit sliding
length, is proportional to the real contact area A

W ∝A

For moderate normal pressures a plastic contact of asperities occurs and the
real contact area is proportional to the ratio of the normal load FN and the
hardness H of the softer body

A=
FN
H

(6.1)

So the wear rate per unit sliding length can be expressed as

W =
KFN
H

(6.2)

where K is the wear coefficient expressing the probability of a single asperity
adhesive contact to form a wear particle. According to Archard’s theory, the
wear rate is proportional to the ratio of the normal force (or pressure) to the
hardness of the softer material. The Archard’s wear law is simple and turned
to be valid in a wide range of sliding conditions, not only those characterized
by asperity adhesion (see the following section).
In continuummechanics it is more convenient to use rates with respect to

time rather than to the sliding length and to use local variables (as the contact
pressure pN) rather than the global ones (as the normal force FN). So the
Archard’s wear law (6.2) can be expressed in a local form

ẇ=
kpNv

H
(6.3)

where ẇ denotes the depth of wear in unit time and v is the relative sliding
velocity. Remind now that for the range of pressures for which the real area of



698 J.Rojek et al.

contact is proportional to normal pressure, cfEq. (6.1), similar proportionality
exists between the friction stress pT and normal pressure pN. The wear law
(6.3) can be further rewritten to yield

ẇ=
k̃Ḋ

H
(6.4)

where Ḋ= pTv is the frictional dissipation rate and k̃=µk is the modified
wear coefficient (µ denotes the friction coefficient).

6.2.2. Delamination theory of wear

It is easy to imagine material removal due to asperity shear-off or micro-
cutting by a ploughing asperity. However, in many contact conditions those
twomechanisms donot occur but the surfaces are still worn.Thedelamination
theory of wear, originally developed for metals, explains another mechanism
of material removal from contact surface during sliding, cf Suh (1973, 1977),
by assuming void and micro-crack nucleation at some characteristic distance
below the surface, as a result of plastic contact of single asperities, followed by
crack propagation parallel to the surface and delamination of flake-like wear
particles.
There are several facts justifying thismechanism. Themechanism of crack

nucleation and propagation at some characteristic depth below the surface is
explained by high compressive stresses which occur just below the contact re-
gions. This depth is usually of the order of asperity size. Plastic deformation
caused by a single pass of an asperity accumulates due to repeated action
during sliding. Thus the voids are created around inclusions or due to dislo-
cation pile-ups. Fracture mechanics analysis confirmed preferred direction of
crack propagation in the plane parallel to the surface. Finally, the flake-like
wear particles (with aspect ratio greatly different from unity) have been obse-
rved experimentally. An overview of the delamination theory of wear can be
found in Suh (1977).
In order to show that the delamination theory of wear leads also to the

wear law of the Archard’s type we shall adopt the model developed by Lim
and Ashby (1987). This model is quite simple and reflects the main ideas of
Suh’s model although it differs in some details.
Voids nucleate at inclusions at a characteristic depth X0 due to repeated

frictional traction on the sliding surface. The area fraction of voids fA on a
plane parallel to the sliding surface (at the depth X0) is

fA =(2riNv)(2ri)(2riγ)≈ 2fvγ (6.5)



Contact problems with friction... 699

where ri is the inclusion radius, Nv is the number of inclusions per unit
volume, γ is the cumulative plastic shear strain and fv is the volume fraction
of the inclusions. A critical value of area fraction fA = f

∗

A is assumed to exist
at which fracture causes a flake-like wear particle to break off. The critical
shear strain associated with this fracture γ∗ is then

γ∗ =
f∗A
2fv

(6.6)

and the number of passes n∗ needed to build up γ∗ follows from the plastic
shear strain γ0 accumulated per pass

γ∗ =n∗γ0 (6.7)

The thickness of the wear flake that is created after n∗ passes is X0, so
the wear rate is

W =
X0An
n∗d

(6.8)

where An is the nominal area of contact and d is the asperity spacing. Using
Eqs (6.6) and (6.7), the wear rate can be expressed as

W =
2γ0fvX0
f∗Ara

FN
H

(6.9)

where the distance between asperity contacts d has been replaced by the
equation

ra
d
=
Ar
An
=
FN
H

(6.10)

Here ra is the radius of asperity contact, Ar and An denote real and nominal
contact areas, respectively, FN is the normal force and H is the hardness.No-
te that Eq. (6.10) is valid only for long wedge-like asperities (two-dimensional
asperity model), cf Lim and Ashby (1987).

Finally, assuming that X0 = ra, what expresses the fact that the characte-
ristic depth of crack nucleation and propagation is of the order of the asperity
radius, the wear rate can be rewritten in the final form of Archard’s law

W =
2γ0fv
f∗A

FN
H

(6.11)

where the wear coefficient has been predicted as K =2γ0fv/f
∗

A.



700 J.Rojek et al.

Remark 6.1. Strömberg et al. (1996a,b) developed a continuum thermody-
namicmodel for interfacial phenomena includingunilateral contact, fric-
tion and wear, cf also Mróz and Stupkiewicz (1994), Johansson (1992),
Strömberg (1997). Themodel is confined to small displacements, imply-
ing small slip. A numerical model allowing for large slip was elaborated
by Agelet de Saracibar and Chiumenti (1999), cf also Agelet de Saraci-
bar (1997), Heegaard (1993), Heegaard andCurnier (1993), Laursen and
Simo (1993). Only theArchard’s law for adhesive wearwas incorporated
into themodel developed byAgelet de Saracibar andChiumenti (1999).
A restricted class of initial-boundary value problems with friction and
wear was studied by usingmethods of variational inequalities in Garčıa
et al. (2000), Kuttler and Shillor (1998), Shiller and Sofonea (1999).We
observe that in none of the above approaches, the adhesionwas included.

Remark 6.2. The study of contact problems in orthopaedic biomechanics re-
quires models incorporating friction, adhesion and wear. However, wear
debris are produced between different components of a prosthesis and
different from the bone-implant interface; for instance between the tibial
and femoral components of condylar TKR. In other words, one has to
take into account the transport of wear debris to the interface where ad-
hesion takes place. On such an interface a possible transport equation is

∂C(x,t)

∂t
=∇·

[
a(x)∇C(x,t)−C(x,t)v(x,t)

]
+R
(
C(x,t)

)
(6.12)

Here C(x,t) denotes the concentration of wear debris on the interface
bone-implant, a is the matrix of diffusion coefficients, v is the (small)
velocity of wear debris, R denotes the reaction term and t denotes
time. Equation (6.12) has to be supplemented by the boundary and
initial conditions.

7. Final remarks

Longevity of prosthesis depends on the quality of bone-implant interface.
Though a lot of studies have been concerned with this subject, yet no general
biomechanical model is available. Such a model should incorporate adhesion,
friction and influence of wear debris. Then the density of adhesion would
depend on the debris concentration C(x,t), appearing in Eq. (6.12).We hope
to cope with this problem in the future.



Contact problems with friction... 701

Bergmann and his coworkers developed experimental procedure for the
determination of forces and friction-induced temperature in implanted hip
joints, cf Bergmann et al. (1993, 1995, 1997, 1999), Graichen et al. (1999) and
the references therein. Similar procedures has not yet been applied to knee
joint after arthroplasty (G.Bergmann, private communication, 2000).

The review paper by Feeny et al. (1998), containing a comprehensive tho-
ughnot exhaustive list of references, provides an interestinghistorical overview
of structural and mechanical systems with friction, beginning from the anti-
quity. However, striving for understanding the friction phenomenon in human
and animal joints was overlooked.

Recently, a new idea of modelling friction was elaborated byTworzydlo et
al. (1998). These authors proposed an approach for constructing constitutive
models of frictional interfaces,whichprovides a linkbetweenmicroscale pheno-
mena andmacroscale phenomenological models at the interface. Tworzydlo et
al. (1998) claim that theirmicro-macro approach is based on statistical homo-
genisation. Unfortunately, this point remains unclear since no homogenisation
problemwas formulated.

Hip and knee implants are partially in contact with cancellous bone,which
is a cellular solid, cf Gibson andAshby (1988), Tokarzewski et al. (1999). For-
tes et al. (1999) investigated the static contact of a cellular solid with open
cells with a compact counter-surface. Such a model seems to be applicable
to cancellous bone of rather low density in contact with a prosthesis, though
obviously then the contact is much more complicated. Namely, according to
Maher and McCormack (1999), the mechanical interlock created at the ce-
ment/cancellous bone interface in cemented femoral reconstruction is crucial
to the longevity of the replacement.This interlock is created throughfingers or
pedicles of bone cement which protrude through cancellous trabecular spaces.

Another interesting idea is due to Zhang and Tanaka (1997), cf also the
references therein. These authors studied themechanisms of friction andwear
on the atomic scale through the investigation of a typical diamond-copper
sliding system with the aid of molecular dynamics analysis. However, one
still lacks a theory bridging the gap between atomic analysis and that using
continuumanalysis. Sliding parts in a prosthesis usually consist of a pairmade
of UHMWPE-metal or UHMWPE-ceramics. Then it seemsmore appropriate
to characterise the UHMWPE by its molecules.

Acknowledgment

Theauthorswere supportedbytheStateCommittee forScientificResearch(KBN,

Poland) through the grant No. 8T11F01718. Thanks are given to Prof. Oñate from



702 J.Rojek et al.

the International Centre of Numerical Methods in Engineering in Barcelona for pro-

viding the numerical code Simpact.

References

1. Agelet de Saracibar C., 1997, ANew Frictional Time Integration for Lar-
ge Slip Multi-Body Frictional Contact Problems, Comput. Meth. Appl. Mech.
Eng., 142, 303-334

2. AgeletdeSaracibarC.,ChiumentiM., 1999,On theNumericalModeling
of FrictionalWearPhenomena,Comput.Meth. Appl.Mech. Eng., 177, 401-426

3. Archard J.F., 1953, Contact and Rubbing of Flat Surfaces, J. Appl. Phys.,
24, 981-988

4. Atkinson J.R., Dowson D., Isaac J.H., Wroblewski B.M., 1985, La-
boratory Wear Tests and Clinical Observations of the Penetration of Femoral
Heads into Acetabular Cups in Total Replacement Hip Joints III: TheMeasu-
rement of Internal Volume Changes in Explanted Charnley Sockets after 2-16
Years in Vivo and the Determination ofWear Factors,Wear, 104, 225-244

5. Bergmann G., Graichen F., Rohlmann A., 1993, Hip Joint Loading Du-
ring Walking and Running, Measured in Two Patients, J. Biomechanics, 26,
969-990

6. Bergmann G., Graichen F., Rohlmann A., 1995, Is Staircase Walking a
Risk for the Fixation of Hip Implants? J. Biomechanics, 28, 535-553

7. BergmannG.,GraichenF., RohlmannA., 1999,Friction InducedTempe-
rature Increase ofHip Implants, In: L. Sedel andG.Willmann,Edit.,Reliability
and Long-Term Results of Ceramics in Orthopaedics, 91-95, Georg Thieme-
Verlag, Stuttgart-NewYork

8. Bergmann G., Graichen F., Rohlmann A., Linke H., 1997, Hip Joint
Forces During Load Carrying,Clinical Orthop. Related Res., 335, 190-201

9. ChenotJ.-L,BayF., 1998,AnOverviewofNumericalModellingTechniques,
J. Mater. Proc. Technology, 80-81, 8-15

10. Dowson D., Edit., 1998, Advances in Medical Tribology, Mechanical Engine-
ering Publications, London

11. Evans C.H., 1983, Application of Ferrography to the Study of Wear and Ar-
thritis in Human Joints,Wear, 90, 281-292

12. EvansC.H., BowenE.R., Bowen J., Tew W.P.,Westcott V.C., 1980,
Synovial Fluid Analysis by Ferrography, J. Biochem. Biophys. Methods, 2,
11-18



Contact problems with friction... 703

13. Evans C.H., Mears D.C., McKnight J.L., 1981, A Preliminary Ferrogra-
phic Surveyof theWearParticles inHumanSynovialFluid,Arthritis andRheu-
matism, 24, 912-918

14. EvansC.H.,MearsD., StanitskiC.L., 1982,FerrographicAnalysis ofWear
in Human Joints: Evaluation by Comparison with Arthroscopic Examination
of Symptomatic Knee, J. Bone Joint Surg., 64-B, 572-578

15. Feeny B., GuranA., Hinrichs H., PoppK., 1998,AHistorical Review on
Dry Friction and Stick-slip Phenomena,Appl. Mech. Reviews, 51, 321-341

16. Fortes M.A., Colaço R., Vaz M. F., 1999, The Contact Mechanics of
Cellular Solids,Wear, 230, 1-10

17. Garćıa J.R.F., Han W., Shillor M., Sofonea M., 2000, Recent Advan-
ces in Variational and Numerical Analysis of Quasistatic Contact Problems,
Preprint, Departamento de Matematica Aplicada, Universidad de Santiago de
Compostela

18. Gibson L., Ashby M., 1988, Cellular Solids: Structure and Properties, Per-
gamon, Oxford

19. Graichen F., Bergmann G., Rohlmann A., 1999, Hip Endoprosthesis for
in Vivo Measurement of Joint Force and Temperature, J. Biomechanics, 32,
1113-1117

20. Heegaard J.-H., 1993, Large Slip Contact in Biomechanics: Kinematics and
Stress Analysis of the Patello-Femoral Joint, ThèseNo 1113, Ecole Polytechni-
que Fédérale de Lausanne

21. Heegaard J.-H., Curnier A., 1993, AnAugmented LagrangianMethod for
Discrete Large-Slip Contact Problems, Int. J. Num. Meth. Eng., 36, 569-593

22. InghamE., Fischer J., 1999,Wear ofUltraHighMolecularWeightPolyethy-
lene in Total Hip Joints,Biological and Mechanical Mechanisms, (Preprint)

23. Johansson L., 1992, Elastic and Thermoelastic Problems with Friction and
Wear, Linköping Studies in Science and Technology, Dissertations, No. 266

24. Joshi M.G., Santare M.H., Advani S.G., 2000, Survey of Stress Analyses
of the Femoral Hip Prosthesis,Appl. Mech. Reviews, 53, 1-18

25. KlapperichC.,KomvopoulosK.,PruittL., 1999,TribologicalProperties
and Microstructure Evolution of Ultra-High Molecular Weight Polyethylene,
ASME J. Tribol., 121, 394-402

26. Kuttler K.L., Shillor M., 1998, Dynamic Contact Problem with Normal
Compliance Wear and Discontinuous Friction Coefficient, Preprint, Dept. of
Mathematical Sciences, Michigan Technological University

27. Laffargue P., Breme H.J., Helsen J.A., Hildebrand H.F., 1998, Re-
trievalAnalysis, In: J.A.Helsen andH.J. Breme,Edit.,Metals as Biomaterials,
467-501, JohnWiley & Sons, Chichester



704 J.Rojek et al.

28. Laursen T.A., Simo J.C., 1993, A Continuum-Based Finite Element For-
mulation for the Implicit Solution of Multibody, Large Deformation Frictional
Contact Problems, Int. J. Num. Meth. Eng., 36, 3451-3485

29. Lim S.C., Ashby M.F., 1987, Wear-Mechanism Maps, Acta Metall., 35, 1,
1-24

30. Mackerle J., 1998, A Finite Element Bibliography for Biomechanics, Appl.
Mech. Reviews, 51, 583-634

31. Maher S.A., McCormack B.A.O., 1999, Quantification of Interdigitation
atBoneCement/CancellousBone Interfaces inCementedFemoralReconstruc-
tions,Proc. Inst. Mech. Engrs., 213, Part H, 347-354

32. Mears D.C., Hanley E.N., Rutkowski R., 1978, Ferrography: Its Appli-
cation to the Study of Human JointWear,Wear, 50, 115-125

33. Mróz Z., Stupkiewicz S., 1994, An Anisotropic Friction and Wear Model,
Int. J. Solids and Structures, 31, 1113-1131

34. Prendergast P.J., 1997, Finite Element Models in Tissue Mechanics and
Orthopaedic Implant Design,Clinical Biomech., 12, 346-366

35. Rakotomanana L., 2000, The Long-Term Behaviour of Human Joints with
Orthopaedic Prosthesis: Finite Element Models, In: European Congress on
Comput. Meth. in Appl. Sci. Eng., ECCOMAS 2000, Barcelona, 11-14 Sep-
tember, 2000

36. Revell P.A., Al-Saffar N., Kobayashi A., 1997, Biological Reaction to
Debris inRelation to JointProstheses,Proc. Instn.Mech. Engrs., PartH,211,
187-197

37. Rockafeller R., Wets R.-B., 1998,Variational Analysis, Springer, Berlin

38. Rojek J., Oñate E., Postek E., 1998, Application of Explicit FE code
to Simulation of Sheet and Bulk Metal Forming Processes, J. Mater. Proc.
Technology, 80-81, 620-627

39. Rojek J., Telega J.J., 2001,ContactProblemsWithFriction,Adhesion and
Wear inOrthopaedicBiomechanics.Numerical implementationandApplication
to ImplantedKneeJoints.Part I:GeneralDevelopments,J.Theor.Appl.Mech.,
39, 3

40. Saikko V.O., 1993,Wear of Polyethylene Acetabular Cups Against Alumina
Femoral Heads,Acta Orthop. Scand., 64, 507-512

41. Saikko V.O., Ahlroos T., 1999, Type of Motion and Lubricant in Wear
Simulation of Polyethylene Acetabular Cup, Proc. Inst. Mech. Engrs., 213,
Part H, 301-310



Contact problems with friction... 705

42. Saikko V.O., Ahlroos T., 2000, Wear Simulation of UHMWPE for Total
Hip Replacement with a Multidirectional Motion Pin-on-Disk Device: Effects
ofCounterfaceMaterial,ContactArea, andLubricant,J. Biomed.Mater. Res.,
49, 147-154

43. Sasaki S., 1992, Effects of Environment on Friction and Wear of Ceramics,
Bull. Mech. Eng. Laboratories, 58

44. Shillor M., Sofonea M., 1999,A Quasistatic Viscoelastic Contact Problem
with Friction, Preprint, Dept. ofMathematical Sciences, OaklandUniversity

45. Strömberg N., 1997, Thermomechanical Modelling of Tribological Systems,
Linköping Studies in Science and Technology, Dissertations, No. 497

46. Strömberg N., Johansson L., KlarbringA., 1996a,Derivation andAna-
lysis of a Generalized Standard Model for Contact, Friction and Wear, Int. J.
Solids and Structures, 33, 1817-1836

47. Strömberg N., Johansson L., Klarbring A., 1996b, AGeneralized Stan-
dardModel forContact,FrictionandWear, In:M.Raous,M.Jean,andJ.J.Mo-
reau, Edit.,Contact Mechanics, 327-334, PlenumPress, NewYork

48. Suh N.P., 1973, The Delamination Theory ofWear,Wear, 25, 111-124

49. SuhN.P., 1977,AnOverview of theDelaminationTheory ofWear,Wear, 44,
1-16

50. Tokarzewski S., Telega J.J., Gałka A., 1999, A Contribution to Evalu-
ation of EffectiveModuli of Trabecular Bonewith Rod-LikeMicrostructure, J.
Theor. Appl. Mech., 37, 707-728

51. Tworzydlo W.W., Cecot W., Oden J.T., Yew C.H., 1998, Computa-
tional Micro- and Macroscopic Models of Contact and Friction: Formulation,
Approach and Applications,Wear, 220, 113-140

52. Underwood P.G., 1983, Dynamic Relaxation, In: T. Belytschko and T.J.R.
Hughes, Eds, Computational Methods for Transient Dynamic Analysis, Am-
sterdam, North-Holland

53. Walker P.S., Blunn G.W., Broome D.R., Perry J., Watkins A., Sa-
thasivam S., Dewar M.E., Paul J.P., 1997, A Knee Simulating Machine
for PerformanceEvaluation of Total KneeReplacements, J. Biomechanics, 30,
83-89

54. Wang A., Essner A., Polineni V.K., Stark C., Dumbleton J.H., 1998,
Lubrication and Wear of Ultra-High Molecular Weight Polyethylene in Total
Joint Replacements,Tribol. Int., 31, 17-33

55. Zhang L., Tanaka H., 1997, Towards a Deeper Understanding of Wear and
Frictionon theAtomicScale – aMolecularDynamicAnalysis,Wear,211, 44-53

56. Zienkiewicz O.C., Taylor R.L., 1989, The Finite Element Method,
McGraw-Hill, London



706 J.Rojek et al.

Zagadnienia kontaktowe z tarciem, adhezją i zużyciem w biomechanice

ortopedycznej. Część II – Implementacja numeryczna i zastosowanie do

stawów kolanowych po implementacji

Streszczenie

Niniejsza praca stanowi część drugą pracy, której autorami są Rojek i Telega
(2001). Do analizy interfazy kość-implant zastosowano alternatywny opis adhezji.
Opracowano algorytm numeryczny, który zastosowano do analizy stawu kolanowego
po endoprotezoplastyce. Zbadanowpływproduktów zużycia na analizowaną interfazę
oraz przedstawiono aktualnie stosowanemodele zużycia.

Manuscript received February 2, 2001; accepted for print May 23, 2001