Mixed finite element method for contact problems

compurers Printed

& Srrucrures

in Great

Vol.

18, No.

5, pp. 937-945,

1984 C

Britain.

MIXED

FINITE ELEMENT METHOD CONTACT PROBLEMS

004s7949p 1984

$3.00 + .oo

Pergamon PressLtd.

FOR

WIESKAWOSTACHOWICZ Technical University of GdaIisk, 1l/12 Majakowskiego Street, 80-952 Gdansk, Poland (Received

23 March

1983; received for publication 10 May 1983)

Abstract-In the paper the model of finite elements for elastic contact problems was used. Real structures are modeled by finite elements and rigid finite elements. We calculate stresses in the surface of two substructures using Coulomb model of friction. The method given here is an iterative procedure which is planed to incorporate this technique in the system allowing for incremental elastic solution. The computer program is adapted to solving spatial problems.

1. INTRODUCTION In the paper the model of finite eiements for elastic contact problems was used. The method of calculation of pressures in the surface of two substructures is shown. We assume that a real structure is modeled by finite elements and rigid finite elements[l]. A contact problem for elastic bodies modeled by finite element method among other in the paper of Parsons and Wilson[2], is presented. Chan and Tuba[3], and Ohte[4], used this method for solving plane contact problems. In the paper of Schafer[S], the special contact element is presented. This element in the model of friction was used. The discrete model of real system which consists of finite elements and rigid finite elements makes possibility analysis in the case of great difference of stiffness among these elements. The mixed finite element method is better in these cases from point of view of accuracy of calculation[l]. We assume linear stiffness of structures and we calculate stresses in the surface of two substructures using Coulomb model of friction. We must consider three cases which may appear in this problem: adhesion, slide and broken contact. The method given here is an interative procedure which is planed to incorporate this technique in the system allowing for incremental elastic solution. The computer program is adapted to solving spatial problems. 2. DESCRIPTIONOF THE PROBLEM Consider a discrete model of elastic bodies when a possibility for contact exists (Fig. 1). The discrete model consists of rigid finite elements (RFE) and finite elements (FE). Rigid finite elements are connetted by elastic elements, which are presented in the paper[l]. Finite elements are connected each other, as we know, in their nodes. We assume a possibility of connection among rigid finite elements and finite elements. Special type of elements, GAP, represents two surfaces usually on different substructures, which may or may not maintain in contact. Consider two nodes p and r at surfaces of two other elements (RFE or FE). If a possibility of separation or sliding exists we assume that nodes p and r are

connected by an element GAP (Fig. 2). For each GAP three various contact conditions are expected: adhesion, slide and broken contact. In the case of adhesion or slide an element GAP are independently deformated in directions 2, S and 7 (Fig. 3), so any GAP is represented by three springs which connect above mentioned points p and r. The axes of these springs are determined by ii, i and i: Next we assume linear properties of these springs. The vector {Au} defines displacements of springs {Au} = co1 [Au!, Aa,, Au,],

(2.1)

Fig. 1. A discrete model of elastic bodies when a possibility for contact exists.

RFE RFE P

0 &

0 n

Pig. 2. An interface element (GAP) which connects two other elements (RFE or FE).

W.

OSTACHOWICZ Fn

AUn

Fig. 3. Local coordinate orientation.

where Au, and Au, describe tangential displacements in iand f direction, Au, in fi direction. Describe by k,, k, and k, coefficients of stiffness of these springs in tangential and normal direction, respectively. GAP may resist shear in tand S direction and compressive load only in the normal direction K 2.1 Normal forces Consider a contact problem at the points p and r (Fig. 2). When the normal force R is negative (compressive) the normal displacements remain in contact. In this case possibility of adhesion or slide appears. When the normal force fi is positive the contact is broken and no force is transmitted. The deformations of a spring in fi direction denote following equation Au, = u,,~- u,, + NT,

Fig. 4(b). Load-deflection

curve.

Positive

GAP (opening).

I

Fn

Fig. 4(c). Load-deflection

curve. Negative

GAP

(2.2)

where unp and u,, denote displacements of points p and I in n direction, NT is a GAP opening (if NT > 0) or preload (if NT < 0). The load-deflection curve for k,, is then represented by Fig. 4. 2.2 Tangential forces Tangential forces in GAP, T and $ are defined only for m < 0. When the absolute values of forces ?= and 3 are less than plfll (where p-the coefficient of friction) GAP does not slide. It is case of adhesion. The tangential forces deform springs in ?or f direction. The displacements of springs in the tangent directions are Au, = u,~ - u,~- T SLIDE, Au, = z+, - u,, - S SLIDE,

(2.3)

where u,~, uspand u,,, u,, denote displacements of the points p and r in S and i direction, respectively. T SLIDE and S SLIDE are the accumulated amounts of sliding in iand S direction, respectively. The load-deflection curve for k, is shown in Fig. 5(a). In this case the accumulated amount of sliding, T SLIDE, equals zero. This same curve for T SLIDE # 0 is pictured in Fig. 5(b). The accumulated values of T SLIDE, we calculate as follows

T SLID& = (utp- utr)i- (PINi(/kJ,

4kJn

Fn Fig. 4.(a).

Load-deflection

curve for k,.

(2.4)

where i denotes a number of step of iteration. The analogous equations describe displacements of GAP in S direction.

F?l

-dun

(preload).

Fig. 5(a). Load-deflection

curve for k,

939

Mixed finite element method for contact problems

local coordinates. Consider the points p and r (Fig. 6). The system of local (c S, 2) and global (X,, X,, X,) coordinates is shown there. We assume that directions bf axes i; 5, fi are the same at the point p and r. In this case matrices of directional coefficients are the same for each point @ and I). The equations of equilibrium of GAP in local coordinates have the following form

(&I = Kl(Au&

(3.2)

where

rs;?

Fig. 5(b). Accumulated slip USLIDE.

[ZJ tptuspyunp iitn riWIuw]I

(A@$= The slide S SLIDE, has the following form S SLIDE, = (c - u,Ji - (p/N&k&

k, 0

(2.5)

For each structural interface element, stiffnesses and forces are updated based upon current displacement values in order to predict new displacements. The iteration is continued until GAP status or sliding force remains unchanged (within previously done some tolerance) for two succesive iterations.

[KJ=

_; -

0’ 0

0 k,

0 0

-k, 0

0 -k,

0 0

;

t

;

;

-k, 0

0 -k,

0’ 0

k, 0

-2 i” _

The eqn (3.2) we may write in global coordinates. From Fig. 6 we obtain (3.3)

3. STIFFNESS MATRIX

The global stiffness matrix we obtain as a sum of strain energy of elastic elements which connect rigid finite elements, strain energy of finite elements and GAPS. The several finite elements (RFE, FE and GAP) are described in local coordinates of these elements. In the case of rigid finite elements the local coordinates are central axes of inertia. The stiffness matrix in this case we form as in the paper [6]. Local stiffness matrices of finite elements we form according to the paper of Bathe and Wilson[7]. The discrete model of structures we transform to the global system of coordinates X,, X,, X,. In these coordinates we describe the generalized displacements of rigid finite elements, displacements of nodes of finite elements and displacements of GAP’s nodes. We bilaterally multiply the local stiffness matrix of directional coefficients. In the case of rigid finite elements method matrix of directional coefficients has the following form

I e

rl I

I%]= --

b2

eR2

er1310 er2,I 0

0

0

0

0

0

0

__!%“__ $i Sym. I

.

4i2 ai, e r22 e,,

I

where

-eilI oil2 e 52

[ej=

__

_

-

413

I0

0

0

en3

I 0

0

0

0

0

h3 , 0 -,T,-e._ey

,

Ill SW.

1

rl2 ea2

I

I13 h3 ei33 j

where

eilg= cos (iTX/J), eiza= cos (5, q,

eisg =cos(6 .&I>, /I = 1,2,3.

, (3.1)

e r33

where: B,, = cos (X,e, X6), t(, /r’= 1, 2, 3; X,, are the axes of local coordinates; and X, are the axes of global coordinates. In the case of finite elements method we bilaterally multiply the local stiffness matrix by the matrix of directional coefhcients, as in the paper[7]. The stiffness matrix of GAP we form on the base of equations of equilibrium which are described in

Fig. 6. Global-local coordinate orientation.

W. OSTACHOWICZ

940

The vector of’nodal forces in global coordinates has the following form {FfJ = ]V(FLj>

(3.4)

where {F,}r = Kp, Fzp,F3prF,,, F2,, F3rl. The vector {FL}is shown in (3.2). Multiplying left side of the eqn (3.2) by the matrix [QJr we obtain [&I’{&} = N~K~{A~L}~

Elements of the vector {q,} define the generalized coordinates of rigid finite elements, XLs and b,& = 1,2, 3) are axes of local coordinates of RFE and coordinates of any point (p or r), respectively (Fig. 7). If we consider a contact problem among a rigid finite element and a finite element we assume that one of two points (p and r) is in the node of the finite element. Displacements of this point are the same as displacements of the node. The global matrix of stiffness we form as a sum of local matrices of stiffness, which previously were transformed to the global coordinates. This procedure is well known (see for instance, [ 1,7]).

and next according to (3.4) {&} =

FVKl{~~~}.

(3.5)

If we put (3.3) to the eqn (3.5) we obtain {F,} = kl’[&l[&l{~~& or {FG}= Kcl{A~c). where [K,J is the global matrix of stiffness of each GAP. The matrix [K,] has the following form &I =

Ml’l&lM.

If any point of GAP (p or r) is on the surface of a rigid finite element the vector of displacements of this point we describe as a function of generalized coordinates of this RFE. The vector of displacements has the following form {A+} = [e,,]]B,,]{qJ,

Xr2

(3.6)

Fig. 7. The axes of local coordinates of RFE. The coordinates of any point (p or r).

(3.7) 4. SLIDE MODEL OF GAP If sliding of any GAP appears we introduce special

where

model of stiffness of this element. In this case we modify the load-deflection curves (see Figs. 4 and 5). It is necessary because of accuracy of numerical calculations. If any GAP is open we may not assume zero stiffness of this GAP-the reason is instability of iterational procedure. When a normal force is positive we assume a little stiffness of spring (relatively to stiffness k,). Similar, if tangential forces are positive we assume a little stiffness of tangential springs (relatively to stiffness k, and k,). The load-deflection curves for these springs in Figs. 8 and 9 are pictured. In this case normal and

{4JT = [CC,,. . . . . . ?$61, where oris,, = cos (xLB,r), &$,Z= cos (XL&% S), e,,,, = cos (xLa, n),

-R

for Fig. 8. Modified load-deflection curve for k,.

941

Mixed finite element method for contact problems 6. NUMERICAL

T

_---

-

-pl~I

q&z?

Fig. 9. Modified load-deflection curve for k,.

tangentiai forces which load any GAP have the following form N = k, . Au,,, for

Au, < 0, (4.1)

N = RIGID.

k, 1Au,,

for

Ay 2 0,

and T =O,

for

T = k, . Au,, T=sign

Au,, 20 for

jk,. Au,/< PIN],

,ulNI+RIGID.k;

(4.2)

( IAUlI-- F)],

I where sign = + 1,

for

Au, 2 0,

sign=-1,

for

Au,
and S = 0,

for

.

S = k, Au,, S =sign

Au,>O, for

jk,. Au,~/< PINI,

p[NI+RIGID+k;

( IAa, I- F)]>

where sign = -t 1, for

Au, 2 0,

sign = - 1,

Au, < 0.

for

5. CONVERGENCE OF MEASURES Convergence of the solution technique may be measured by monitoring the status of each element and change occuring between successive iterations and by measuring the percent change in friction force plN[ between successive iterations. The status I FLAG of an interface element may be I FLAG = 1, for rigid contact, I FLAG = 2, for sliding contact, I FLAG = 0, for free body. Iteration is continued until the status of all elements does not change between two successive iterations. For structures having complex loading histories or much relative motion an incremental application of load may be necessary in order to obtain convergence.

7. EXAMPLE

(4.3)

[

PROCEDURE

The scheme of computational procedure in Fig. 10 is shown. The program consists of 9 subroutines. It was used to analysis of structures modeled by a rigid finite element and finite elements (plates). The subroutines in Fortranlanguage were written In the main program (PLARFE) all problem parameters such as convergence tolerance number of load steps and maximum allowable iterations per step must be supplied. Atso each interface substructure must be assigned as coefficient of friction @), initial conditions for each interface element (IFLAG), initial GAP signs, interface tangent and normal stiffnesses (k,, k, and k,) and a rigid body stiffness multipler (RIGID). The local matrices of stiffness of finite elements we form in the subroutine PLATE. The same matrices for rigid finite elements we form in the subroutine STIRFE. The subroutine PLATE was written similar to procedure which in[8] is presented. The subroutine ASSEMB forms the global matrix of stiffness. As a next step of calculations we solve linear system of algebraic equations. We use Gaussian elimination procedure to this problem (subroutine SOLVG). Having results of subroutine SOLVG we use the subroutine LOGIC which determines displacements and forces in GAPS. The subroutine LOGIC changes stiffness of GAPS according to functions presented in Figs. 8 and 9. It should be noted that the gap closure/finite friction problem is a nonlinear one and that solution convergence and unusual results may be sensitive to many paramters such as boundary conditions (displacements and loads), load slip size coefficient of friction and interface tangent and normal stiffnesses. The analysis procedure for this class of problems may involve several analysis to assure invariance of results.

Determine displacements of nodes of finite elements, rigid finite element and GAPS for the structure which in Fig. 11 is shown. Also, determine normal and tangential forces at A 1, . . . I A 6 points and their status (IFLAG parameter). Dimensions of finite elements and rigid finite element in Fig. 11 are shown. Other data: --length of FE ---width of FE -thickness of FE -length of RFE -width of RFE -height of RFE

a=O.l m b=O.l m h=O.OlSm a, = 0.2 m b, =O.l m h, = 0.1 m

--Young’s

E= 0.21 x 1012$

-Poisson’s

modulus ratio

v = 0.3.

The finite elements have the same dimensions and elastic properties. We used rectangular elements which have 4 nodes (3 degrees of freedom per node). These elements in[9] are described. The rigid finite element is loaded by two forces, Er in tangential and pN in normat direction, respectively. We assumed constant value of normal force,

W. OSTACHOWICZ

942

fTnpuL,

1

data

I form

the

lncal

matrices

of

stiffness

for

FE

of

stiffness

for

FIFE

I Form

the

local

matrices

1

0 Form

the

global

matrices

of

stiffness

T Introduce

Solve

the

fDr

stiffness

interface

of

GHD

substructure

displacements

Llraate

interface

elements

status

forces

Uutcut

edit

interface

and

continuum

substructures

0 5TOP

Fig. IO. Computational

procedure for problem solution.

FN = 2.10’N. The values of tangential Table 3 are shown. The elastic properties

forces F, in of GAPS are

following: N k N = 0.2 lo’“m' IQ= 0.2. 105 m The coefficient

of friction of all GAPS we assumed

as

/I =O.l. Fig. 11. A finite element representation of a rigid body and plate.

The results of calculations presented. The displacements

in Tables l-3 are of nodes of finite ele-

943

Mixed finite element method for contact problems Table 1. The displacements of nodes of finite elements for Fr = 18000 N I

I I

I lSt iteration

jth

iteration

_______:_______+_______+______i__ ________i__________i_________-_ i ; cm ‘I md i red i cm I rsd : red

I I

I I

1

i 0,000

; 0.0565

, I

I I I

; -0,054a

; 0,000

I

I I

0,0501

;

I I

f -0,0854

2

, I 0,483

I 0,029l

, -0,0303 I

i 0,430

i 0,0260

, I -0,0594

3

I 0,480

I ; 0.0299

I -0.0337

I 0,769

i 0,0221

1 -0,0584

4

I 0,754

i 0.0175

; -0,0128

i 0,989

i 0,0147

i -0.0419

5

' 0,483 I ; I 0,755 I 0,000

j-o,0291

i -0,0303 1 -0,0128 : -0,0548

I 0,430 I ! 0,989 1 I 0,000

i -0,026o

;-0,017' I i-o,0565 j-O,0299

9

I 0,480 I ; 0,636

I 0,OlOl

I ' -0,0337 I 1 0,0013

I ; 0,769 I : 1,060

i -0,0594 I i -0.0419 1 ; -0p354

10

I ; 0,741

I ; 0,0075

I ;

0,003O

! ; 1,170

11

I I 0,741

j-o.0075

i 1,170 I ' 1,060 I

6 7 8

i -0,0147 I -0,050l I

I ; OS0079

1 -0,0584 I 0,0018 I

; 0,0066 I I -0,0006

I : ;

0,0039

i -0,0079

i

0,0018

; I 0,0176

i I

0,0584

: I I

0,0455

; -0,022o

12

! 0,636

I-0,OlOl I

i 0,003o I 1 0,0013

13

1 0,463

i 0,0265

I i

0.0336

i 0,746

14

i 0,706

i 0.0155

i

0,0156

; 0,924

i

15 16

a-o,0155 I j-o.0265

I :

O,a56, 0,0336

I 0,924 i 0,746

17

I , 0,706 : 0,463 I 1 0,000

I 0,0497 ;

I ;

0,0523

I 0,000 i

, I -0,012l I I -0,0176 I 0,0411 ;

i I 1

18

I 0,425

I 0,0256

I I

0.0308

! 0,354

19 20

I 0,425 I 1 0,000

~-0.0256 I f :-0,049

: I i

0,030~:

0,354 I ; 0,000

,

I

I

I 0.0214 I : -0,0214 I 1 -0,041l I

I I : I ; I

0.0523

,

0,0121

0,0039

O&55

0.0584 0,082o 0,0601 0,0601 0,082o

Table 2. Displacements of the rigid finite element

I

f

I i I Iteration

I

Displacements

I I I

I I

lSt

I ,

o,Ol40

i i :______________*__________

I

I I

gth

: , I

0,0188

ments for FT = 18 . 103N in Table 1 are shown. Numbers of these nodes in Fig. 12(a) are shown. Results of the first step of the iteration in the first column and in Fig. 13 are presented. Results of the fifth step of the iteration in the second column and in Fig. 14 are presented.

I ! 0,735O

i -t---------[ 0,9590 I

I

I !

O&O246

I !

i i j____________’

I

i

0,00326

; I

Displacements of the rigid finite element in Table 2 are shown (for first and fifth step of iteration). The displacements qi(i = 1,2,3) are assumed according to Fig. 12(b). Results for 5 variants of calculations in Table 3 are shown. Numbers of GAPS in Fig. 11 are signed. The

W. OSTACHOWICZ

944

Table 3. The normal and tangential forces of GAPS for five variants of calculations 4 I

I

Number

II

i

II

I

I

I

1

I

I I of 1 A2 j A3 1 A4 ! A5 1 A6 f I GAP I Al I : I I :______,_______!______,__‘-________C___________________________________~ = 15000 N I:_____‘__‘___t___‘___‘--____‘-” FT------~~~~~-~~~~--~--~~-_~~____~ ! I ; IN 1 0 1 51500 i 51500 1 0 i 48500 I 48500 1 k______’ “_______j________ L______&______;________: :----+------1 i T I Ii 0 i 0 i jOO0 I 3000 i i 3000 1 3000 ; ____,,L_____f_______

+_______+_______

y_____

+_______

:---------:

0 I;110 1 ; 1 I I I r__________-__-__-__-__________-'___-______________________________~ I = 14000 N I % l _-_-_-__-__-__-__+ -------~----_-__+_______r--_____r----__~~~~_~~~~~_~~~~ I iijNI 1 0 \ 48250 i 40250 51750 i 51750 ; 0

IFLAG

: ------e..----A i T’ I N

_______:________t_______ 3500 ; 0 i 3500 ;

+_____ :_____+_______‘________ t-------IFLAG 0

:______ + _______ +_______:

i 0 j 3500 ; 3500 L______*----_--_~__--_-__~

I IO I I I I I ___‘____J __ - __‘_ __ _I________1_______*___1___,____l____l r-----------= 16000 N I FT p--_________ -r_______~________t______-r______T-_____-_~________~

: 1 I

I

i I ,

I

;j f N : 52000 : 52000 1 0 : 0 i 46000 : 48000 I L ______:_____A_______:________&______&_____A _______:________: I i :Nj 4000 i 4000 i 0 / 0 i 4000 : 4000 i L______:_____*_______L______+______ *______+_______:________~ t IFLAG tI I 0 I I i 0 1 1 i 1 1 I i :____________*_______~________~______*______~_______~________~ I

FT= 18000 N

:r’_‘__‘r_‘“‘_T__“-__,__-~~__‘-I-

I

i ! N 1 52250 : 52250 ! 0 I 0 I 47750 i 47750 I L__=__&____ &______ &______ &_____ _c______ +_______ I :--------: 0 I 0 I 4500 ’ 4500 I i T 4500 1 4500 I I’ N : _______:________&__ +_____+_____I ____~______~_______&______~ 1Ii~0~0~1~1~ 1 IFLAG I L____,__,__,_l_,_____~____~___---lI I N PT - 20000 &Z___r_____-r_______,________r--______ r______,________7________j : 0 : N : N : 52500 ! 52500 ! 0 I 47500 I 47500 I :__Z___:_____ +_______ &______- +______ +_____ +_______ +_______; 0 IO IO IO IO IO I : N t i T ________&______A _____%_______:________: L ______:_____I_______: I 2; 2 IO 10: 2; 2 1 I IFLAG 1 i

I

Fig. 12(a). Numbers of nodes of the finite elements.

illi 4

-,_-_-.

&-+c__

/’ (?f I

/’

43

--

9

t

Fig. 12(b). Generalized coordinates of a rigid finite element.

Fig. 13. Displacements of the nodes for FT= 18000 N. Results of the first step of iteration.

Mixed finite element method for contact problems

2. 3. 4. Fig. 14. Displacements of the nodes for Fr= 18000N. Results of the fifth step of iteration.

5.

6.

normal force and tangential force of GAP we signed by fl and T, respectively. Status of GAP (parameter IFLAG) is defined as in Section 5. The results which in Table 3 are presented we obtained after 5 steps of iterations.

7.

8. 9.

945

REFERENCES J. Kruszewski, The Rigid Finite Elements Method (in Polish). Arkady, Warsaw (1975). S. Parsons and E. A. Wilson, Finite element analysis of elastic contact problem using differential displacements. ht. J. Num. Meth. Engng 2, 387-395 (1970). S. K. Chan and T. S. Tuba, A finite element method for contact problems of solid bodies-1 and II. ht. J. Mech. Sci. 13, 61.5639 (1971). S. Ohte, Finite element analysis of elastic contact problems. Bull. JSME 16, 797-804 (1973). H. Schafer, A contribution to the solution of contact problems with the aid of bound elements. Comput. Meth. Appl. Mech. Engng 6, 335-354 (1975). J. Kruszewski, HESAS-System of structural computation based on methods of finite elements. Comput. Mech. (in Polish) 1, 81-100 (1978). K. J. Bathe and E. L. Wilson, Numerical Methods in Finite Element Analysis. Prentice-Hall, Englewood Cliffs, New Jersey (1976). C. S. Desai, Finite Element Methodr. Prentice-Hall, Englewood Cliffs, New Jersey (1979). 0. C. Zienkiewicz, The Finite Element Method in Engineering Science. McGraw-Hill, London (1971).