A boundary element approach for contact problems involving large displacements

PERGAMON

Computers and Structures 68 (1998) 315±324

A boundary element approach for contact problems involving large displacements A. Lorenzana *, J.A. Garrido Department of Strength of Materials, Escuela TeÂcnica Superior de Ingenieros Industriales, University of Valladolid, Paseo del Cauce s/n, 47011 Valladolid, Spain Received 20 January 1997; accepted 16 February 1998

Abstract A new algorithm is presented for the boundary element analysis of the two-dimensional contact problem between elastic solids involving large displacements. The contact constraints are not applied node-on-node but node-onelement, using the element shape functions to distribute the geometry, displacements and tractions on each element in the contact zone. Thus, the discretizations performed along the two surfaces in contact do not necessarily have to be the same. The solution procedure is based on the updated Lagrangian approach and the resulting method is incremental. At least one increment of load must be done for each node that changes its boundary conditions, the geometry being updated after each increment. Additional updatings must be done in order to note any signi®cant changes of the geometry that might occur during the loading process. The algorithm guarantees equilibrium and compatibility at the nodes in the ®nal deformed con®guration and allows us to deal with problems undergoing large displacements without being necessary to change the initial discretization of the boundary of the bodies. Only the frictionless static problem is dealt with and the proposed algorithm is applied to two representative cases. The results obtained, when the displacements are small, are in good agreement with the analytical solution, if possible, or with the solutions obtained by other authors. When large displacements are considered, another non-linearity appears and its in¯uence is going to be shown. # 1998 Elsevier Science Ltd. All rights reserved. Keywords: Boundary element analysis; Large displacements; Contact constraints

1. Introduction As it is known, there are mechanical engineering applications where the stress analysis at the interface between two bodies in contact is essential in the structural design, because the response of the structure depends on it. A contact problem occurs when not mechanically joined bodies touch each other without becoming rigidly attached. In several cases, high stress concentrations are developed in the contact areas. This fact and the presence of friction and wear often cause crack initiation and fretting fatigue. Examples of this nature are those related with roller and rubber-lined

* Author to whom correspondence should be addressed.

bearings, connecting rods, gear teeth, joint and support elements, soil-structure interaction members, etc. The analytical solution of these problems is not easy and only available for idealised models. For the general study, it is necessary to develop numerical techniques based on the ®nite element method (FEM) or on the boundary element method (BEM), capable of handling arbitrary geometries and including frictional e€ects at the interface, plasticity, possibility of large displacements and large strains, etc. Most of these techniques require, in order to apply the contact conditions between corresponding nodes, dividing each body into elements and nodes, with identical disposition along the surfaces of both bodies, not only at the contact zone, but also at the zone expected to come into contact. First works with this aim are the

0045-7949/98/$19.00 # 1998 Elsevier Science Ltd. All rights reserved. PII: S 0 0 4 5 - 7 9 4 9 ( 9 8 ) 0 0 0 7 4 - 1

316

A. Lorenzana, J. Garrido / Computers and Structures 68 (1998) 315±324

ones presented by Chan and Tuba [1], Francavilla and Zienkievicz [2] and Okamoto and Nakazawa [3] using the FEM, followed by the papers of Anderson [4, 5] and Paris and Garrido [6±8] using the BEM, presenting a general formulation for the two-dimensional case. The requirement of identical discretizations not only takes away freedom when discretizing, but also con®nes its use to the case of the small relative displacement at the interface. This is due to the fact that, in the general case of arbitrary geometry of the bodies and arbitrary elastic properties, placing nodes in advance at the indeformed contact region which occupy the same position in the deformed contact con®guration is not possible, especially if large displacements with large slipping at the interface have taken place. In the present paper this restraint is removed and one node-on-element contact algorithm based on the BEM is set. This algorithm uses a length mapping procedure to place the appropriate point a node is in contact with at the ®nal con®guration of each increment. The solution is obtained through a system of equations of ®xed size, depending on the total number of nodes. Each increment of load is determined by two possible motives: either a node changes some of its boundary conditions, or the geometry has su€ered substantial modi®cations that force to obtain again the new coecients of the matrix involved. Some algorithms have already been used removing the node-on-node requirement, not only by means of the FEM (Bathe and Chaudhary [9] and Kodikara and Moore [10]), but also with the BEM (Olukoko et al. [11] and Paris et al. [12]). Also noteworthy are the BEM±FEM approach to contact problems presented by Liolios et al. [13] and the thorough formulation in the ®eld of contact mechanics that can be found in the work of Antes and Panagiotopoulos [14]. The next sections establish the most relevant aspect of the proposed algorithm and in the last section two examples are presented showing the applicability and reliability of the method.

2. Boundary formulation Before developing a method for the analysis of the interaction of two bodies, the BEM is reviewed brie¯y. It is known that one appropriate formulation for the static analysis of tractions and displacements in a solid body under homogeneous, isotropic and linear elastic behaviour, assuming the hypothesis of small strains and in the absence of body forces, can be done by means of integral equations extended only at the boundary of the body using Somigliana's identity …   Cij …P†ui …P† ˆ U ij …P; Q†pi …Q†ÿT ij …P; Q†ui …Q† dSQ tG

…1† This equation is applied at any instant of time t, being P and Q two points of the current boundary tG and ui and pi the displacement and traction vector at these points. Cij(P) is the free term and U ij (P, Q) and T ij (P,Q) are displacements and stresses associated to the fundamental solution (Kelvin problem [15]). It should be noticed that the integral is extended along the current boundary, that is unknown as long as ui remains unknown. To solve this problem, an updated Lagrangian approach is adopted. Thus, the integration boundary taken is not the current one tG, but some other previously known t-DtG. Eq. (1) provides two subsequent equations associated to any point P belonging to tG. In order to evaluate the boundary integrals and formulate the problem in a discrete form it is necessary to replace the original boundary by the discretized one, using elements and nodes in the usual way and taking the nodes as points P where Eq. (1) is collocated. After performing the numerical integration over the boundary elements, a set of linear algebraic equations is obtained. They can be written, with the usual terminology of the BEM, in the form Hu ˆ Gp

…2†

where the vectors u and p are formed by the components of the displacement and traction vectors at the nodes. The coecients of the matrix G and H, obtained by appropriate Gaussian quadrature, depend

Fig. 1. Element de®nition.

A. Lorenzana, J. Garrido / Computers and Structures 68 (1998) 315±324

on the fundamental solution, the geometry and also on the coordinate system chosen to express each variable ui and pi, as well as on the interpolation procedure selected. In this paper isoparametric quadratic elements (three-noded) have been chosen. The procedure to interpolate any variable ci (representing the components of the displacement vector ui, the traction vector pi or Cartesian coordinates xi) is Ci ˆ

3 X

j N…x† j Ci

;

i ˆ 1; 2

…3†

jˆ1

where x is the intrinsic coordinate (x = [ÿ1., + 1.]) and jci are the values of those variables at the nodes of the elements. The three nodes are placed along the element at 1x = ÿ 1., 2x = k and 3x = + 1., as it is shown in Fig. 1. It should be noticed that the central node is not placed at x = 0.; thus, the nodes of any element can be used with good results (Jacobian of the coordinate transformation almost constant along x) even at advanced stages of load, where the displacements might have been large. The quadratic shape functions jN(x) for each node, showed in Fig. 1, are given by …x ÿ 1†…x ÿ k† ; 2…1 ‡ k† …x ‡ 1†…x ÿ k† 3 N…x† ˆ 2…1 ÿ k† 1 N…x†

ˆ

2 N…x†

ˆ

…x ÿ 1†…x ‡ 1† ; …k ÿ 1†…k ‡ 1† …4†

The determination of the k parameter cannot be direct using the expression   La 1 …5† kˆ2 ÿ Lb 2 because the evaluation of La and Lb, by Gaussian quadrature, needs the shape functions jN(x), that depend on k. Nevertheless, the parameter k can be approached iteratively taking as ®rst approximation for La and Lb the length of the straight segments between nodes. Usually no more than three iterations are needed to determine k with a tolerance under 10ÿ8.

Fig. 2. General de®nition of the problem.

317

In the points of the boundary where a continuous normal vector cannot be accepted it is necessary to place two nodes, one for each element concurring at that point, and deal with them properly [16]. 3. The contact algorithm 3.1. Boundary and contact conditions Fig. 2 presents a general situation of two bodies A and B in contact, bearing some external actions (prescribed tractions and/or restraints in the displacement at the boundaries GA and GB) proportional to a certain parameter P. Let nT k=nCk+nLk be the total number of nodes of body k (k = A,B), with nCk of them placed at the contact interface GC. In the part of the boundary not a€ected by the contact, two variables of the four ones associated to any node (ui, pi, i = 1,2) are the boundary conditions 9 pi ˆ pi in Gkpi = k ˆ A; B; i ˆ 1; 2 …6† ui ˆ ui in Gkui ; so that 2nLA+2nLB unknown variables remain. Gkpi is the part of the boundary of the body k where the traction pi is known and Gkui is that part where the displacement ui is known. The correct partition of the whole boundary Gk in Gkpi and Gkui cannot, in the general case, be de®ned before hand, because it depends on the load applied P. This generates an implicit non-linearity that a€ects the solution technique of the problem. In the nodes at the interface GC between the two bodies, it is known, for the frictionless contact problem under study, only the displacement or the traction in the tangential direction (direction t in Fig. 2). Thus, three of the four variables associated to each node are unknown (3nCA+3nCB). Eq. (2), applied in incremental form to each body, provides 2nTA+2nTB equations that can be expressed in a unique system of equations, namely 2 2 3 32 3 2 3 . . A A A . A . 6 H . 0 7 Du 6 G . 0 76 Dp 7 6 76 7 7 6 …7† 6 . . . . . . . . . 74 . . . . . . 5 ˆ 6 . . .. . .. . . 74 . . . . . . 5 4 4 5 5 B .. B . B Dp Du 0 . H 0 .. GB . . where Du ˆ ‰DuA .. DuB ŠT and Dp ˆ ‰DpA .. DpB ŠT are vectors formed by incremental components of displacements and tractions at the nodes. By means of the contact conditions that will be established next, the number of unknown variables for each node in contact is going to be reduced, enforcing the contact constraints between each contacting node and its corresponding point in the other body. If the discretizations performed on both bodies in the contact region (initially in contact or expected to

318

A. Lorenzana, J. Garrido / Computers and Structures 68 (1998) 315±324

tion and/or displacement does no appear in Eq. (7). To illustrate the procedure, suppose that the traction at the node j of body A ( j < nTA) is expressed as a function of the tractions at nodes s1, s2 and s3 of body B, then, the row i in Eq. (7), that formerly was Fig. 3. Length mapping procedure.

come into contact) are matching so that the contact is just between nodes, the contact constraints can be imposed only a€ecting variables of the system of Eq. (7), being in this case immediate their application. Otherwise, one node (m) of one of the bodies (master body M) will come into contact at a certain point (s) of an element of the other body (slave body S). The contact constraints, enforcing equilibrium and compatibility between node m and point s, are imposed in the local coordinate system by means of  …a† m Du1 ‡ s Du1 ˆ 0: …8† …b† m Dp1 ˆ s Dp1 the subindex 1 a€ecting the normal direction (n in Fig. 2). The point s is determined by its corresponding intrinsic coordinate sx within the element where s belongs to, whose nodes are s1, s2 and s3, being s1 and s3 its extremes. The evaluation of sx is made with the length mapping procedure sketched in Fig. 3 in the form   L2 ÿ Lm x ˆ 1 ÿ 2 …9† s L2 ÿ L1 where Li is the length of the boundary between the reference point Q and the nearest (i = 1) or farthest (i = 2) extreme of the element where point s belongs to, and Lm is the length between point Q and node m. 3.2. The general solution procedure The contact constraints are implicitly included in the system of Eq. (7). Once sx is known, sDu1 and sDp1 could be expressed in function of the corresponding displacement and traction components at the nodes s1, s2 and s3 using Eq. (3). The traction and/or the displacement of each node m in the contact zone can be then related by Eq. (8) with the tractions and/or displacements of the nodes s1, s2 and s3. Thus, this trac-

Fig. 4. Node m coming into contact.

A A A A A HA i1 1 Du1 ‡ H i2 1 Du2 ‡ H i3 2 Du1 ‡


A A A ‡ HA i…2jÿ1† j Du1 ‡ H i…2j† j Du2 ‡


A A A ˆ GA i1 1 Dp1 ‡ Gi2 1 Dp2 ‡


A A A ‡ GA i…2jÿ1† j Dp1 ‡ Gi…2j† j Dp2 ‡




now changes into A A A A A HA i1 1 Du1 ‡ Hi2 1 Du2 ‡ H i3 2 Du1 ‡


A A A ‡ HA i…2jÿ1† j Du1 ‡ H i…2j† j Du2 ‡


A A A ˆ GA i1 1 Dp1 ‡ Gi2 1 Dp2 ‡


! 3 X A B ‡ Gi…2jÿ1† i N…s x† si Dp1

‡

iˆ1 A Gi…2j† j DpA 2 ‡






where jDpA 1 does not appear. The procedure works in the same way for j>nTA as for interpolating displacements. This process must be done for every node at GC belonging either to body A or B. After this, the contact constraints are established. It can be arbitrarily chosen to express the tractions and displacements at the nodes of body k(k = A,B) as functions of those at the nodes of the opposite body. As the number of nodes of the contact zone nCk of both bodies is not necessarily the same, in order to ensure the same number of equations than unknowns, if the displacements of the nodes of body A at GC are expressed as functions of the displacements at the nodes of body B (not necessarily all these nodes belonging to GC), then the tractions of the nodes of body B at GC must be expressed as functions of those of body A. Thus, nCA unknowns on displacements and nCB on tractions come out of the equation system (or vice versa if the choice had been the opposite) and the number of unknowns (2(nLA+nLB) + 3(nCA+nCB) ÿ (nCA+nCB)) agrees with the number of equations (2nTA+2nTB). Once the contact constraints have been applied in that form, it is necessary to apply the prescribed boundary conditions as well as rearrange the system Eq. (7) in the usual way, so that it can be solved. When the contact zone is not known because it depends on the applied load P, an incremental procedure is performed so that the nodes out of the initial contact zone are put into it. In order to ®nd out the fraction of load l that makes a new node (of some of both bodies) to come into contact, the starting point is to calculate the gap g in the reference con®guration (time increment t ÿ Dt) between master node m and its corresponding slave point s by the

A. Lorenzana, J. Garrido / Computers and Structures 68 (1998) 315±324

scalar product tÿDt

g ˆ tÿDt n
tÿDt r

…10†

The point s is determined again by the length mapping procedure (making Lm=Ls) sketched in Fig. 4. t ÿ Dt n is the normal vector at node m, t ÿ Dtr is the Cartesian vector that joins node m and point s, namely tÿDt

r ˆ tÿDts x ÿ tÿDtm x

…11†

and t ÿ Dtix is the position vector of node i at time t ÿ Dt. When Eq. (7) is solved for the load P, the displacements Du obtained allow us to know the provisional deformed con®guration at instant t. In this con®guration, the gap tg, evaluated again as the scalar product between tn and tr, where tr is now given by     t r ˆ tÿDts x ‡ s Du ÿ tÿDtm x ‡ m Du …12† can be positive (not contact yet) or negative (interpenetration). The goal is to determine the load parameter l that makes tg = 0. Despite the fact that every load increment is linear, the determination of l can not be done straight away because sx corresponding to point s is still unknown, as it depends on l (Fig. 5). This drawback is solved scaling the incremental solution found in the following form  t tÿDt m xi ˆ m xi ‡ l m Dui …13† t tÿDt sj xi ˆ sj x0i ‡ l sj Dui and imposing that node m and point s do occupy the same geometric position (tmxi=tsjxi for the gap to be zero), so that, using Eq. (3), the following expression is obtained tÿDt m x0i

‡ lm Dui ˆ

3 X

j N…s x†

ÿtÿDt

sj x0i

‡ l sj Dui




…14†

jˆ1

which becomes the following non-linear system of two

319

equations and two unknowns, sx and l, easily resoluble by any numerical procedure " # tÿDt tÿDt tÿDt s1 xi s2 xi s3 xi 2 ‡ ‡ ‡ sx 2…1 ‡ k† …k ‡ 1†…k ÿ 1† 2…1 ‡ k†  ‡s x2 l

s1 Dui

2…1 ‡ k†

" ‡s x

tÿDt s3 xi

2

ÿ

‡

s2 Dui

…k ‡ 1†…k ÿ 1†

tÿDt s1 xi

2

#

‡

s3 Dui

2…1 ÿ k†

" ‡ zl

s3 Dui

2



ÿ

s1 s1 Dui

2

‡ # ‡

…15†

  ks1 Dui ks3 Dui s2 Dui ÿ ÿ ÿm Dui ‡ ‡l 2…1‡k† …k ‡ 1†…kÿ1† 2…1ÿk† "

# tÿDt ktÿDts1 xi ktÿDts3 xi tÿDt s2 xi ‡ ÿ ÿ ÿ m xi ˆ 0: 2…1‡k† …k‡1†…kÿ1† 2…1ÿk† The procedure chosen has been Newton±Raphson, which starts with a ®rst estimation for l and sx calculated expecting a linear behaviour, given by   tÿDt L2 ÿ Lm g l ˆ tÿDt sx ˆ 1 ÿ 2 L2 ÿ L1 g ÿt g Three or four iterations are usually enough to achieve the solution with a tolerance under 10ÿ8. If sx is close to the limits of its domain [ÿ1., + 1.], this procedure might produce values out of the valid domain. This indicates that the initial slave element, to which sx refers, is not the correct one. Then, sx must be found in the adjacent slave element and the iterative procedure described started anew. Once l has been found out, the solution Du and Dp is scaled and the incremental result obtained is accumulated to that known at instant t ÿ Dt, namely 9 t u ˆ tÿDt u ‡ lDu = t …16† p ˆ tÿDt p ‡ lDp ; t x ˆ 0 x ‡t u where 0x is the position vector of all nodes in the initial con®guration. Interpenetrations at the nodes between the bodies are prevented by the construction of the algorithm itself, but it could be necessary to check if any traction (not compression) appears after the application of each increment of load at the nodes of the contact zone, in order to avoid that incompatibility [Eq. (8b)] by leaving the node outside the contact zone.

4. Applications

Fig. 5. Procedure to ®nd out the load parameter.

Two examples have been selected to show the applicability of the proposed method for the analysis of contact problems under in®nitesimal and ®nite displa-

320

A. Lorenzana, J. Garrido / Computers and Structures 68 (1998) 315±324

Fig. 6. De®nition of the geometry, properties and boundary conditions.

cement. The problems deal with compression of a rectangular punch and a roller on a foundation and are analysed under plane strain condition. 4.1. Compression of a rectangular punch on a foundation Fig. 6 represents the geometry, the boundary conditions and the properties of both bodies, being Ek the modulus of elasticity of body k and nk the Poisson's ratio. In this kind of contact problem, when displacements are in®nitesimal, the length of the contact zone does not depend on the amount of load applied. Therefore, the problem is linear and only one increment is required for its solution. A detailed study of this problem can be found in Paris and Garrido [6]. However, when displacements are large so that the solution changes slightly from the linear one, as it is going to be shown, the problem is no longer linear. Due to the symmetry of the problem, only half of the geometry needs to be modelled and an implicit symmetry procedure is used [7] so that no boundary elements are set in the symmetry axis. Twenty elements have been employed for body A and 21 for body B. In order to take into consideration that the contact stress at point D is singular, the distribution of the nodes is thicker as that point is approached. An increment of load is made each time the displacement of any node exceeds 0.25% of the greatest length of the body it belongs, being this percentage established beforehand. It has been veri®ed that any percentage lower than that does not a€ect the solution substantially, although a larger number of increments are necessary. This reveals that convergence is achieved, that is, the geometry is updated as many times as necessary so that

the hypothesis on which the updated Lagrangian approach is based is ensured. In the context of the large displacement theory exposed, the geometry should be updated and the coef®cients of the matrices H and G computed at each increment. This would require an important amount of computational time, comparable with that necessary for the resolution of the system of equations. From a numerical point of view, in order to reduce computational time, generating those coecients every increment could be avoided, as they are not expected to vary a great deal. To show the in¯uence of the latter, the solution using both possibilities (updating and not updating) is presented. The distribution of the contact stresses along the contact zone for four di€erent loads is shown in Fig. 7 (corresponding to q1=175, q2=300, q3=700 and q4=1400. N/mm2). It is important to notice that the tractions are plotted in their corresponding deformed con®guration. Thus, although the original length of the contact zone is 50 mm, when the load applied is 1400 N/mm2, the length in this deformed con®guration is almost 60 mm. At the ®rst stages of the load application, when displacements are small, the solution agrees with the linear one obtained by Garrido [17]. From a certain load (say about q1) onwards, the di€erences between the linear solution of reference and the solution obtained by the proposed method become signi®cant. The one obtained by the proposed method presents higher stresses around the point D (just in the point D and under the assumption of elastic behaviour, the stress should be in®nite), therefore the stresses in the central area are slightly lower. It is interesting to observe that the traction solution obtained without updating the coecients of G and H is almost the

A. Lorenzana, J. Garrido / Computers and Structures 68 (1998) 315±324

321

Fig. 7. Normal stress along the contact zone.

same as that of reference. This e€ect can be explained bearing in mind that only the geometry is updated and it does not su€er any qualitative changes (relating the relative position of nodes and their normal vectors) that may a€ect the incremental process in comparison with the linear situation. Fig. 8 shows the modulus of the displacement of point A (Au) and point D (Du) versus the applied load. Again the solutions obtained di€er from the linear ones when the bodies su€er no in®nitesimal displacements. For the same load applied, when the large displacement theory is considered, the displacements of both points are smaller than those evaluated assuming linear theory. If only the geometry is updated, not G and H, the displacement solution obtained di€ers from the linear one, the former being between the linear and the fully updated. 4.2. Compression of a roller on a foundation Fig. 9 represents the geometry, the boundary conditions and the properties of both bodies. The contact zone depends on the load, which is applied with a pre-

scribed displacement u at the base of the foundation. This problem, when the displacements are suciently small, has an analytical solution because it is a Hertzian problem. Otherwise, a non-Hertzian problem, without analytical solution, appears. Because of the symmetry, only one quarter of the geometry is employed for the discretization and again implicit symmetry is used. 20 elements have been employed for body A and 26 for body B. The prescribed displacement is progressively applied until the deformed con®guration at the contact region is the one presented in Fig. 10. An increment of load is made every time a node comes onto the contact zone. The geometry is updated after the application of each increment of load, as well as the computation of the coecients of H and G. No additional updatings are required, the ®ne discretization used ensures that the geometry is updated as often as necessary. Again, in order to reduce computational time and to observe the in¯uence of not updating G and H, the solution with this possibility is also presented.

Fig. 8. Displacements versus total load applied.

322

A. Lorenzana, J. Garrido / Computers and Structures 68 (1998) 315±324

Fig. 9. De®nition of the geometry, properties and boundary conditions.

It can be observed (Fig. 10) that when H and G are not updated, and those obtained in the initial con®guration are employed, the deformed con®guration (dotted line) is slightly di€erent, with the contact zone slightly larger for the same applied load. This e€ect is also observed in Fig. 11, which shows the length of contact versus the applied load. This load has been evaluated by integration of the resultant stress distribution at the contact zone of either of the bodies. Fig. 12 shows the stress distribution for ®ve di€erent lengths of contact. When this magnitude is under, let us say 10.0 mm, the solution almost coincides with the analytical one by Hertz. Otherwise, the di€erences are gradually bigger as the displacements grow. If the Hertz solution were applicable (that is not the case because the contact is not Hertzian any more) for the

same length of contact, the predicted applied load and also the stress distribution would be smaller than those calculated without updating G and H. Moreover, the one calculated updating G and H is smaller than the former. 5. Conclusions An incremental procedure for two-dimensional frictionless contact problems, under static and proportional loading conditions, with the possibility of large displacements has been presented. The node-onelement algorithm of contact has been implemented using the boundary element techniques, probably the most appropriate for a contact problem, where the

Fig. 10. Con®gurations at the contact zone.

A. Lorenzana, J. Garrido / Computers and Structures 68 (1998) 315±324

323

Fig. 11. Applied load along the contact zone.

Fig. 12. Contact stresses along the contact zone.

magnitudes of interest are on the boundary and where domain loads usually do not appear. The method has shown its reliability when the displacements are small, getting numerical results in good agreement with the analytical ones, if possible, and with those obtained by other authors. Three-noded isoparametric elements are employed in the discretization. With these types of element, even when the displacements at the contact zone are large, the algorithm guarantees compatibility (no interpenetrations) at the nodes and equilibrium (same distribution of contact stress) along the contact zone. Nevertheless some problems of ill conditioning of the system of equations might appear, that can be avoided changing slightly the discretization performed. The method has a very wide ®eld of application, but a future extension to some aspects like dealing not only with advancing contact but also receding contact, frictional e€ects at the contact zone, inclusion of large

strain theory, possibility of plastic zones, three-dimensional contact problems, etc. is necessary. Through the examples treated, it is observed that even for the classic contact problems there are signi®cant di€erences whether they are studied under the ®nite displacement theory or the in®nitesimal. This justi®es the development of contact techniques such as the one shown in this paper in order to achieve the right approach to the problem.

References [1] Chan SK, Tuba IS. A ®nite element method for contact problems of solid bodiesÐI. Theory and Validation. Int J Mech Sci 1971;13:615±25. [2] Francavilla A, Zienkiewicz OC. A note on numerical computation of elastic contact problems. Int J Numer Meth Engng 1975;17:913±24.

324

A. Lorenzana, J. Garrido / Computers and Structures 68 (1998) 315±324

[3] Okamoto N, Nakazawa M. Finite element incremental contact analysis with various frictional conditions. Int J Numer Meth Engng 1979;14:337±57. [4] Anderson T. The Boundary element method applied to two-dimensional contact problems with friction. In: Boundary element methods. Springer, Berlin, 1982:239± 58. [5] Anderson T, Fredriksson B, Allan Person BG. The boundary element method applied to two-dimensional contact problems. In: Brebbia CA, editor. New developments in boundary element methods. CML Publications, 1980:247±63. [6] Paris F, Garrido JA. An incremental procedure for friction contact problems with the boundary element method. Engng Anal Boundary Elements 1989;6:202±13. [7] Paris F, Garrido JA. On the use of discontinuous elements in two dimensional contact problems. In: Brebbia, Maier, editors. Boundary Elements VII. Springer, Berlin, 1985;Vol. II:1327±30. [8] Garrido JA, Foces A, Paris F. B.E.M. applied to receding contact problems with friction. Math Comput Modell 1991;15:143±54. [9] Bathe KJ, Chaudhary A. A solution method for planar and axisymmetric contact problems. Int J Numer Meth Engng 1985;21:65±88.

[10] Kodikara JK, Moore ID. A general interaction analysis for large deformations. Int J Numer Meth Engng 1993;36:2863±76. [11] Olukoko OA, Becker AA, Fenner RT. A new boundary element approach for contact problems with friction. Int J Numer Meth Engng 1993;36:2625±42. [12] Paris F, Blazquez A, CanÄas J. Contact problems with nonconforming discretizations using boundary element method. Comput Struct 1995;57:829±39. [13] Liolios A, Tsalkatidis S, Logothetidis S. A BEM±FEM approach to unilateral contact problem of dynamic soil± pile interaction. In: Brebbia CA et al., editors. Boundary elements XVII. Computational Mechanics Publications, Southampton, 1995:387±92. [14] Antes H, Panagiotopoulos PD. The boundary integral approach to static and dynamic contact problems. Equality and inequality methods. Birkhauser, Basel, 1992. [15] Brebbia CA. The boundary element method for engineers. Pentech Press, London, 1978. [16] Banerjee PK, Butter®eld R. Boundary element method in engineering science. McGraw-Hill, London, 1981. [17] Garrido JA. The contact problem in elasticity by using integral equation methods. PhD Thesis. Univ. of Las Palmas, Spain, 1986.