Mathematical models of friction for contact problems in elasticity

Wear, 113 (1986)

61

61 - 77

MATHEMATICAL MODELS OF FRICTION PROBLEMS IN ELASTICITY*

FOR CONTACT

J. J. KALKER Department Netherlands)

of Mathematics

and Znformatics,

Delft

University

of Technology,

Delft (The

Summary

In this paper a number of mathematical models of friction are reviewed for contact problems between three-dimensional elastic bodies. These models are (1) simplified theory (Kalker), (2) the theory of Duvaut and Lions and (3) plasticity type theories.

1. Introduction

The purpose of this paper is to review a number of models of frictional contact between three-dimensional elastic bodies with special attention to numerical aspects. These models are (1) simplified theory (Kalker 1973, 1982) [ 1, 21, (2) the theory of Duvaut and Lions and its generalizations (1972 and later) [ 3,4] and (3) plasticity-type theories (initiated by Frederiksson in 1976) [ 51. To that end we consider two elastic bodies V, (a = 1, 2) which are rigidly mounted on rigid axles (see Fig. 1). They are pressed together so that a contact area forms between them. Two Cartesian coordinate systems 3t,i are introduced (i = 1, 2, 3 like all small Latin indices with the exception of a and n) which are fixed to the axles of the bodies. At the contact area local coordinates are defined, i.e. the single normal coordinate x, (normal (n)) along the outer normal on body 1 and two orthogonal tangential coordinates x, (T = 1, 2 as all Greek indices; Greek indices indicate tangential components) in the local tangent plane of the contact area. We adopt the summation convention for Greek and Latin indices. The index a is omitted when it is clear which body is intended. The bodies undergo an elastic deformation. Connected with such a deformation is an (undeformed) reference state which is obtained by removing all body forces and surface loads from the bodies and enforcing static equilibrium while the axles are kept in the same position as in the deformed state. The particles of the bodies are labelled by their position *Paper presented at the Seminar on Friction and Contact Noise, Delft Technology, Delft, The Netherlands, June 20 - 21,1985. 0043-1648,/86/$3.50

0 Elsevier Sequoia/Printed

University of

in The Netherlands

62

Fig. 1. Two bodies pressed together. Fig. 2. The potential contact area 3 V,.

in the undeformed state; their position in the deformed state is X,i. X aoi The displacement is defined by = xk

%itxaof~

-x6zoi

(1)

The displacement components u,~ are assumed to be small in the sense that

(2) where bi is any vector and D is the characteristic diameter of body Q at xooi, and the displacement gradients U,i,j are likewise small ahi Ual,j = axooj lu,iJ

(3)

Q 1

IbjjI

I=

(bfjbij)l'*

where b,j is any tensor element. The (linearized) strains are defined by def

Q

= $(Uf,

j

+

Lkj,

eij = ejf leijl

j)

(4)

4 1

Stresses urj and a body force per unit volume fr act upon the particles of the bodies. When a particle of the body is not in static equilibrium it undergoes an acceleration Ci where (5) is the material differentiation (where t is the time). The equations of equilibrium acting in V, are crij,j +

fi -j+&=()

(61

where p is the mass density, iii = 0 at static equilibrium and ii, + 0 in a dynamic system.

63

In elasticity, CT ii =

the stress may be derived from a potential

aw(ekk)

(7)

aeij

where eii = eji etc. (see eqn. (4)) and W is the strain energy. As a consequence of the symmetry of the strains the stresses are likewise symmetric uij =

Uji

(8)

According

to Hooke’s law for anisotropic

W = +Eijhkeijehk

__f

bodies W is quadratic

Dij = Eijhkehk (9)

Eijhk

= Ehkij

= Ejtkk

where EUj,k etc. are elastic Ej,hkehk may be inverted. We obtain

constants.

The

stress-strain

relations

eij =S ijhk ghk W(eij)

= iEijhkeijehk

Uij =

(10) = iSijhk(JIj(Jhj

= u(Otj)

where S is the tangential component of the velocity of body 1 with respect to body 2 at the point xi E aV, and time t. Equations (7), (9) and (10) are called the constitutive relations. The boundary V, of each body consists of three parts, i.e. ax,, aI&, and aV,. In aI$, the displacement ui is prescribed as iii ui = iii at

av,,

(11)

In particular, the surface of the axle of body u belongs to aI&, and a, = 0. In a V,, the surface traction is prescribed as Pr pj=pi

at av aP

The surface traction Pi = Uijnj

(12)

pr is connected

to the surface stress by the relation (13)

where nj is the outer normal on aV. dV, is called the potential contact area. It consists of the contact area and its immediate neighbourhood. We assume that the displacements, displacement gradients, surface loads and normals on ax at two opposing points x,~ of the bodies at aV, (see Fig. 2) may be identified with the same quantities at the point xi of aV,. In particular, owing to Newton’s third law, we have Pi = Pli

= -PZi

(14)

while we also assume

ni = TZli= - nzi

(15)

64

At the potential contact we assume contact conditions. These are separated into normal conditions which are connected with the contact formation and tangential or frictional conditions. We consider the normal conditions as follows. (1) Let d be the distance between the opposing surfaces in the deformed state. We must have that d(Xi) > 0

(16)

i.e. no penetration d=O contact. d may be expressed in terms of the distance h in the undeformed and the displacement difference Aui = Uii - u2i (see Fig. 3)

state

(17) Aui = Uii -Uzi (2) We also assume that on 8 V, , outside tion vanishes Pitxj)

=

the contact,

the surface trac-

(18)

O

on aV, outside the contact (d > 0). (3) Inside the contact we assume surface traction pn is compressive

that

component

of the (19)

Pntxj) G O

inside the contact (d = 0). There are models in which isfied.

the normal

eqns. (16),

(18) and (19) are not all sat-

Fig. 3. Construction of the deformed distance.

1.1. Normal adhesion Equation (19) does not hold when the surfaces of the contacting bodies attract each other. Let a, be the force of attraction that the surfaces exert on each other. a, depends crucially on d d>O ifd>Oifd=O

Pn = a, -pn
(20)

65

Equation (20) means that a compressive traction is superimposed on the force of attraction where there is contact. If a, depends only on position and d, and when it is a monotonically decreasing function of the latter, then we can, if d > 0, express d(xi) as a monotonically decreasing function of p,&) d(xi)

=

d{Pn(xi)l

ford

> 0. It should be noted that eqns. (18) and (19) are violated

(21)

by eqn. (20).

1.2. The model of Oden and Martins Oden and Martins [6] regarded the surface zone, which is dominated by asperities, as a layer of finite thickness; by means of this layer they attempted to model asperity behaviour. They assumed the following relations: - h)+‘“n = - c,(-

pn = - c,(Au,

d)+“n

(22)

where f+ = 0 if f < 0, f, = f when f > 0 and the material constants m, are greater than zero. This relation replaces the non-penetration tion (eqn. (16)) and gives a definite connection between d and pn .

c, and condi-

1.3. Slip In all the laws a crucial role is played by the slip S. It is given by S, = Au, + w,

(23a)

(see at position .lcfE i3V, and time t, where . is the material differentiation eqn. (5)) and wI is the velocity of the coordinate system-reference state of body 1 at (3cI, t) with respect to that of body 2. We can discretize the slip. To that end we can assume that the calculation has proceeded till time t’ and that we are engaged in calculating the displacement-strain-stress field at time t > t’. Then we have, approximately, (t -

t’)S, = Au, - AU; + (t -

t’)w,

(23b)

at position xi E aV, and time t, with Au: the surface displacement difference at position 3ti and time t’ < t. Further frictional conditions will be introduced with the models that are reviewed. 2. Simplified

theory

Simplified U&I)

= LP&l)

theory

is characterized

by the relations (24)

(no sum over a; xi E i3V,) where u,, is defined in eqn. (23a), L, is a material constant called the flexibility and par is the tangential surface traction. This model is isotropic and it may be visualized as a brush of stiff needles (see Fig. 4). It is seen from eqn. (23) that we are interested in Au, and therefore using eqn. (14)

66 “2r

“21

-.

P2r

P2r

/’

~

Fig. 4. Simplified

theory.

Au, = u17 - u2T= Lp,

(25)

with L = L1 + L2. All models of friction may be combined using the brush model of simplified theory. Here we consider the following macroscopic friction model, where “macroscopic” indicates that we do not operate on the asperity level. IPTI
=dxis

Is71,P"e)

CW

if IS,1 # 0, p = -@‘,/I&l )gand where g=g(Xi, &I, p”,) is the traction bound. p”, is a so-called mollification of pn

f%W = JJ C(r) exp lY,K

PAX,

+

Y,) dyl dy2 if r > 0

I”

ifr=O.

= P”(X,)

(27)

where r is some positive constant to be chosen K?r)F”

= JJ exp{r2/(y,y, lu,K

- r2)I dyl dy2

r>O

r

C(r) is a normalization constant that ensures that j&(x,) = p,(x,) as r + 0. for r > 0 has been used by Oden and Pires [7] to establish the existence of the elastic field (equal to displacement, stress and strain) due to combined normal and frictional conditions (eqns. (16) - (19), (23), (26) and (27)) in the case of linear elasticity. We wish to use the model (eqn, (26)) to find the frictional traction at time t; pr when the normal traction pn and hence its mollification & are known, together with traction pi, slip SC and displacement difference AU:, all at time t’. We consider only the case that the traction bound is independent of the slip; the case that the traction bound depends on the slip is under investigation. We have using eqn. (23b)

&, for r > 0 is infinitely many times differentiable. &

(t -

t’)Sr(xi,

t) = AU~- AU: + (t - t’)~,

= UP, -pi) We definep,,

by

+ (t -

t’)w,

(28)

67

t’)w7

(t -

0 = L@HT

-p:)+(t--')w,-PH7=P:-

L

(29)

If

Indeed, for this value of p7 the slip (t - t’)& vanishes. If

lpH,I >g-P,

Wb)

= -

while (t - t’)s, = Lcp, -p:> = L@, =

L@,

+ (t - t’)W,

-PHT)

+ L@H~

-PH7)

= LP,

---Pi>

(

1 -

+ tt

-

t’)w,

$d)

from which it is seen that S, is precisely opposite the slip (lpn,l > g). This is the Fastsim algorithm (Kalker [ 21). It is very fast in operation.

3. The Panagiotopoulos process When the normal pressure is unknown, both in simplified theory as well as when the complete elastic equations are taken into account, we work in the following way. (1) Assume that the tangential tractionp,(0) vanishes (let I = 0). (2) We determine the normal traction pntn with prtn as the tangential traction. Next, we determine the tangential traction p7(‘+‘) with pnQ as the normal traction. If p7(‘+l) is close enough to p7(0 we stop, otherwise we set I = I + 1 and restart at “I”. This is the Panagiotopoulos process (1975). It was used by Oden and Pires [7] to prove the existence of the elastic field in the elastostatic case and by myself to perform calculations in the half-space elastostatic case. I found that the process sometimes converged to give a reasonable solution; however, in other cases it diverged. Kikuchi [8] had the same experience. We believe that the reason is the non-uniqueness of the solution. The existence of the elastic field has been proved but uniqueness applies only for small coefficients of friction. A Panagiotopoulos process terminated after one step was introduced by Johnson [ 91 as an approximation.

68

4. Symmetry

(De Pater [lo])

In general, the Panagiotopoulos process requires an infinite number of steps to converge. It is of interest to consider cases where the Panagiotopoulos process terminates after a finite number of iterations. One iteration can suffice in the case of symmetry. It takes place when dV, is a plane about which the bodies 1 and 2 are geometrically symmetric and also elastically symmetric which means that the bodies 1 and 2 have equal elastic constants. Consider a loading of the bodies which is mirror symmetric with respect to the plane aV,. Then, by Newton’s third law, pli =-p2i and by symmetry (see Fig. 5(a)) P, = Pin = -Pzn PI7 = P27 = -P27

= 0

u, = Uln = -U2n

(31) Au, = 24,

U17 = + u27 -

Au, = 0 Consider now a loading which is mirror antisymmetric with respect the plane a&. Then, by symmetry and Newton’s third law (see Fig. 5(b)) -pzn

Pln = P2n =

PT = Plr Uln

=

to

= 0

= -p27

(32)

U2n

u, = U17 = -

u27

-Au,=0

Au, = 24 of the fields of eqns. (31) and The total field follows by superposition (32). It is seen that the tangential traction p7 does not influence the normal quantity Au,. So we can calculate the symmetric field first without reference to the tangential traction. p,, is then known and the tangential traction can be calculated. A technical application is the wheel-rail problem in railways where both wheel and rail are regarded on elastically symmetric half-spaces when performing elasticity calculations or where simplified theory is applied to the tangential problem (eqn. (32)). It should be noted that simplified theory cannot be applied to the normal problem (31)!

UT=",T=U>T (a)

I

Pl” = P”

PT=o

------)

PT =Plr=-Pzr

u1n=u, c

LIZ"=-U"

P,,=

’ P” = 0 4

-P,

3°C

* "T=U,T=-U2~

@I

Fig. 5. Geometric and elastic symmetry. antisymmetric loading.

! u,=u,"=u~, I

a",

t

(a) Mirror symmetric

loading and (b) mirror

69

5. The principle of virtual work and its dual We will state the principles of complementary virtual work formation is described by eqns. (26). They are, in a weak sense, inequalities :

equations of equilibrium

uij, j + fi -pi&=0

(Jfj = ehk

ujf9

= Sfjhk

= Efjhk uij,

ehk = Eijkk uh, k

Efjhk

= Ehkfj

=_

‘2 @h,

Uf = iii

stress-strain relations

= Ejfkk

similar

S fjhk ehk

ofj

of virtual work and its dual, the principle for bodies in contact, where the contact (16), (18) and (19) and friction by eqn. equivalent to the following equations and

@W valid in v,

Wb)

I

compatibility equations :

k + uk, h)

(33c)

in a&,

Pi = Oijnj, Pi = Pi d = h -Au,

(33d) in aV,,

(33e)

> 0,6 = -@Au,)

non-penetration

Pn GO

compressive normal traction

p,d=O

outside contact, pn = 0 traction bound

IPA
(t - t’)S, = (t - t’)~, + Au, - AU; ISI f 0 -p7=-

slip definition

SS,

in area of slip I

Is,I

contact formation (a&) frictional conditions (aI%)

(33f) (33g) (33h) (33i) (33j) (33k)

In addition, very often the following relation is valid for g: (34)

g = g@f)

is a known function of position. The principle of virtual work reads (Duvaut-Lions [ 31, in classical notation) SQ
v 6Uf

(35a)

substitute eqns. (33b), (33c), (33d), (33f) and (33j)

SQ = .z,

${-UfjSUf,j ’

+

(ft

- piif)6uf}

dV + $ PfGuf dS - $ a”,,

[ “a

I

gSlS,l dS

avc Wb)

= Z2

-

J (ufj,i

+

ff -piii)SUf

i v,

s IPnNAu,) + (t - t’)@,&

av,

@f -pf)bUf dV + J av,, + gSlS,lN

@

dS t (35c)

70

where dV is element of volume and dS is element of area. If we substitute eqns. (33b) and (33~) into eqn. (35b), and consider the static case I&i = 0 without body force (Ifi = 0) and observe eqn. (34), it can be shown that eqn. (35) is equivalent to sup Q= U

c 0=1.2

is

LE s ijkk

Ui,fUh,

k

&ut dS - f

dV + s

Vcr

1

av,,

(t - t’)glS,f dS

av,

(36)

substituting eqns. (33b), (33d), (33f), (33j) and (34). The function Q of eqn. (36), called the work function, is concave and it can be shown (Duvaut and Lions [ 31) that the solution exists and is unique. The variational inequality (eqn. (35)) has a dual which is called the complementary virtual work. It is also equivalent to eqn. (33) and reads 6C>O

voa,,

(37a)

substitute eqns. (33a), (33b), (33e), (33g), (33i) and (33j)

Sg = 0 if eqn. (34) holds. Just like the virtual work variational inequality (eqn. (35)) the complements virtual work inequality (eqn. (37)) can be regarded aa the extremality condition of a function C in the case that Zi = 0, fr 5 0 and Sg = 0 infC

(384

Oii

substitute eqns. (33a), (33e), (33g) and (33i)

+

s [-

hp, + ((t - t’)i.uf - AZ&,]

dS

av,

The function C, called the complementary

work function, is convex.

(38b)

71

6. Numerical methods based on the problem formulation of Section 5 Using eqn. (33b) we have +E fjhkeijehk

= f%,,Utj

= +%jhkUijUhk

substituting eqns. (33b) and (33~). Hence we can integrate the volume integrals of eqns. (36) and (38b) partially

=-

s

iUlUfj,j

dV + J

V,

ill#,

dS

%

substituting eqns. (33b) and (33~). If we build up our solution by superimposing elastic fields that satisfy the static equations of equilibrium without body force ((T~],~ = 0) and the constitutive relations (eqn. (33b)) then the volume integrals of Q and C can be written as a surface integral LE s 2

ijhkeifehk

dV

va

=

s

+&hkUiphk

V,

dV

= s

$wi

m

a&

and the principles of maximum work and minimum complementary work

dS

+

s [(+n

- h)p, + (+Au,

- Au; + (t - t’)w,)p,]

PW

ds

W’b)

av,

The importance of the so-called surface mechanical formulation (eqn. (39)) is that the dimensionality of the problem is decreased by unity. The difficul-

72

ty is how to obtain the basic solutions of the elasticity equations with which the problem is solved. In the technically important case of the half-space geometry (wheel-rail contact) such solutions can be obtained from the Boussinesq-Cerruti integral representation (see Love [ 111 for this representation and Kalker [ 41 for the basic solutions). In order to solve eqns. (39a) and (39b) they must be discretized. We will not enter here into a discussion of that problem. Then, in the case of a twodimensional elastic problem, both representations lead to so-called quadratic programming problems, that is, minimization problems consisting of a quadratic objective function (i.e. a function to be minimized) under linear equality and inequality constraints. In the case of a three-dimensional problem the objective function of the complementary work C remains quadratic but one set of non-linear constraints (p,p,
IS,1 = (s,s,)“*

= $le

“2’(S,S, + ep* def

z ti

=-

2E

t&s, 2E

= I&I -(E/2)

if IS,] < e ) foracertaine>O if IS,] > E i

(40)

It should be noted that tile is often infinitely differentiable in the argument S, while GZe is continuously differentiable. With these regularlzations the maximization of the work function Q for a fixed regularization parameter E under linear equality and inequality constraints can be performed fairly easily. Strictly speaking, E should then be diminished and Q should be maximized again with the solution corresponding to the previous E used as a starting value, and so on, until convergence occurs. This is called a sequential process [ 121. However, in the case of concavity of the function Q e can be calculated beforehand in such a manner that the error in Q does not exceed a pre-assigned value (Fiacco and McCormick [ 121). If the error in the work function is regarded as a measure of the error in the elastic field this implies that in our case the maximization need be performed for a single value of E. This value may, however, be so small that the convergence of the maximization is endangered.

73

A Panagiotopoulos process is still required while also the method is not suitable for dynamic problems (]&I + 0). Oden and Martins [ 61 have proposed a way of dealing with the dynamic contact problem which also does away with the need for a Panagiotopoulos process. The basis of this method is the constrained variational inequality of virtual work (SQ < 0, see eqn. (35)) which they reduce to an unconstrained variational equality as follows. They replace IS,1 by ele (see eqn. (40)), they set ufi = Ellhk~h, k (see eqn. (33b)) whereby they take account of eqns. (33b) and (33c), they substitute gi for U[ in aV,, (see eqn. (33d)) and they substitute S, according to eqn. (33j). As eqn. (33f), they state that the surface aV, should be regarded as a surface zone, i.e. a surface layer of small but finite thickness, to take asperities into account. Equations (33f) and (33i) are replaced by constitutive relations of the form pn = -t&(--d)

mn

g = c+~)+~T+

(41)

d=h-Au, where f+(x) = f(x) when f(x) > 0, f+(x) = 0 when f(x) G 0 and c,, cT* m, and MT are positive material constants. Then the principle of virtual work becomes SQ=O substituting

tid6ui eqn. (33j), eqn. (41) substituted

(42) which upon discretization in the sense of the finite element (FE) method yields a system of non-linear ordinary differential equations which have been solved numerically by Oden and Martins in a number of cases. As the problem (eqn. (42)) is not a concave maximization they should apply a sequential process of diminishing E until convergence occurs. They did not do so, however. 7. Plasticity models of friction Plasticity-based models of friction have been considered by Frederiksson [ 51, Michalovski and Mroz [ 131 and Cumier [14]. We have not been able to incorporate these theories in the scheme of Section 6. As an example, we present the theory of Frederiksson which despite errors (e.g. the proof of his uniqueness theorem is false) is considered a pioneering work.

74

The theory is incremental and the increment of the quantity a is denoted by Da. The theory holds for slip hardening which is defined by the following two properties. (1) Positive work is done by the load increment during application (work is defined as positive when fed into the system studied). (2) The net work performed by the load increment during a cycle of application and removal is non-negative. We assume that the reference states do not move with respect to each other w, = 0

(43)

Then, (1) and (2) imply that (t - t’)S, = (t - t’)w, + Au, - Au:, = DAu, - (t - t’)&Dp,

= - Dp,DAu,

> 0

(44) (sign, see (3) below)

Also, the following basic assumptions are made. (3) The slip increment DAu,of body 1 with respect to body 2 depends linearly on the load increment D~,.Jexerted on body 1 DAu, = - h, @Dpp (slip rule)

(45)

h,p is a material constant. (4) We call g&) the slip function. We have g@t)GO g&)=0 lDAu,l f 0g@)

3s -pi>0 apt

(46)

= 0 is called the slip surface. Then it can be shown that

(47) L > 0, g(pl) determined by the material. In the absence of slip hardening the material constant L is a positive indeterminate function of position which is determined by the problem. We obtain a numerical FE method by observing that the stress increment Du,~ satisfies the equations of equilibrium, Doui + Dfi = 0 in V,, a = 1, 2 and that Duf are the corresponding displacements. They may be formulated in a weak sense as a variational equation

75

0 = C

J (DUij,j + Dfi)GDUi dV

a=1,2

Vo

Now, we assume that Doij and DuI,j satisfy the constitutive relations of linear elasticity as follows: on SV,, , 6DUi = 0; On 6 I&, , Dpi = Dlj,, prescribed; on aVK (the contact area which we consider prescribed) 6DAu, = 0; on aVA (the adhesion area) 6DAu, = 0, if we consider the area of adhesion prescribed. In 6 Vs , the area of slip, the slip increment DAu, satisfies eqn. (45). This gives

+

J

Dp,&DAu, dS

a%

(48)

substituting DUij = Eijhk DU,J in V, DAu

=

cl

_

1 miaPmglaPp) L miap,)(wap,)

We solve eqn. (48). Then we must determine the contact area and the area of slip and adhesion. To that end we start by guessing which FE nodes are in these regions. Then If xi E area of slip, and g < 0, then Xi is transferred to the area of adhesion. If Xi E area of adhesion, and g > 0, then x3 is transferred to the area of slip. If Xi E contact area, and p3 > 0, then contact is lost. If Xi is outside the contact area, and d < 0, then contact is established.

l(49)

We used a variant of the algorithm (eqn. (49)) for three-dimensional elastic contact problems (1981) but found that in some cases of steady state rolling with spin the method did not converge. Frederiksson applied his theory to the essentially simpler twodimensional and axisymmetric geometries and he reports no failures. This throws doubt on the general validity of the criteria (eqn. (49)) which is lent colour by the circumstance that Frederiksson’s uniqueness proof breaks down precisely because the variation in the areas of contact, adhesion and slip is not taken into account.

76

For example

we consider

Coulomb’s

law. We take (50)

Then, the associated DAu, = -

slip rule (eqns. (45) and (47)) becomes

1 PaPp

- dPP L IPA2 L > 0, a material constant. In the case of no slip hardening minate-positive and the associated slip rule is

(51) L is indeter-

DAu, = - Xp, X20 xg=O x

=

(52)

1 PPd& L

IPA

which is Coulomb’s law. Frederiksson applied problems.

his theory

to twodimensional

and axisymmetric

8. Conclusion The static Coulomb-type friction law of Duvaut and Lions can be treated both with the non-variational simplified theory (large computing velocity) and with a complete variational theory (small computing velocity). A proof of the numerical methods can be given. The dynamic Coulomb theory has been treated by Oden and Martins by reducing the constrained virtual work variational inequality to an unconstrained variational equation by the use of regularization and penalty methods. This necessitates a sequential method but this aspect is omitted by Oden and Martins. A proof of the method is not provided. Frederiksson gives a plasticity-based theory and a numerical method for slip hardening which he applies to two-dimensional and axisymmetric problems. He gives an attractive method for the determination of the areas of contact, adhesion and slip which is successful in all cases he considered but which fails in some cases of three-dimensional steady state rolling with spin.

References 1 J. J. Kalker, Simplified theory of rolling contact, Delft Progretw Rep. Ser. C 1, 1973, pp. 1 - 10 (Delft Universityof Technology, Delft, The Netherlands).

2 J. J. Kalker, A fast algorithm for the simplified theory of rolling contact, Veh. Syst. Dyn., 11 (1982) 1 - 13. 3 G. Duvaut and J. L. Lions, Les Zn6quations en M&unique et en Physique, Dunod, Paris, 1972. 4 J. J. Kalker, The computation of three-dimensional rolling contact with dry friction, Znt. J. Numer. Methods Eng., 14 (1979) 1293 - 1307. 5 B. Frederiksson, On elastostatic contact with friction, Ph.D. Thesis, Linkiiping, 1976. 6 J. T. Oden and J. A. C. Martins, Models and computational methods for dynamic friction phenomena, Comput. Methods Appl. Mech. Eng., 52 (1985) 527 - 634. 7 J. T. Oden and E. B. Pires, Non-local and non-linear friction laws and variational principles for contact problems in elasticity, J. Appl. Mech., 50 (1983) 67 - 76. 8 N. Kikuchi, personal communication, 1982. 9 K. L. Johnson, Tangential traction and microslip in rolling contact. In J. B. Bidwell (ed.), Rolling Contact Phenomena, Elsevier, Amsterdam, 1962, pp. 6 - 28. 10 A. D. De Pater, On the reciprocal pressure between two elastic bodies. In J. B. Bidwell (ed.), Rolling Contact Phenomena, Elsevier, Amsterdam, 1962, p. 33. 11 A. E. H. Love, A Treatise on the Mathematical Theory of Elasticity, Cambridge University Press, London, 2nd edn., 1926. 12 A. V. Fiacco and G. P. McCormick, Non-linear Programming, Wiley, New York, 1968. 13 R. Michalovski and Z. Mroz, Associated and non-associated sliding rules in contact friction problems, Arch. Mech., 30 (1978) 259 - 276. 14 A. Curnier, A theory of friction, Znt. J. Solids Struct., 20 (1984) 631 - 647.