Influence of a curvature discontinuity on the pressure distribution in two-dimensional frictionless contact problems

International Journal of Engineering Science 38 (2000) 197±205

www.elsevier.com/locate/ijengsci

In¯uence of a curvature discontinuity on the pressure distribution in two-dimensional frictionless contact problems C. Mechain-Renaud *, A. Cimetiere Universit e de Poitiers et E.N.S.M.A., L.3M.A., S.P.2M.I., B.P. 179, 86960 Futuroscope Cedex, France Received 23 April 1998; received in revised form 4 August 1998; accepted 16 December 1998

Abstract The purpose of this paper is to discuss the in¯uence of a curvature discontinuity in two-dimensional Hertz contact problems. We analyse here the frictionless unilateral contact problem between an elastic cylinder and a rigid plane support. At the ®rst contact point, the surface pro®les present a curvature discontinuity. The problem is solved in the context of the boundary element method by a new numerical method based on a regularity result for the contact area boundary. The numerical results show that the curvature discontinuity strongly in¯uences the contact pressure distribution which signi®cantly di€ers from the classical Hertz theory one. Ó 1999 Elsevier Science Ltd. All rights reserved.

1. Introduction The study concerns two-dimensional Signorini contact problems between an elastic cylinder and a rigid plane foundation. In usual frictionless Signorini problems, the normal contact pressure distribution is given by the well-known Hertz theory. One of its main assumptions is to consider that the boundaries of the bodies coming into contact are locally regular (C2 regularity: the pro®les of the surface are continuous up to their second derivatives); in particular in the twodimensional case, the curvature of each pro®le is assumed continuous in the vicinity of the ®rst contact point. The aim of the paper is to show that the Hertzian contact pressure distribution is radically modi®ed in the case of two-dimensional frictionless contact problems when a curvature discontinuity exists at the ®rst point of contact. The study is carried out in the context of the usual linear elasticity, with plane strains assumption.

*

Corresponding author. Tel.: +33-5-4949-6796; fax: +33-5-4949-6791. E-mail: [email protected] (C. Mechain-Renaud)

0020-7225/00/$ ± see front matter Ó 1999 Elsevier Science Ltd. All rights reserved. PII: S 0 0 2 0 - 7 2 2 5 ( 9 9 ) 0 0 0 1 8 - X

198

C. Mechain-Renaud, A. Cimetiere / International Journal of Engineering Science 38 (2000) 197±205

Such problems are non-linear. The numerical resolution of unilateral frictionless contact problems is usually achieved through one of the following approaches. The ®rst one concerns iterative methods. The structure is loaded by successive increments and the unilateral contact condition is veri®ed at each step. This method was used by Chan and Tuba [1], Francavilla and Zienkiewicz [2] and Gaertner [3] who developed modi®ed ®nite element methods, working on sti€ness matrix analysis. Andersson et al. [4] applied the boundary element method to two-dimensional problems, using a ``trial and error'' method which is also an iterative one. Namely, the structure is loaded at one time and the contact area is found iteratively by testing nodes, according to the normal surface traction signs obtained at the previous step. Then, the new contact zone being known, normal surface tractions are recalculated and so on, until the solution is reached. A second approach also used in unilateral frictionless contact problems is based on mathematical programming techniques. Displacements have to satisfy the unilateral contact conditions. The problem is equivalent to the minimization of the potential energy under unilateral constraints on the displacements (see Ref. [5]) and so can be regarded as a mathematical programming problem. Fremond [6] studied the Signorini problem under a variational form and reduced it to researching a saddle-point of a Lagrangian. Panagiotopoulos [7] generalized this approach to the inelastic foundation and presented some dual forms of the variational inequalities governing the contact problems. Meanwhile, Conry and Seireg [8] developed a linear programming procedure using a simplex-type algorithm in order to solve Hertzian problems. In a more original way, Kalker and Van Randen [9] preferred to use the enthalpy minimum principle and formulated the problem as a boundary element problem. They then solved it by a standard quadratic programming method. An appropriate linearization of the contact condition was adopted by Hung and De Saxce [10] which enabled to reduce the problem to a quadratic programming one. The works of Barbosa and Feijoo [11] compared the Lemke algorithm solving the dual problem with the Uzawa algorithm used to solve saddle-point problems. The ®nite element method is adapted to contact problems through a Lagrangian formulation in which the constrained problem is transformed into an unconstrained one by introducing additional Lagrange multipliers. See for instance Simo et al. [12] who developed a perturbed Lagrangian procedure including a penalty term. For incompressible materials contact problems, Glowinski et al. [13] introduced a di€erent formulation in which the contact constraints are taken into account by a penalty method and the penalized problem is turned into a saddle-point form in order to take into account the incompressibility condition. They used an augmented Lagrangian formulation of the penalized problem. On the other hand, penalty methods transform the constrained problem into an unconstrained one without additional variables and the constraint is approximately satis®ed for ®nite values of the penalty parameter. That is what Kikuchi and Song [14] did when they discretized the variational formulation of the penalized problem by ®nite element methods. The resolution of the model problem considered in the paper will be performed by a new numerical method based on the strong formulation of the problem which has already successfully worked for several classical Hertz problems [15] as well as for the obstacle problem of an elastic membrane coming into contact with a rigid parallel plane [16]. The continuous problem is discretized by the boundary element method.

C. Mechain-Renaud, A. Cimetiere / International Journal of Engineering Science 38 (2000) 197±205

199

2. Resolution method In this new approach, the position of the contact area boundary appears as the main unknown of the problem. The method is based on the following fundamental remark: both the gap between the two solids and the normal contact stress are simultaneously null when the boundaries of the two bodies are regular in the vicinity of the contact zone. The problems considered here are two dimensional; the contact area is supposed to be an interval, whose limits M1 and M2 are, respectively, located by q1 and q2 (Fig. 1). The method can be extended to problems with a multiply connected contact area. We note d the gap between the elastic body and the rigid support and rn the normal contact stress, decomposed in the form r
n ˆ rn n ‡ rt :

…1†

The unilateral contact conditions are d
rn ˆ 0;

d P 0; rn 6 0:

…2†

At both extremities M1 and M2 , the following conditions hold dˆ0

and rn ˆ 0:

…3†

In the problem, q1 and q2 are two additional unknowns but Eq. (3) gives two complementary equations rn …M1 † ˆ 0;

…4†

rn …M2 † ˆ 0:

…5†

The boundary element method allows to take easily into account both prescribed displacements and forces conditions (3) which have to be satis®ed simultaneously on the contact zone boundary. We denote by X the domain which represents the cross section of the linear elastic in®nitely long

Fig. 1. Contact zone description.

200

C. Mechain-Renaud, A. Cimetiere / International Journal of Engineering Science 38 (2000) 197±205

asymmetrical cylinder. Its equilibrium is given by the following two-dimensional integral equations written for any point x on the boundary oX 1 ui …x† ˆ 2

Z

uij …x; y†pj …y† dSy ÿ

oX

Z

pij …x; y†uj …y† dSy ‡

oX

Z

uij …x; y†fj …y† dV ;

…6†

X

where pi ˆ rij nj , uij are the fundamental displacements and pij are the fundamental stresses associated with the elementary Kelvin solution for two-dimensional elastostatic problems [17]. In the following, the body forces f are assumed small enough to be neglected. The boundary C of the domain X is divided into M segments C ˆ C1 [


[ C M :

…7†

On each segment Cm , displacements and stress vectors are approximated by constant values taken at the middle point am uj …y†  uj …am † ˆ ujm ; pj …y†  pj …am † ˆ pjm ;

y 2 Cm ;

…8†

y 2 Cm :

After discretization, the boundary integral Eq. (6) reads M X 2 M X 2 X X Aijmk pjm B0ijmk ujm ˆ 0; mˆ1 jˆ1

i ˆ 1; 2; k ˆ 1 . . . M:

…9†

mˆ1 jˆ1

Aijmk and B0ijmk result in the integration of uij …ak ; y† and pij …ak ; y† on each boundary element [15]. They depend on q1 and q2 in a non-linear way through the discretization points. Consequently, the boundary integral equations lead to a non-linear system of 2M equations with 2M unknowns at which we have to add the two unknowns q1 and q2 and the two supplementary Eqs. (4) and (5). This non-linear system is solved by the standard Newton±Raphson method. 3. The test problem In this problem, we consider the contact between a very long cylinder and a ®xed rigid plane support. The cross section is composed of a left circular part and of a parabolic part on the other side (Fig. 2). The loading is imposed through a vertical displacement UD along the upper plane part C4 . So any rigid displacement is impossible. The parts C1 and C3 of the boundary oX are not loaded. Thus, p ˆ r
n ˆ 0. The contact interval C2 is limited by M1 and M2 . In this interval, the vertical displacement is known. The solid pro®les are de®ned by: q y1 …x† ˆ ÿ R21 ÿ x2 ;

…10†

C. Mechain-Renaud, A. Cimetiere / International Journal of Engineering Science 38 (2000) 197±205

201

Fig. 2. The solid diagram.

y2 …x† ˆ

1 2 x ÿ R1 ; 2R1 a

…11†

where y1 de®nes the circular part of radius R1 and y2 the parabolic part. So both pro®les join up with each other in (0,ÿR1 ) with a horizontal tangent. At this ®rst point of contact, the curvature radius of the parabolic part is R2 ˆ aR1 . The parameter a measures the curvature discontinuity between the two pro®les at the joining point. Now di€erent examples are presented in the case of a weak sti€ness: R1 ˆ 100 mm, E ˆ 4000 MPa and m ˆ 0:3. We obtain similar results for greater sti€nesses. The imposed vertical displacement UD is equal to 0.5 mm. In each case, q1 is initialized to R1 =2, q2 to a R1 =2 and the other unknowns to zero. The calculus process gives the solution with a maximum system size of 300. 4. Results Several results are presented here for di€erent values of a showing the in¯uence of the pro®le curvature discontinuity on the contact stress distribution. When the curvatures of both pro®les are the same …a ˆ 1† we note (Fig. 3) that the numerical contact pressure agrees with the Hertz theory q

one [18], i.e.: p…x† ˆ p0 1 ÿ …x=q†2 and that the contact zone is symmetric although the structure is not perfectly symmetric. For the following comparisons, we will take the Hertz pressure as reference. We are going to show that the presence of a curvature discontinuity signi®cantly modi®es the Hertz contact pressure distribution. At ®rst, we consider two cases with a < 1 (Fig. 3). The contact zone side which is parabolic is discretized into 15 linear elements whereas the circular part is divided into 30 elements, in order to take into account the di€erence in size of the two pro®les. The number of elements on each boundary part is ®xed at the beginning of the calculus. The mesh changes with q1 and q2 values at each iteration but, in the example, the

202

C. Mechain-Renaud, A. Cimetiere / International Journal of Engineering Science 38 (2000) 197±205

Fig. 3. Contact pressure distribution for a < 1.

contact zone is discretized into 45 elements up to the end of the calculus. For R2 smaller than R1 , the contact area is wider for the circular side and the contact pressure maximum moves towards the parabolic side. A pressure peak appears, far from the reference pressure, and the result becomes all the more remarkable as a goes away from 1. Moreover as shown in Table 1 six iterations only are needed to reach the solution. Table 1 gives the evolution of the contact zone during the successive iterations in the case of a very asymmetrical situation: a ˆ 1=8. The calculus was performed on an Origin 200 server within 48 s CPU time. As shown in Fig. 4, the inverse phenomenon happens when the curvature radius R2 is greater than R1 (for a > 1). Indeed for the two cases considered, we can note that the contact pressure maximum moves towards the circular part and remains less than the reference pressure. These results highlight the point that a curvature discontinuity leads to a pressure peak on the contact area unpredicted by the Hertz theory. Contrary to what could be expected, the maximum of the Table 1 Contact zone evolution according to the successive iterations Number of iterations

q1 (10ÿ3 m)

q2 (10ÿ3 m)

ÿ 1 2 3 4 5 6

50.00 25.188 13.360 7.942 5.893 5.446 5.420

6.250 3.376 1.990 1.449 1.386 1.386 1.388

C. Mechain-Renaud, A. Cimetiere / International Journal of Engineering Science 38 (2000) 197±205

203

Fig. 4. Contact pressure distribution for a > 1.

pressure does not coincide with the ®rst point in contact. It is always located on the greatest curvature side (Fig. 5). This seems paradoxical but can be justi®ed through the analysis of the normal surface forces distribution corresponding to the imposed displacement on the upper part

Fig. 5. Contact pressure distributions.

204

C. Mechain-Renaud, A. Cimetiere / International Journal of Engineering Science 38 (2000) 197±205

Fig. 6. Normal surface forces.

C4 (Fig. 6). This distribution must prevent the solid from rolling on the smaller curvature side and hence is necessarily greater on the other side. Thus for a > 1 the maximum of the surface forces moves towards the circular part whereas for a < 1 it moves towards the parabolic part. We can notice that for a ˆ 1, the distribution of the normal surface forces is practically symmetrical and its maximum is just vertically above the initial point in contact.

Fig. 7. Circular pro®le with an average curvature.

C. Mechain-Renaud, A. Cimetiere / International Journal of Engineering Science 38 (2000) 197±205

205

5. Conclusion In conclusion we note that the Hertz theory gives the contact pressure distribution only when the surface pro®les are C2 in the vicinity of the contact area. The study of the two-dimensional asymmetric problem with imposed displacement considered in the paper has shown that a curvature discontinuity of the pro®le generates a pressure distribution far from the Hertz theory one. A pressure peak appears on the side whose curvature radius is the smaller. As shown in Fig. 7 this peak cannot be described by considering a circular pro®le with an average curvature. Moreover if we consider a whole parabolic body whose curvature is the smaller of the two previous ones, we take into account the peak but we make an important error on the contact area extension. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18]

S.K. Chan, I.S. Tuba, Int. J. Mech. Sci. 13 (1971) 615. A. Francavilla, O.C. Zienkiewicz, Int. J. Num. Meth. Engng. 9 (1975) 913. R. Gaertner, Comput. & Structures 7 (1977) 59. T. Andersson, B. Fredriksson, B.G.A. Persson, New developments in Boundary Element Methods, C.A. Brebbia, 1980, p. 247. G. Duvaut, J.L. Lions, Les inequations en mecanique et en physique, Dunod, Paris, 1972. M. Fremond, Etude des structures visco-elastiques strati®ees soumises  a des charges harmoniques et de solides elastiques reposant sur des structures, These, Universite de Paris, 1971. P.D. Panagiotopoulos, Ingenieur-Archiv 44 (1975) 421. T.F. Conry, A. Seireg, J. Appl. Mech. 38 (1971) 387. J.J. Kalker, Y. Van Randen, J. Eng. Math. 6 (1972) 193. D.H. Hung, G. De Saxce, Comput. & Structures 11 (1980) 55. H.J.C. Barbosa, R.A. Feij oo, Variational Methods in Engineering, C.A. Brebbia, 1985, pp. 6±53. J.C. Simo, P. Wriggers, R.L. Taylor, Comput. Meth. Appl. Mech. Engng. 50 (1985) 163. R. Glowinski, M. Vidrascu, P. Le Tallec, Finite element methods for nonlinear problems, in: Europe ± US Symposium, Tronheim, Norway, 1985, 745. N. Kikuchi, Y.J. Song, Quart. Appl. Math. 39 (1981) 1. A. Cimetiere, C. Mechain, C. R. Acad. Sci. 323 (1996) 455. A. Cimetiere, T., Texier, in: Second International Symposium of Contact Mechanics, Raous-Jean-Moreau, Carry Le Rouet, France, 1994. C.A. Brebbia, J.C.F. Telles, L.C. Wrobel, Boundary Element Techniques, Springer, Berlin, 1984. K.L. Johnson, Contact Mechanics, Cambridge University Press, Cambridge, 1985.