INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING, VOL 38,2603-2617

(1995)

A NOVEL FINITE ELEMENT FORMULATION FOR FRICTIONLESS CONTACT PROBLEMS PANAYIOTIS PAPADOPOULOS, REESE

E. JONES AND JEROME M. SOLBERG

Department of Mechanical Engineering, University of California at Berkeley, Berkeley, C A 94720, U.S.A.

SUMMARY This article advocates a new methodology for the finite element solution of contact problems involving bodies that may undergo finite motions and deformations. The analysis is based on a decomposition of the two-body contact problem into two simultaneous sub-problems, and results naturally in geometrically unbiased discretization of the contacting surfaces. A proposed two-dimensional contact element is specifically designed to unconditionally allow for exact transmission of constant normal traction through interacting surfaces. KEY WORDS:

frictionless contact; large deformations; finite elements

1. INTRODUCTION Finite element methods are used extensively in the solution of contact problems. From a purely computational standpoint, detection of contact and subsequent satisfaction of the impenetrability constraint are the two key issues to be addressed in the development of a general algorithmic framework. Numerous methodologies have been proposed in the literature of computational contact mechanics since the early works of Conry and Seireg,' and Chan and Tuba.' A fairly comprehensive survey on the topic is found in Reference 3. The present work is concerned with the development of finite element methods suitable for the solution of two-body contact problems in the presence of large motions and deformations. This class of problems is of particular significance in numerous practical applications, such as metal forming processes and vehicular crash analyses. Various commercial and research computer codes employ algorithms for the solution of such problems. Lagrange multiplier m e t h ~ dand ,~ their regularizations (penalty and augmented Lagrangian methods596)are typically used in enforcing impenetrability. The choice of integration method for the work-like term associated with the contact tractions in the weak form of linear momentum-balance plays a pivotal role in the construction of contact elements. Use of nodal quadrature involving the contacting nodes of one of (resp. both) surfaces yields the standard one-pass (resp. two-pass) node-on-surface algorithm^.^ Other integration rules are also applicable, provided there exists a continuous discretization of the contact interfa~e.~' The main contribution of the present paper is in the identification of a general procedure according to which the two-body problem is approached as a sequence of two simultaneous sub-problems. As in the traditional two-pass algorithms, the surfaces of both interacting bodies are used in the analysis without need for introduction of an (often arbitrarily chosen) CCC OO29-5981/95/152603-15 0 1995 by John Wiley & Sons, Ltd.

Received 1 Jury 1994 Revised 29 October 1994

2604

P. PAPADOPOULOS, R. E. JONES AND J. M. SOLBERG

intermediate contact surface. The main advantage of the proposed approach over the two-pass node-on-surface algorithms is that it allows for a straightforward interpretation of the integration rules used on the contacting surfaces and, for appropriate choices of admissible fields, permits the exact transmission of constant pressure from one body to another. In the spirit of the patch test originating in the work of Irons," and its subsequent generalizations, capability for exact representation of constant pressure (in both magnitude and direction) is viewed as a necessary condition for robustness and convergence of the overall contact algorithm. A brief exposition to contact mechanics is presented in Section 2, with particular emphasis on formulations to be used in the ensuing algorithmic developments. A two-dimensional contact element is proposed and analysed in Section 3, while the results of selected numerical simulations using this element are presented and discussed in Section 4. Concluding remarks are given in Section 5. 2. THE TWO-BODY CONTACT PROBLEM Consider bodies W", a = 1,2, identified with open and connected sets R" in linear space R3, equipped with canonical basis (El,E2,E3)and the usual Euclidean norm. At least one of the bodies is assumed to be deformable. A typical material point of W"in the reference configuration is algebraically specified by vector Xu.The boundary of each body is represented by set an" and possesses a unique outer unit normal N" at each of its points. The motion of W"is described by the mapping x" on R" u 8R", so that the position vector x" of material point X" in the current configuration is given by x u = x"(X",t ) at each time t, and the displacement vector u" is defined according to U"(X", t ) : = xyxa, t ) - X" The mapping x" is assumed smooth throughout its domain and invertible at least on nu.Body a and its boundary in the current configuration (at time t) are identified with 0;and aCJy, respectively, hence a;:= x"(n",t) and an;:= x"(an", t). Also, the outer unit normal to 8 0 ; is denoted by n'. The motion of any system of bodies (including a single body) is subject to the principle of impenetrability of matter, as stated by Truesdell and Toupin in Reference 11 (p. 244). This implies for the two-body problem that at all times n R: = 8

(1)

At any given time, the two bodies are said to be in contact along a subset C of their boundaries if, and only if, an: n aa::= c # 8 It follows from the above definition that the boundary of each body can be uniquely decomposed into three mutually exclusive regions according to

an; = r: u r; u c where Dirichlet and Neumann boundary conditions are enforced on r: and r;, respectively. Although not explicitly noted, it should be clear from the above that r:, r; and C generally depend on time. Gap functions g'"), possibly multi-valued, can be defined on the boundary of each body as follows: for each x2 E an:, g('):dR: x an: H R is given by g(I):= (x2

- x').

n1

(24

FRICTIONLESS CONTACT PROBLEMS

2605

Figure 1. Definition of gap functions

where x' = x'(x2; n') is such that (x2 - x') x n' = 0, see Figure 1. Convexity of an' renders g") single-valued, although such a restrictive geometric condition will not be imposed at the outset. A completely analogous definition for g''):dQ: x 80: H 53 yields

-

g(2)., (x 1 - x2) n2

(2b)

where, again, for each x1 E an:, x2 = x2(x'; n2)satisfies (x' - x2)x n2 = 0. Defining equations (2a) and (2b)imply that gap functions g'') and 9") are identically equal to zero on C, namely that g y c ; an:) = g'Z'(C;

an:):= g(C) = 0

(3)

Consequently, impenetrability condition (1) can be rewritten in terms of the above gap functions as g'".

2 0, g'2' 2 0

At the absence of inertial effects, the local form of the equations governing the motion of each body LY are given by divt" + pub"= 0

in Q:

(44 (4b)

tun" = tyn.,

r: on r;

g'") 2 0

on

(4d)

u"=

Ua

-

on

(44

where t" is the Cauchy stress tensor, pa the mass density, b" the body force per unit mass, U' the prescribed boundary displacement, and the prescribed traction vector on r:. Application of the standard weighted-residual method, in conjunction with the introduction of Lagrange multiplier p 3 0 for the impenetrability constraint, results in the weak form of the equations of motion, which states that the displacement solution u" of equations (4) and the Lagrange multiplier field p satisfy

q,.,

2606

P. PAPADOPOULOS, R. E. JONES AND J. M. SOLBERG

for all admissible function w" and q. Without loss of generality g(2)is used in equation (5a) for the definition of the gap function on C. Displacement fields u" belong to spaces Q" with %" = {u" E H1(fZa)luu= u" on

r:}

and weight functions w" belong to spaces P'defined as W' = {w" E H1(Ru)lwu= 0 on

r:}

The admissible functions q 2 0 are piecewise continuous. To prove that equation (5a) holds, note that the work done by contact forces along C on admissible functions w1 and w2 is given by r

r

At the absence of friction and recalling that n1 implies that t(,l) 1 =

=

- n2 on C, Cauchy's lemma on the stress vector

- t&2):= - Pn'

(7)

With the aid of equation (7), equation (6) is written as n

p(w' - w2).n'dy = -

W, = JC

Jc

p(wl - w2).n2dy,

which shows that the Lagrange multiplier field is naturally identified with normal traction (pressure) along the contact region. The pressure field p is generally assumed to be only piecewise smooth, thus allowing each body to feature material interfaces in the neighbourhood of C. Moreover, since

(8)

pg'") = 0 on

inequalities (5b) follow from (4d) and the assumed non-negativeness of q. In order to further clarify the role of the Lagrange multipliers in the two-body contact problem and provide some motivation for the ensuing numerical approximations, use (2b) and (8) to rewrite the integral in (5b) on surface C as

r

r

The notation C(") is employed to merely emphasize that the contact surface C can be selectively viewed as part of aQ;, as indicated in (3). The integral expression (9) suggests that, given a properly defined Lagrange multiplier field p(")on the boundary of each body, integration of (5b) can be performed separately on each of the contacting surfaces. However, it is clear that fields p(') should satisfy balance of linear momentum on the (common) contact surface C@),namely that p(l)=

P(2)

The above observation will be exploited in the approximate solution of (5a) and (5b).

2607

FRICTIONLESS CONTACT PROBLEMS

Following the procedure used in the derivation of (9),equations (Sa) and (5b) are rewritten as

2

r

for all admissible wa and q'"), where q'l) = q('). The classical penalty regularization of (Sa) and (5b) is obtained by setting P=&<-g) in the last integral term of (Sa), where is the Macauley bracket and E > 0 is a constant penalty parameter. Constraint conditions (4d) are relaxed (thus permitting controlled penetration to take are determined so that place) and displacement fields u: E ( a )

(11) for all w' E W'. Under restrictions such as formal equivalence of the original unconstrained problem to a convex minimization problem, the sequence of solutions u: as E + co can be shown to converge to the solution of the constrained problem (S).lz In practice, penalty methods are used successfullyeven when the equilibrium state is not associated with minimization of the total potential energy. Due to the occurrence of penetration for finite values of the penalty parameter E, a unique definition of the contact surface C and gap function g is not readily available for use in (1 1). A crucial step in the proposed formulation is the introduction of a well-defined penalty regularization for (lOa) and (lob), such that p'") = &( - g'"')

(12)

With the aid of (12), equation (10a) becomes

Each of the last two integrals of (13) is identical in form to the boundary integral emanating from penalty regularization of a Signorini problem. Although in the present context these two integrals are clearly coupled by the definitions (2a) and (2b) of gap functions g@),this decomposition has important computational implications, as will be discussed in the next section. The developments presented here can be easily extended to encompass the multi-body contact problem by reducing it to a series of coupled two-body problems. 3. APPLICATION TO TWO-DIMENSIONAL CONTACT The remainder of this article is devoted to the discretization of two-dimensional contact problems based upon the penalty formulation and the identification of gap functions on both surfaces, as

2608

P. PAPADOPOULOS, R. E. JONES AND J. M. SOLBERG

Figure 2. Determination of contact elements by projection

suggested in Section 2. The two key issues to be addressed here are the geometric construction of contact elements and the choice of admissible fields for the finite element approximation.

3.1. Discretization of the contact surfaces A continuous discretization of the contact surfaces is advocated. Surfaces C(')and C2)are uniquely determined as those on which g") < 0 and g'2) < 0, respectively. Setting aside implementational details, contact surface C") is discretized by a series of normal projections (not necessarily closest-point) from an: to an:, as shown in Figure 2. An analogous procedure is followed for the discretization of C2).Consequently, each contact element C!.) relates a single spatial element edge on surface C(') to the opposite surface. The apparent non-smoothness of the discrete boundaries results in discontinuity of unit normals and gap functions. No attempt is made here towards circumventing this problem, although its effect might not be negligible, especially in problems of rolling contact. For a special discretization that results in smooth boundary representation in two-dimensions, see Reference 13.

3.2. Finite element jieids A specific contact finite element is suggested in this section, based upon equations (loa) and (12). The choice of finite dimensional fields is guided by previous works especially on the Signorini p r ~ b l e m . ' ~ In . ' ~this context, it has been formally demonstrated that appropriate choices of the admissible displacement and pressure fields lead to unilateral contact elements that are able to stably replicate the impenetrability condition (i.e. they satisfy the underlying LBB condition) and are accurate. Such convergence analysis, although not currently available for the kinematically non-linear two-body contact problem addressed here, provides a guideline for the selection of admissible fields. Numerical integration is typically employed for all boundary terms on the contact surface. Discounting discretizations that use straight-edge spatial elements (e.g. three-node triangles and four-node quadrilaterals), the gap functions g(') within a single contact element are non-polynomial with respect to local boundary co-ordinate systems. Thus, numerical integration on C(') is generally inexact and all integration rules introduce errors that directly influence the formulation. The proposed contact element is based on displacement fields emanating from standard Q9 (nine-node isoparametric) elements and employs Simpson's integration rule.

FRICTIONLESS CONTACT PROBLEMS

2609

Element (C3-2) has a quadratic assumed pressure field defined by means of the isoparametric co-ordinate system tt)on the boundary edge of the spatial element, namely

1

p("' =

L2*i(tk"')pl"'

i = I. c, r

where Lz,iare the Lagrangian interpolation functions of degree 2, and subscripts 1, r and c denote the ordinates of a field at the left, right and center of element Cf), respectively.Use of (12) implies that

i =I,c,i

By construction, the element preserves continuity of the pressure fields assuming smooth boundaries and constant penalty parameter E. For each element C?), the contribution of contact forces to balance of linear momentum is obtained from (10a) as 3

where 3

W" =

1 N,(~F')w: I=1

In the above interpolatior, NI=[" 0

"1

NI

and N I , I = 1, 2, 3, are the standard one-dimensional quadratic isoparhmetric shape functions. The proposed identification of the original two-body problem with a series of two simultaneous sub-problems presents two definite advantages to traditional formulations: first, it renders the distinction between 'master' and 'slave' surfaces obsolete by removing the geometric bias in the representation of the contact surface without introducing an interpolated intermediate surface. Second, it allows for a clear analysis of the constant normal traction patch test suggested in Reference 16, regardless of the spatial discretization of the interacting flat surfaces.Element (C3-2) unconditionally passes the aforementioned test, since Simpson's rule is exact for the (at most cubic in 5:)) integrand of (14), given that both p(")and nu are constant. Contrary to the traditional implementations (e.g. that of NIKE2D, Reference 9,the proposed framework yields two-pass algorithms that require no superposition of the response in the two passes, but, instead, naturally recovers the overall response in a two-step process. One of the key characteristics of the proposed formulation is that it produces pressure fields that are not a priori equal and opposite (pointwise or in resultant form) on both contacting surfaces. However, pressure fields p(l) and p(') converge to a common limit with mesh refinement, as the two surfaces conform more accurately to each another. An algorithm based on the preceding formulation is symbolically summarized in Box I. Newton method is employed for the iterative solution of the global non-linear problem. The

2610

P. PAPADOPOULOS, R. E.JONES AND J. M.SOLBERG

determination of the contact stiffness matrix which is consistent with (14) is discussed in the appendix. 1. 2. 3. 4. 5. 6. 7. 8. 9.

Loop over a Loop over Ct) Compute g(Ol) If g‘“’ > 0 go to 8 If g(’) < 0 then Form the element residual vector and stiffness matrix and append to the corresponding global matrices Endif Endloop Endloop

4. NUMERICAL SIMULATIONS The proposed formulation has been implemented within the environment of FEAP, a generalpurpose non-linear finite element program documented in Reference 17 and 18. Selected numerical simulations have been conducted using a compressible isotropic neo-Hookean hyperelastic model for all deformable bodies. The strain energy functional Wof the model is defined by W=

f p ( ~ C-

3) - pln111;”

+ $ A ( z z z ~-’ ~1)2

where C is the right Cauchy-Green deformation tensor, ZC = tr C,ZlZc = det C and A, u , are the Lame constants. Nine-node quadrilateral isoparametric elements are used for the spatial discretization of all interacting bodies.

4.1. A patch test infinite deformations

A system of two rectangular blocks with identical elastic properties is subjected to constant distributed load jj in the reference (undeformed) configuration, as in Figure 3. This problem is specifically designed to test the capability of the proposed contact elements to exactly transmit constant normal traction regardless of the spatial discretization and the (constant) value of the penalty parameter. Indeed, due to the nature and magnitude of the uniaxial load and the assumed constitutive response, the contact surface should remain flat and the computed traction field constant for any fixed value of the applied load, even at finite strains. Element (C3-2) unconditionally passes the above test, see Figure 3 for an application to an unsymmetrically distorted mesh. As argued in Reference 16, traditional node-on-surface algorithms typically fail the test (except for special cases), owing to lack of accuracy in the integration along the contact surface. 4.2. Contact between deformable cylinders

Two cylinders of radii R1 = 8.0 and R2 = 18-0in plane strain are pressed against each other by two flat rigid surfaces, as in Figure 4. A uniform vertical displacement t i z @ )= t is imposed on the upper rigid surface, while the lower rigid surface remains stationary. The material properties of

261 1

FRICTIONLESS CONTACT PROBLEMS

.............................

................................ P .

................................

-

I '*

I

"

P

0

I

Figure 3. Finite deformation constant pressure patch test

Figure 4. Contact between deformable cylinders

b th cylinders re

, I= 2.88 x lo2, p = 1-92 x 10'

The pressure fields p(') developed between the two cylinders at times t = 3-0, 5.0 and 11.0 are plotted in Figure 5. The penalty parameter is set to E = lo5 on all slidelines. The computed pressure fields are smooth and nearly identical on the two surfaces. For reference, the finite element solutions are compared to the approximate linear elastic solution due to Hertz," as

2612

P.PAPADOPOULOS, R. E.JONES AND J. M.SOLBERG

"I-'

"

''

I_

- - - C3-2 Hne

40-

20

-

x,-coordinate pmsun on Iowa cylinder at k3.O

c3-2

~

8040

-

20

-

0.0 0.5

1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5

x,-eoordinate.

xlQoordinate

prusue on upper cylinder 01 t=5.O

pressure on lower cylhdu at 1~5.0

I

- c3-2 -

I

I

I

I

I

- - - HerU

120 -

120 -

100

100 -

140

-

80-

80

-

W -

80

-

40-

40-

20-

20

0.0

1.0

2.0

3.0

4.0

5.0

6.0

-

7.0

x,-coordinate pressure on uppac y l i e r at t= l I .O

x,-cwrdinatc pnssun on lower cylinder at tnl1 .O

Figure 5. Contact between deformable cylinders; pressure fields

2613

FRICTIONLESS CONTACT PROBLEMS I

I

I

I

Figure 6. Contact between deformable cylinders; deformed mesh at t = 11.0

Figure 7. Pinching of an elastic ring

2614

P. PAPADOPOULOS. R. E. JONES AND J. M. SOLBERG

Figure 8. Pinching of an elastic ring: deformed mesh at G2 = 6.2

2 16000

I

g

a

I

I

I

I

I

I

g 18000

-

-

12000 -

-

14000

-

8000 6000 4000 2000 -

10000

0

I

I

I

I

I

I

v)

I

I

I

I

I

I

14000 12000 -

-

-

8ooo-

-

60oo-

-

loo00

4000

-

2000

-

-

-

.

pressure distribution on punch

x,-coordinate pressure distribution on interior contact surface

Figure 9. Pinching of an elastic ring: pressure fields

FRICTIONLESS CONTACT PROBLEMS

2615

presented in Reference 20. The plots demonstrate that the deviation of the computed pressure fields from the respective Hertzian fields is increasing with ti2, as expected. The current configuration of the contacting bodies is shown in Figure 6 for tiz = 11.0.

4.3. Pinching of an elastic ring

An elastic ring of inner radius Ri = 6.0 and outer radius R, = 8.0 is pressed by means of prescribed vertical displacementsti2 on the outer edge of two stiff elastic plates, as in Figure 7. The material properties of the ring and plates are

A, = 5.77 x 103, pr = 3.84x 103 and

A, = 5-77x 104, p, = 3-84x 104 respectively. Plane strain conditions are enforced and a penalty parameter E = lo6 is used on all slidelines. The ring is deformed until its interior surface is subject to self-contact. The current configuration of the bodies for tiz = 6.2 is depicted in Figure 8, while plots of the pressure distribution on the interacting surfaces are included in Figure 9. It is observed that the pressure fields are quite smooth, despite the near singularity at the edges of the two plates and the relative coarseness of the discretization. 5. CONCLUSIONS

The geometrically non-linear two-body contact problem has been formulated as a series of two simultaneous sub-problems. It is shown that this decomposition has important algorithmic implications, as it constitutes the basis for contact finite elements which retain the overall structure of elements used for unilateral problems. These elements are easily constructed in a geometrically unbiased fashion, allow for exact transmission of constant normal traction and, generally, exhibit excellent numerical performance. ACKNOWLEDGEMENTS

The research work was supported by the National Science Foundation under contract Nr. MSS-9308339 with the University of California at Berkeley. R. E. J. was also supported by a National Science Foundation Graduate Research Fellowship. APPENDIX: CONSISTENT TANGENT STIFFNESS The tangent stiffness consistent with application of classical Newton method is derived for the integral terms of (13) associated with the contact surfaces. In order to simplify the derivation, assume that the motion takes place on the plane defined by vectors El and E2,and consider a representative integral term from equation (13)written as

Zl:=

jctl)&( -

g('))w'

an1

dy

Admitting a convected clock-wise surface parameterization of C(")by tangent to surface C(') as

(15) define a vector s"

2616

P. PAPADOPOULOS, R. E. JONES AND J. M. SOLBERG

and note that dy = where

I(0 (I is the 12-vectornorm on

(1 S" (1 2 dt'"'

R3. Moreover, E3 XS" na = -

I1so ll 2

so that with the aid of (17) and (18),the integral in (15) takes the form I' =

jC(,,E(

g('))w'.(E3 xs')d((')

Consequently, the Giteaux derivative DA,[Z,] of Il in the direction Au is written as DAuCI11 = jc(,,&DAu[( - g'")lw''(E3

+ Jc,l)~(- g('))w'-(E3

xs')d<'"

XDA,[S'])~((''

Use of (2a) leads to DA,[( - g('))] = - 9( - g('))(DBu[x2 - x ' l - n '

+ g(')n'.DA,[n'])

(21)

where 9(-) is the Heaviside step function. The second term inside the parenthesis on the right-hand side of (21) vanishes identically, given that n' is of constant magnitude (hence n 1 DAu[n'] = 0). In addition, it follows from the definition of g(') that

-

(x' - x ' ) d = 0 so that for fixed

<(')

DAu[(x2 - x').s']

= (Au'

+ A((')s2

- Au').~'

+ (x'

- x')*As'

=0

(22)

where from (16) AS':= DA,[S'] = c?Au'/c?<''' From (21) it follows that DAo[( - g('))] = - &'( - g('))(Au2

where (22) yields A((') =

+ A((')s2

- Au').n'

-

(Au' - Au2) s' - g(')n' As' s'

'S2

Therefore, substituting (24) to (20) results in DA,[I~]= Jc(,'~&'(- g"))((Au'

L

E(

+

- g'")W'

+ A((%' - Au')-n')w'.(E3 '(E3 x As')d<"'

where As' and A<(') are given by (23) and (25), respectively.

xs')d((')

(24)

FRICTIONLESS CONTACT PROBLEMS

2617

The consistent tangent stiffness matrix for element C3-2 is easily constructed by identifying parameters <(a) with the natural co-ordinates (t)employed in the parametric definition of element boundaries of C(I). REFERENCES

1. T. F. Conry and A. Seireg, ‘A mathematical programming method for design of elastic bodies in contact’, J . Appl. Mech., 38, 387-392 (1971). 2. S. K. Chan and I. S. Tuba, ‘A finite element method for contact problems of solid bodies-part I. Theory and validation’, Int. J . Mech. Sci., 13, 61%25 (1971). 3. Z.-H. Zhong and J. Mackerle, ‘Static contact problems-a review’, Eng. Comp., 9, 3-37 (1992). 4. T. J. R. Hughes, R.L. Taylor, J. L. Sackman, A. Curnier and W. Kanoknukulchai, ‘A finite element method for a class of contact-impact problems’, Comput. Methods Appl. Mech. Eng., 8, 249-276 (1976). 5. J. 0. Hallquist, N I K E 2 D A Vectorized, Implicit Finite Deformation, Finite Element Code for Analyzing Static and Dynamic Response of 2-D Solids, Lawrence Livermore National Laboratories, University of California, UCID-19677 (rev. 1) edn, 1986. 6. J. A. Landers and R. L. Taylor, ‘An augmented Lagrangian formulation for the finite element solution of contact problems’, SESM Rep. 85/09, University of California at Berkeley, 1985. 7. J. 0.Hallquist, G. L. Goudreau, and D. J. Benson, ‘Sliding interfaces with contact-impact in large scale Lagrangian computations’, Comput. Methods Appl. Mech. Eng., 51, 107-137 (1985). 8. J. C. Simo, P. Wriggers and R. L. Taylor, ‘A perturbed Lagrangian formulation for the finite element solution of contact problems’, Comput. Methods Appl. Mech. Eng., 50, 163-180 (1985). 9. P. Papadopoulos and R. L. Taylor, ‘A mixed formulation for the finite element solution of contact problems’, Comput. Methods Appl. Mech. Eng., 94, 373-389, (1992). 10. B. M. Irons, ‘Numerical integration applied to finite element methods’, in Proc. Cont on Use ofDigita1 Computers in Structural Engineering, University of Newcastle, 1966. 11. C. Truesdsell and R. A. Toupin, ‘The classical field theories’, in S. Fliigge (ed.), Handbuch der Physik I I I / l , Springer, Berlin, 1960, pp. 226-793. 12. P. G. Ciarlet, Introduction to Numerical Linear Algebra and Optimization, Cambridge University Press, Cambridge, 1989. 13. A. L. Eterovic and K.4. Bathe, ‘An interface interpolation scheme for quadratic convergence in the finite element analysis of contact problems’, in P. Wriggers and W. Wagner, (eds.), Nonlinear Computational Mechanics, Springer, Berlin, 1991, pp. 703-715. 14. Y. J. Song, J. T. Oden and N. Kikuchi, Discrete LBB-conditions for RIP-finite element methods’, TICOM Rep. 80-7, The University of Texas at Austin, 1980. 15. J. T. Oden, ‘Exterior penalty methods for contact problems in elasticity’, in E. Stein K.-J. Bathe and W. Wunderlich (eds.), Europe-US Workshop: Nonlinear Finite Element Analysis in Structural Mechanics, Springer, Berlin, 1980. 16. R. L. Taylor and P. Papadopoulos, ‘A patch test for contact problems in two dimensions’, in P. Wriggers and W. Wagner, (eds.), Nonlinear Computational Mechanics, Springer, Berlin, 1991, pp. 690-702. 17. 0.C. Zienkiewicz and R. L. Taylor, The Finite Element Method: Basic Formulation and Linear Problems, Vol. 1,4th edn, McGraw-Hill, London, 1989. 18. 0. C. Zienkiewicz and R. L. Taylor, The Finite Element Method; Solid and Fluid Mechanics, Dynamics and Non-linearity, Vol. 2, 4th edn, McGraw-Hill, London, 1991. 19. H. Hertz, ‘Uber die beriihrung fester elastische korper’, J. Rein. Ang. Math., 92, 156171 (1882). 20. K. L. Johnson, Contact Mechanics, Cambridge University Press, Cambridge, 1987.