Computational Mechanics 31 (2003) 390–399 Ó Springer-Verlag 2003 DOI 10.1007/s00466-003-0441-4

BEM and FEM analysis of Signorini contact problems with friction C. Renaud, Z.-Q. Feng

390 Abstract In this paper, we present a comparative study of the boundary element method (BEM) and the finite element method (FEM) for analysis of Signorini contact problems in elastostatics with Coulomb’s friction law. Particularities of each method and comparison with the penalty method are discussed. Numerical examples are included to demonstrate the present formulations and to highlight its performance. Keywords Contact, Friction, Finite element method, Boundary element method

1 Introduction The analysis of contact problems with friction is of great importance in many engineering applications. The numerical treatment of the unilateral contact with dry friction is certainly one of the non smooth mechanics topics for which many efforts have been made in the past. In the literature, many attempts have been developed to deal with such problems using the finite element method (FEM), these include the penalty function method [1–4], the flexibility method [5, 6], the mathematical programming method [7–9], the Lagrangian multiplier method [10, 11] and the augmented Lagrangian method [12–17]. A large literature base is available for a variety of numerical algorithms (See the monographs of Zhong [18] and Wriggers [19]). However, the modeling of frictional contact problems by the boundary element method (BEM) remains limited. At first a trial and error algorithm has been developed [20]. Then numerical procedures based on an incremental approach [21, 22], on the flexibility method [23] or on mathematical programming [24] have been proposed. For elastoplastic materials in the presence of frictionless contact conditions, a formulation characterized by two variational principles based on boundary integral equation method has been presented [25]. The aim of the present paper is to present a comparative study of BEM and FEM for analysis of Signorini contact problems in two-dimensional elastostatics with Coulomb’s Received: 17 December 2002 / Accepted: 18 March 2003

C. Renaud, Z.-Q. Feng (&) Centre d’E´tudes de Me´canique d’Iˆle de France, Laboratoire de Me´canique d’Evry (CEMIF-LME), Universite´ d’E´vry-Val d’Essonne, 40, rue du Pelvoux, 91020 E´vry cedex, France e-mail: [email protected]

friction law. After a definition of contact kinematics, we outline the penalty method for comparison reason. We present then the BEM and FEM developed for contact modeling. Two numerical examples are performed in this study to show the validity of the developed models. The performance of present methods is reported as compared to the general purpose finite element code ANSYS in which the penalty method is used for contact modeling.

2 Contact kinematics In this section, the geometric and kinematic quantities found suitable for describing the contact compatibility of deformable bodies are defined. First of all, some basic definitions and notations are set up. For the sake of simplicity, let us consider contact between two bodies X1 and X2 , one of which may be a rigid foundation. In order to state the contact constraints, we have to find the minimum distance of a point P of one body with respect to the other one. The displacements of the particles of X1 and X2 being respectively u1 and u2 , the relative displacement is: u ¼ u1
u2 . Let r be the contact traction acting at P from X2 onto X1 . Then X2 is subjected to the traction
r, acting from X1 . Let n denote the normal unit vector at the projection point P0 to the bodies, directed towards X1 , and Tðt1 ; t2 Þ denotes the orthogonal plane to n in R3 (Fig. 1). Any element u and r may uniquely be decomposed in the form: u ¼ ut þ un n; ut 2 T; un 2 R r ¼ rt þ rn n; rt 2 T; rn 2 R

ð1Þ

Classically, a unilateral contact law is characterized by a geometric condition of non-penetration, a static condition of no-adhesion and a mechanical complementarity condition. These three conditions are so-called Signorini conditions written in terms of the signed contact distance xn and the normal contact force rn :

xn  0;

rn  0 and

rn xn ¼ 0

ð2Þ

where xn denotes the magnitude of the gap between the contact node and the target surface and is a violation of the contact compatibility. With the approximation of small displacement, xn is given by

xn ¼ g þ un

ð3Þ

with the initial gap:

g ¼ ðx1
x2 Þ  n

ð4Þ

Fig. 3. Coulomb’s friction law Fig. 1. Projection and gap vector

method or a combined penalty plus Lagrange multiplier method. The penalty method approximately enforces compatibility by means of a contact stiffness Kn (i.e., the penalty factor). The normal contact force is calculated by




rn ¼

Kn xn 0

if xn  0 if xn > 0

ð8Þ

For the combined method, the Lagrange multiplier component of force is computed iteratively. It is expressed as

Fig. 2. Unilateral contact law

Classically, a rate independent dry friction law is charac- rn ¼ minð0; Kn xn þ kiþ1 Þ terized by a kinematic slip rule. Let Kl denote Coulomb’s where
cone:

Kl ¼ fr 2 R3 such that jrt j  lrn g

ð5Þ

where T is equipped with the Euclidean norm j  j. l is the friction coefficient and Kl is a closed convex set. The complete contact law is a complex non smooth dissipative law including three statuses: no contact, contact with sticking and contact with sliding. The resulting analytical transcripts yields two overlapped ‘‘if. . .then. . .else’’ statements:

if rn ¼ 0 then xn  0 ! no contact else if r 2 Kl then u ¼ 0 ! sticking else ðrn > 0 and r 2 Kl Þ;
 rt ! sliding xn ¼ 0 and 9k > 0 such that jut j ¼
k jrt j ð6Þ An alternative statement is the inverse law:

if xn > 0 then r ¼ 0 ! no contact else if u ¼ 0 then r 2 Kl ! sticking ut else u 2 T; rn > 0 and rt ¼
lrn jut j

ð7Þ ! sliding

For plane problems, the contact law can be represented through two diagrams (Figs. 2 and 3). Nevertheless, take care of the fact that the diagram of Fig. 3 is misleading, because, it should be considered as a family of multivalued mappings controlled by a scalar parameter rn .

kiþ1 ¼

ki þ aKn xn ki

ð9Þ

if jxn j  e if jxn j < e

ð10Þ

a is an internally computed factor ða < 1Þ and e is a userdefined compatibility tolerance. It is very known that the choice of the penalty factor is a delicate task for the user. Low factor leads to important contact penetration and high factor results in numerical stability problems and high numerical cost. The calculation of tangential contact forces is even more difficult. It has been shown [26] that the stiffness matrix for contact is formed by the outer product of the interpolation vectors Nn ; Nt and the penalty factor Kn ; Kt :

8 T T > : 0

if stickingcontact if sliding or frictionlesscontact if open contact ð11Þ

This means that, in the iterative solution procedure, the global stiffness matrix will be modified by introduction of contact elements. The reader is invited to refer to the theory manual of ANSYS [27] for further information.

4 The boundary element method (BEM)

4.1 Principle For the comparative study we use a BEM resolution based on a recent approach which consists in determining contact zone sizes (slip and stick areas) at the same time as the 3 usual unknowns (displacements and surfaces forces in Penalty algorithms BEM). This method has already been tested on some exAs a comparison, we outline here the method used in the general purpose finite element program ANSYS for amples like the obstacle problem of an elastic membrane analysis of contact problems. Two methods of satisfying coming into contact with a rigid parallel plane [28], Signorini frictionless hertzian problems [29, 30] and other contact compatibility are used in ANSYS: a penalty

391




jrt j 
lrn ut ¼ 0

whereas on the sliding subareas, we have:
jrt j ¼
lrn ut ¼
krt ðk  0Þ

Fig. 4. Contact zones division example


classical frictional problems [31]. We recall the principle of the resolution method. The position of the contact zones is defined through geometrical parameters. Introduction of these extra unknowns requires additional equations in the same number. Under the assumption of regular bodies boundaries and loads in the vicinity of the contact zone, additional equations are given by limit contact and limit slip conditions. This allows to reduce the unilateral contact problems to a nonlinear equations system.

4.2 Limit contact and limit slip conditions To define these conditions we consider the two dimensional example of a frictional contact between an elastic body and a rigid plane foundation as described in Fig. 4. Let us suppose that the contact zone is divided into an internal adhesion area surrounded by two sliding zones. M1 and M2 are the extremities of the contact zone and I1 ; I2 the stick area limits. Here the contact conditions are given writing that on the boundary between the contact zone and the free of load edge, the gap between the two bodies and the contact pressure are simultaneously zero [28, 29]. Considering that M1 and M2 belong to the contact zone, the conditions between the initial gap and the normal displacements are automatically satisfied. The limit contact conditions are then given by: rn ðM2 Þ ¼ 0

ð12Þ

where rn is the normal stress. Thanks to overabundant conditions (12), the parameters a and b can be introduced as additional unknowns. In the same way, additional equations at I1 and I2 points are written in order to integrate the parameters c and d into the problem unknowns. As mentioned in Sect. 2, the Coulomb friction law can be defined using contact reaction forces (Eq. 6). This definition is usually adopted by the finite element method. However, in the case of the boundary element method, the friction law is usually defined in terms of contact stresses:

8 < jrt j  ljrn j jr j < ljrn j ) ut ¼ 0 : t jrt j ¼ ljrn j ) 9k  0=ut ¼
krt

ð15Þ

Under assumptions of regularity in the vicinity of the contact zone, two conditions have to be simultaneously satisfied on the boundaries between adhesion zone and sliding ones:

392

rn ðM1 Þ ¼ 0;

ð14Þ

ut ¼ 0 jrt j ¼
lrn

ð16Þ

Considering that the boundary points I1 and I2 belong to the sliding zones, the equation involving normal and tangential stresses is automatically verified. The additional equations named limit sliding conditions are then given by:

ut ðI1 Þ ¼ 0;

ut ðI2 Þ ¼ 0

ð17Þ

Using the limit contact conditions (12) and the limit sliding ones (17), the problem formulation can take into account the contact parameters a; b; c and d. More generally, if contact zones are defined by P contact parameters qq (q ¼ 1 to P) and adhesion areas
r (r ¼ 1 to Pa ), then overby Pa contact parameters q abundant conditions will be given by the contact limit conditions:

rn ðMq Þ ¼ 0 q ¼ 1 to P

ð18Þ

(where Mq are boundary points between contact zones and no contact ones) and by the limit sliding conditions:

ut ðIr Þ ¼ 0 r ¼ 1 to Pa

ð19Þ

(where Ir are the boundary points between sliding zones and adhesion ones). So the contact and sliding areas positions are defined by P0 ð¼ P þ Pa Þ contact parameters. Limit contact conditions (18) and limit sliding ones give the P0 lacking equations.

4.3 Contact problem solution The usual form of the Coulomb law (13) involves two difficulties, first the non differentiability of the term jrt j, secondly the positivity of the multiplier k. To circumvent these two difficulties, a rewriting of the Coulomb law is introduced. On the sliding zone, the condition: jrt j ¼
lrn

ð20Þ

is replaced by

r2t
l2 r2n ¼ 0

ð21Þ

Then relation (21) is completed by an expression meaning that uT and rT have opposite signs. Replacing the ð13Þ unknown positive multiplier k by the exponential of another variable a, we make sure the multiplier sign is positive and the correspondence between k and a is one to where rt is the tangential stress. In the adhesion zone, here one: the interval between I1 and I2 , the following conditions ut ¼
ea rt ð22Þ hold:

So on the sliding regions, the following boundary conditions hold:

un
g ¼ 0 ut ¼
ea rt

ð23Þ

The Eq. (26) is reduced to an algebraic form using the approximation of the boundary by a set of N elements:

oX ¼

r2t ¼ l2 r2n

ð27Þ

cm

m¼1

where g is the initial gap and rn ; ut anda are three unknowns. The boundary oX of the domain X is divided into four parts. Cu is the part where displacements are given, CF the part where forces are imposed, Ca and Cs are respectively the adhesion subareas and the sliding subareas of the contact zone. The formulation of the unilateral frictional contact problem whose contact state can be described by the P0 contact parameters reads:

  divrðuÞ ¼ 0 in X   u ¼ 0 on Cu   rn ¼ F on CF    un
g ¼ 0  on Ca
fIr ; r ¼ 1 to Pa g  ut ¼ 0  9  u
g ¼0 =  n    r2 ¼ l2 r2 on Cs
Mq ; Ir ; q ¼ 1 to P; r ¼ 1 to Pa  t n ;   ut ¼
ea rt9   un
g ¼ 0 =   ut ¼ 0 at fIr ; r ¼ 1 to Pa g  2  r ¼ l2 r2 ; n  t  u
g ¼ 0  n at fMq ; q ¼ 1 to Pg r ¼0 n

uj ðyÞ ffi uj ðam Þ ¼ ujm pj ðyÞ ffi pj ðam Þ ¼ pjm

y 2 cm

ð28Þ

where the point am is the middle of the segment cm whose extremities are bm and bmþ1 . So the integral equation (26) leads to a usual 2N linear equations system. Since the different contact zones are characterized by contact parameters, the discretization elements will change along with these contact parameters. Finally the integral equation (26) gives 2N non linear equations: For i ¼ 1; 2 and k ¼ 1; . . . ; N N X 2 X


Þpjm
Aijmk ðq; q

m¼1 j¼1

N X 2 X


Þujm ¼ 0 Bijmk ðq; q

m¼1 j¼1

ð29Þ

ð25Þ

ð26Þ

Fi ðXÞ ¼ 0;

Moreover, after resolution we have to control that the solution of the problem (24) satisfies:

on Cs [ Ca :

with


¼ ð

Pa Þ where the with q ¼ ðq1    qP Þ and q q1    q coefficients Aijmk and Bijmk depend on contact parameters [31]. To equation system (29) we have to add extra equations since the discretized version of the problem involves the 2N classical unknowns, the P0 contact parameters and the Ns coefficients a introduced by the Coulomb friction law (22) (Ns being the number of elements on the sliding zones). So, we add to (29) the limit contact conditions (18), the limit sliding conditions (19) and the relation (21) on the sliding zone. The relation (22) is directly taken into account in the discretized integral equations (29). The final form of the nonlinear system is then:

ð24Þ

rn  0

N [

 N 2 N P 2 PP P 

Þujm ¼ 0 A ðq; q Þp
Bijmk ðq; q ijmk jm  m¼1 j¼1 4.4 m¼1 j¼1  Non linear system form   rn ðMq Þ ¼ 0 q ¼ 1; . . . ; P The above formulation contains many boundary condi tions and the boundary integral equations method is very  ut ðIr Þ ¼ 0 r ¼ 1; . . . ; Pa well adapted to take into account this kind of conditions.  r2
l2 r2 ¼ 0 on each sliding node t n So the problem is approximated by the boundary element ð30Þ method. Under small displacements assumption and body forces This non linear system is solved by the Newton Raphson being neglected, the boundary integral equation (15), also method and the final algorithm runs with simple precision known as Somigliana identity, links displacements and reals and with double precision reals as well. tractions along the domain boundary. For any point Details about the discretized form of contact limit and x 2 oX, sliding limit conditions are given in [31]. Z Z In a more general form, this new discretized formulacij ðxÞuj ðxÞ ¼ uij ðx; yÞpj ðyÞdSy
pij ðx; yÞuj ðyÞdSy tion of unilateral contact problems with friction leads to the following nonlinear equations oX oX uij ðx; yÞ

pij ðx; yÞ

where and are respectively the displacements and stresses associated to the Kelvin fundamental solution (Brebbia et al. [32]). The term cij ðxÞ usually known as the free term depends on the local geometry at point x. At a regular boundary point x, its value is 1=2dij .

i ¼ 1; . . . ; 2N þ P0 þ Ns 0

ð31Þ 0

where X is a 2N þ P þ Ns vector composed of P contact parameters, 2N classical unknowns and Ns values of a when sliding is considered. Since only static frictional contact problems are considered, the resolution process scheme is not incremental

393

394

whereas the resolution of Eqs. (31) using the Newton Raphson method is iterative.

The gap vector in the global coordinates system between body X1 and X2 is defined by:

5 The finite element method (FEM) The popular displacement finite element method is largely used to for scientific computing in engineering. Without going into all the details, we present here just the algorithm for contact modeling. After finite element discretization in the context of small displacements, the global set of equilibrium equations of two contacting elastic bodies can be written as

X ¼ uc1
uc2 þ X0

KU ¼ F þ R



1 where S ¼ K cc is the contact flexibility matrix and Uf the gap vector resulting from each pair of contact nodes due to the application of the external loads F. Let Q be the rotation matrix between the local frame ðt1 ; t2 ; nÞ and global one ðX; Y; ZÞ. Let x and r be respectively the gap vector and the contact force vector in the local frame. Equation (38) may be written in the local frame:

ð32Þ

where K denotes the stiffness matrix. U is the displacement vector and F external known forces vector. R is the contact reactions vector. As U and R are both unknown, Eq. (32) cannot be directly solved. In the popular penalty method, the complex tree-like structure (6) or (7) is directly managed. Accordingly, K is modified by introducing contact elements and the global set of equations is solved at each iteration. The resulting numerical algorithms are not very reliable. In this paper, we propose an alternative way. The key idea is to determine iteratively the contact reactions vector R in a reduced system which only concerns the contact nodes. Then, vector U can be computed in the whole structure, using contact reactions as external loading. Consequently, the global set of equations is just solved once. Another advantage is that the stiffness matrix is not changed as opposite to the penalty method or to the Lagrange multiplier method.

5.1 Contact equations The contact phenomenon is characterized by its proper nonlinearities, due to the variation of the contact area with deformation and the irreversibility of frictional effects. In order to treat the irreversibility with the incremental formulation, the load is divided into some incremental steps and the contact equations are solved by an iterative method in each incremental step. To reduce the dimension of the problem before the iterative procedure, one can use the condensation technique. Reformulating the structure equilibrium equations under the following form :  

 Krr Krc F ur ð33Þ ¼ KTrc Kcc R uc where uc is the displacement vector of the contact node and ur for the others. The condensed degrees ur are eliminated by:

ur ¼ K
1 rr ðF
Krc uc Þ

ð34Þ

Replacing (34) into (33), one obtains:


þR
cc uc ¼ F K

ð35Þ

with


cc ¼ Kcc
KTrc K
1 K rr Krc T
1
¼
Krc Krr F F

ð36Þ

ð37Þ

where X0 is the initial gap vector. Substituting (35) into (37), the gap vector can be expressed in terms of contact forces vector and the flexibility matrix such that:

X ¼ SR þ Uf þ X0

ð38Þ

with


Uf ¼ SF

ð39Þ

x ¼ sr þ uf þ g0

ð40Þ

where

x ¼ QG;

r ¼ QR;

uf ¼ QUf ;

s ¼ QT SQ :

g0 ¼ QX0 ;

ð41Þ

5.2 Incremental form and iterative solution procedure Because of the irreversible characteristics of contact phenomena, the contact equation (40) must be written in incremental form. Then the previous mechanical quantities can be decomposed in increments following the scheme: xi ¼ xi
1 þ Dxi ;

ri ¼ ri
1 þ Dri

ð42Þ

As (40) must be satisfied at the beginning and the end of the step, it holds

sDri ¼ b

ð43Þ

with

b ¼ Dxi
Duif
Dgi

ð44Þ

where Duif and Dgi are known values and Dxi and Dri unknown ones. At each iteration, one calculates a new iteration of xi and modifies s and b, satisfying the local contact law (7) in different contact states: no contact, contact with sticking, contact with sliding. Then one can calculate the new values of Dri by (43). This iterative procedure continues until all the conditions of contact are satisfied. When all the steps are performed, the final contact forces r are obtained. Introducing these values into the Eq. (32), the unknown displacements U can be determined and additional quantities such as stresses and strains can be evaluated for the whole structure. The algorithm structure is described in Box 1 and implemented in a code named FER [33, 34].

6 Numerical examples Two typical numerical examples are solved for comparison by the previously described BEM and FEM methods and

395 Fig. 5. Example 1 Fig. 7. BEM normal contact stresses vs. contact angle

according to the friction coefficient. It is noted that these analyses were performed on a PC (Pentium III 733 MHz).

Fig. 6. FEM stress rxx

by ANSYS using the penalty method. The first example is a Signorini’s problem such that a deformable body around a rigid circular cylinder is stretched by a uniform traction. In this frictionless contact problem already treated by Kikuchi and Song [35], the interesting point is the angle of the contact area of the body. The second well known problem deals with the frictional contact between an elastic block and a rigid plate. The purpose of this example is to study the convergence of the different methods

6.1 Example 1. Frictionless contact between a 2D elastic bloc and a rigid cylinder A deformable body which is pivoted by a rigid circular cylinder is stretched uniformly as shown in Fig. 5. A plane strain state is assumed. Because of symmetry, only half of the body is modeled. The characteristics of this example are E ¼ 1000 N/mm
2 , m ¼ 0:3 and P0 ¼ 10 N/mm
2 . In BEM resolution, only the edge of the deformable body is meshed. Normal stresses are initialized to zero and the contact angle to 20 . Normal stresses given by rxx stress with FEM algorithm (FER) in Fig. 6 and by the contact normal stresses with BEM approach in Fig. 7 are very close. The CPU time is 4 s for the BEM and 15 s with FER. The size of the system to solve is 179 with BEM and 1970 with FEM. As shown in Fig. 8, the deformed edges and the contact zone are exactly the same with BEM and

Fig. 8. Zoomed deformed meshes by BEM and FEM

396

Fig. 9. Contact between a 2D elastic slab and a rigid plate (h ¼ 40 mm)

Fig. 11. Distribution of shear stress rxy ðl ¼ 1:0Þ

FEM approaches. The value of the contact angle is given according to the non deformed boundary.

6.2 Example 2. Frictional contact between a 2D elastic slab and a rigid plate This example is often studied for validation of computer codes. It is an elastic slab pressed and pushed on a rigid plate with a given friction coefficient as shown in Fig. 9. This example is of interest because we can have at the same time three different contact areas: non-contact area AB, sliding area BC and sticking area CD. Because of symmetry, only half of the slab is modeled. The characteristics of this example are:  E ¼ 13000 daN/mm2 ; t ¼ 0:2; friction coefficient: l ¼ 1:0  boundary conditions: u ¼ 0 on ED (symmetry)  loadings: f ¼ 5 daN/mm2 on GE, F ¼ 10 daN/mm2 on GA. The half slab is modeled by 561 nodes and 512 linear quadrilateral plane strain elements. Figure 10 shows the initial (BEM) and deformed (BEM and FEM) meshes (with

Fig. 12. Normal and tangential stresses by BEM and FEM (FER)

an amplification factor 500). The two deformed meshes overlap perfectly. Figure 11 shows the distribution of shear stress from which a concentration is observed at the zone of separation between sliding area BC and sticking area CD. Normal pressures and tangential stresses on the AD zone are given for the BEM and FEM (FER) resolutions in

Fig. 10. Initial and deformed meshes by FEM and BEM (amplification factor = 500)

397 Fig. 15. Normal displacements on AD zone

Fig. 13. Normal stresses by FEM (FER, ANSYS) and BEM

Fig. 14. Tangential displacements on AD zone

Fig. 12. Results are in a good agreement. Normal stresses are also compared with the ANSYS ones in Fig. 13. We can see that ANSYS results depend on the penalty coefficient Kn value (unknown a priori). Results are better with Kn ¼ 106 compared with Kn ¼ 105 . Figures 14 and 15 show the tangential and normal (penetration) displacements on AD zone by means of different methods. The remark about the influence of the Kn coefficient is still valid. It is worth noting that with ANSYS, no sticking area can be clearly identified. The contact nodes (if penetrated) are always assumed to be sliding on the target. All these illustrate that, concerning numerical stability and precision, the methods developed in the present paper have better performance than the penalty method implemented in ANSYS. Another advantage of proposed approaches is that no numerical parameters such as penalty factors are chosen by users. In order to study the influence of friction effects on the contact behaviors, this example is also modeled by

Table 1. Contact zones and CPU time l

Method

AB

BC

CD

UA (mm)

VA (mm)

CPU (s)

0.0

BEM FER ANSYS1 ANSYS2

0 0 0 0

40 40 40 40

0 0 0 0

0.025485 0.025846 0.025831 0.025845

0 0
0.3206E-04
0.3133E-05

1 4 16 19

0.4

BEM FER ANSYS1 ANSYS2

0.33 1.25 x x

36.01 35 x x

3.66 3.75 x x

0.019496 0.019122 0.019225 0.019122

0.1640E-04 0.1780E-04
0.2040E-05
0.1544E-06

7 4 17 20

0.8

BEM FER ANSYS1 ANSYS2

2.03 2.5 x x

23.58 22.5 x x

14.39 15 x x

0.015428 0.015574 0.016352 0.015678

0.3506E-03 0.3721E-03 0.2808E-03 0.3373E-03

6 5 19 22

1.0

BEM FER ANSYS1 ANSYS2

2.73 3.75 x x

19.36 17.5 x x

17.91 18.75 x x

0.014454 0.014597 0.015652 0.013387

0.4945E-03 0.5576E-03 0.4360E-03 0.2394E-05

6 6 20 24

1.3

BEM FER ANSYS1 ANSYS2

3.58 3.75 x x

15.16 15 x x

21.260 21.25 x x

0.013483 0.013642 0.015060 0.013548

0.6578E-03 0.8009E-03 0.5844E-03 0.7681E-03

7 10 22 26

ANSYS1 : penalty factor Kn ¼ 105 ; ANSYS2 : penalty factor Kn ¼ 106

398

Fig. 16. Evolution of contact zones vs. friction coefficient

7 Conclusions 1. Read the data: mesh, material properties, boundary The main purpose of this paper is to present a comparative conditions, contact definition . . . study of the boundary element method (BEM) and the finite 2. Compute K and F element method (FEM) for analysis of Signorini contact 3. Modify K and F for essential boundary conditions problems in elastostatics with Coulomb’s friction law. The 4. Compute s 5. Incremental load steps: BEM consists of solving a nonlinear system of equations 5.1. compute b involving displacements, contact tractions and geometrical 5.2. contact iterations: parameters on the boundary. The FEM developed solves the 5.2.1. modify s and b taking into account the contact states contact problem iteratively in a reduced linear system and i
1 5.2.2. solve linear equations Dr ¼ s b computes the displacements in the whole structure, using 5.2.3. update the gap vector and the contact force vector contact reactions as external loading. Two examples have 5.2.4. check contact states, if changed, go to 5.2. been studied in detail. The results have shown that these two 5.3. update step count, if simulation not complete, go to 5. 6. Compute the global contact force vector R from r methods are in good agreement and have better perfor7. Solve KU ¼ F þ R mance as compared to the penalty method used in ANSYS 8. Compute stresses and strains for each element and output. program. It is mainly due to the fact that the proposed methods (BEM and FEM) do not make intervene user-defined parameters, so that the contact conditions such as impenetrability and Coulomb’s friction law are satisfied considering other friction coefficients: l ¼ 0:0, 0.4, 0.8 and accurately. Further studies and software development are 1.3. Table 1 shows the corresponding results obtained by being undertaken for coupling of BEM and FEM. the two methods FEM-BEM and by ANSYS. UA and VA denote respectively the horizontal and vertical displaceReferences 1. Chan SH, Tuba IS (1971) A finite element method for contact ments of point A. As shown in Table 1, the CPU time problems of solid bodies. Int. J. Mech. Sci., 13: 615–639 increases with the friction coefficient in the case of FEM. 2. Tsuta T, Yamaji S (1973) Finite element analysis of contact This tendency is not observed using BEM. problem, Theory and practice in finite element structural As discussed by Klarbring et al. [36, 37] and Martins analysis. Proc. 1973 Tokyo seminar on element analysism, et al. [38], there exists a friction limit to guarantee the pp. 177–194 existence and uniqueness of solutions to quasistatic 3. Fredriksson B (1976) Finite element solution of surface nonfrictional contact problems. In the present case, we linearities in structural mechanics with special emphasis to contact and fracture mechanics problems. Comput. Struct. 6: note that the convergence of the FEM depends effec281–290 tively on the friction coefficient l and the numerical scheme diverges for l ¼ 1:4 or higher. However, in the 4. Kikuchi N, Oden JT (1984) Contact problems in elastostatics. In: Oden JT, Carey GF (eds) Finite Elements, Vol. 5, Prenticecase of BEM, the convergence of solution seems to not Hall, Englewood Cliffs, NJ be affected by the friction coefficient as shown in 5. Francavilla A, Zienkiewicz OC (1975) A note on numerical Fig. 16. For the moment, we have no explication to this computation of elastic contact problems. Int. J. Numer. Meth. difference. Eng. 9: 913–924 Box 1. Algorithm structure for solving the contact problem

6. Sachdeva TD, Ramarkrishnan CV (1981) A finite element solution for the two-dimensional elastic contact problems with friction. Int. J. Numer. Meth. Eng. 17: 1257–1271 7. Nguyen DH, De Saxce´, G (1980) Frictionless contact of elastic bodies by finite element method and mathematical programming technique. Comput. Struct. 11: 55–67 8. Klarbring A, Bjo¨rkman G (1988) A mathematical programming approach to contact problems with friction and varying contact surface. Comput. Struct. 30: 1185–1198 9. Zhong WX, Sun SM (1989) A parametric quadratic programming approach to elastic contact problems with friction. Comput. Struct. 32: 37–43 10. Bathe KJ, Chaudhary A (1985) A solution method for planar and axisymmetric contact problems. Int. J. Numer. Meth. Eng. 21: 65–88 11. Nour-omid B, Wriggers P (1986) A two-level iteration method for solution of contact problems. Comput. Meth. Appl. Mech. Eng. 54: 131–144 12. Jean M, Touzot G (1988) Implementation of unilateral contact and dry friction in computer codes dealing with large deformation problems. J. Theor. Appl. Mech. 7: 145–160 13. Alart P, Curnier A (1991) A mixed formulation for frictional contact problems prone to Newton like solution methods. Comput. Meth. Appl. Mech. Eng. 92: 353–375 14. Simo JC, Laursen TA (1992) An augmented Lagrangian treatment of contact problems involving friction. Comput. Struct. 42: 97–116 15. Klarbring A (1992) Mathematical programming and augmented Lagrangian methods for frictional contact problems. In Curnier A (ed.), Proc. Contact Mechanics Int. Symp., PPUR, pp. 369–390 16. Heegaard JH, Curnier A (1993) An augmented Lagrangian method for discrete large slip problems. Int. J. Numer. Meth. Eng. 36: 569–593 17. De Saxce´ G, Feng ZQ (1991) New inequality and functional for contact with friction: The implicit standard material approach. Mech. Struct. Mach. 19: 301–325 18. Zhong ZH (1993) Finite Element Procedures for Contactimpact problems. Oxford University Press, Oxford 19. Wriggers P (2002) Computational Contact Mechanics. John Wiley & Sons 20. Andersson T, Fredriksson B, Allan Persson BG (1980) New developments in boundary element methods. In: Brebbia CA (ed) The Boundary Element Method Applied to Two Dimensional Contact Problems. CML Publications, pp. 247–263 21. Paris F, Garrido JA (1989) An incremental procedure for friction contact problems with the boundary element method. Eng. Anal. Bound. Elem. 6(4): 202–213

22. Paris F, Garrido JA, Foces A (1994) An incremental procedure for three-dimensional contact problems with friction. Comput. Struct. 50: 201–215 23. Takahashi S, Brebbia CA (1992) A boundary element flexibility approach for solving contact problems with or without friction. Eng. Anal. Bound. Elem. 4: 24–30 24. Kwak BM, Lee SS (1988) A complementarity problem formulation for two dimensional frictional contact problems. Comput. Struct. 28(4): 469–480 25. Polizzotto C, Zito M (1998) B.I.E.M.-based variational principles for elastoplasticity with unilateral contact boundary conditions. Eng. Anal. Bound. Elem. 21(4): 329–338 26. Wriggers P, Van TV, Stein E (1990) Finite formulation of large deformation contact-impact problem with friction Comput. Struct. 37: 319–331 27. ANSYS Theory Reference (1999) Release 5.6 28. Cimetie`re A, Texier T (1995) A new numerical approach to obstacle problems. Proceedings of the Sd Contact Mechanics International Symposium, Plenum Press, New York, Londres, pp. 259–262 29. Cimetie`re A, Me´chain C (1996) Une nouvelle me´thode nume´rique pour les proble`mes de Signorini. C. R. Acad. Sci. 323: se´rie IIb, 455–461 30. Me´chain-Renaud C, Cimetie`re A (2000) Influence of a curvature discontinuity on the pressure distribution in two dimensional frictionless contact problems. Int. J. Eng. Sci. 38: 197–205 31. Me´chain-Renaud C (1998) Une nouvelle strate´gie nume´rique pour les proble`mes de contact unilate´ral – Influence de discontinuite´ de courbure. The`se, Universite´ de Poitiers 32. Brebbia CA, Telles JCF, Wrobel LC (1984) Boundary Element Technique – Theory and Applications in Engineering. Springer-Verlag, Berlin, Heidelberg, New York, Tokyo 33. Feng ZQ, Cros JM (2002) FER/SubDomain: an integrated environment for finite element analysis using object-oriented approach. Math. Modelling Numer. Anal. 36(5): 773–781 34. http://gmfe16.cemif.univ-evry.fr:8080/ feng 35. Kikuchi N, Song YJ (1980) Contact problems involving forces and moments for incompressible linearly elastic materials. Int. J. Eng. Sci. 18: 357–377 36. Klarbring A (1991) Examples of non-uniqueness and nonexistence of solutions to quasistatic contact problems with friction. Ingenieur Archiv 60: 529–541 37. Stro¨mberg N, Johansson L, Klarbring A (1996) Derivation and analysis of a generalized standard model for contact, friction and wear. Int. J. Solids Struct. 33: 1817–1836 38. Martins JAC, Manuel DP, Monteiro M, Gastaldi F (1994) On an example of non-existence of solution to a quasistatic frictional contact problem. Eur. J. Mech. A/Solids 13: 113–133

399